BTL: fix warnings and extend to 5k matrices, update GotoBlas to OpenBlas, etc.

This commit is contained in:
Gael Guennebaud 2014-03-31 10:58:30 +02:00
parent 8d0441052e
commit f603823ef3
19 changed files with 104 additions and 136 deletions

View File

@ -129,13 +129,25 @@ int main(int argc, char ** argv)
int tries = 2; // number of tries, we keep the best
int s = 2048;
int m = s;
int n = s;
int p = s;
int cache_size = -1;
bool need_help = false;
for (int i=1; i<argc; ++i)
{
if(argv[i][0]=='s')
{
s = atoi(argv[i]+1);
m = n = p = s;
}
else if(argv[i][0]=='m')
{
m = atoi(argv[++i]);
n = atoi(argv[++i]);
p = atoi(argv[++i]);
}
else if(argv[i][0]=='c')
cache_size = atoi(argv[i]+1);
else if(argv[i][0]=='t')
@ -155,9 +167,7 @@ int main(int argc, char ** argv)
if(cache_size>0)
setCpuCacheSizes(cache_size,96*cache_size);
int m = s;
int n = s;
int p = s;
A a(m,p); a.setRandom();
B b(p,n); b.setRandom();
C c(m,n); c.setOnes();

View File

@ -104,6 +104,7 @@ add_subdirectory(libs/mtl4)
add_subdirectory(libs/blitz)
add_subdirectory(libs/tvmet)
add_subdirectory(libs/STL)
add_subdirectory(libs/blaze)
add_subdirectory(data)

View File

@ -33,7 +33,7 @@ class Action_axpby {
public :
// Ctor
Action_axpby( int size ):_size(size),_alpha(0.5),_beta(0.95)
Action_axpby( int size ):_alpha(0.5),_beta(0.95),_size(size)
{
MESSAGE("Action_axpby Ctor");

View File

@ -35,7 +35,7 @@ public :
// Ctor
Action_axpy( int size ):_size(size),_coef(1.0)
Action_axpy( int size ):_coef(1.0),_size(size)
{
MESSAGE("Action_axpy Ctor");

View File

@ -4,10 +4,7 @@ if (ATLAS_LIBRARIES)
endif (ATLAS_LIBRARIES)
find_file(ATLAS_LIB libatlas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_LIB atlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_CBLAS libcblas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_CBLAS cblas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_LIB satlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_LAPACK liblapack_atlas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_LAPACK lapack_atlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
@ -22,14 +19,14 @@ find_library(ATLAS_F77BLAS f77blas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
if(ATLAS_LIB AND ATLAS_CBLAS AND ATLAS_LAPACK AND ATLAS_F77BLAS)
set(ATLAS_LIBRARIES ${ATLAS_LAPACK} ${ATLAS_CBLAS} ${ATLAS_F77BLAS} ${ATLAS_LIB})
set(ATLAS_LIBRARIES ${ATLAS_LAPACK} ${ATLAS_LIB})
# search the default lapack lib link to it
find_file(ATLAS_REFERENCE_LAPACK liblapack.so.3 PATHS /usr/lib /usr/lib64)
find_library(ATLAS_REFERENCE_LAPACK NAMES lapack)
if(ATLAS_REFERENCE_LAPACK)
set(ATLAS_LIBRARIES ${ATLAS_LIBRARIES} ${ATLAS_REFERENCE_LAPACK})
endif()
# if(ATLAS_REFERENCE_LAPACK)
# set(ATLAS_LIBRARIES ${ATLAS_LIBRARIES} ${ATLAS_REFERENCE_LAPACK})
# endif()
endif(ATLAS_LIB AND ATLAS_CBLAS AND ATLAS_LAPACK AND ATLAS_F77BLAS)

View File

@ -1,15 +0,0 @@
if (GOTO_LIBRARIES)
set(GOTO_FIND_QUIETLY TRUE)
endif (GOTO_LIBRARIES)
find_library(GOTO_LIBRARIES goto PATHS $ENV{GOTODIR} ${LIB_INSTALL_DIR})
if(GOTO_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
set(GOTO_LIBRARIES ${GOTO_LIBRARIES} "-lpthread -lgfortran")
endif()
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(GOTO DEFAULT_MSG GOTO_LIBRARIES)
mark_as_advanced(GOTO_LIBRARIES)

View File

@ -1,25 +0,0 @@
if (GOTO2_LIBRARIES)
set(GOTO2_FIND_QUIETLY TRUE)
endif (GOTO2_LIBRARIES)
#
# find_path(GOTO_INCLUDES
# NAMES
# cblas.h
# PATHS
# $ENV{GOTODIR}/include
# ${INCLUDE_INSTALL_DIR}
# )
find_file(GOTO2_LIBRARIES libgoto2.so PATHS /usr/lib $ENV{GOTO2DIR} ${LIB_INSTALL_DIR})
find_library(GOTO2_LIBRARIES goto2 PATHS $ENV{GOTO2DIR} ${LIB_INSTALL_DIR})
if(GOTO2_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
set(GOTO2_LIBRARIES ${GOTO2_LIBRARIES} "-lpthread -lgfortran")
endif()
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(GOTO2 DEFAULT_MSG
GOTO2_LIBRARIES)
mark_as_advanced(GOTO2_LIBRARIES)

View File

@ -0,0 +1,17 @@
if (OPENBLAS_LIBRARIES)
set(OPENBLAS_FIND_QUIETLY TRUE)
endif (OPENBLAS_LIBRARIES)
find_file(OPENBLAS_LIBRARIES libopenblas.so PATHS /usr/lib $ENV{OPENBLASDIR} ${LIB_INSTALL_DIR})
find_library(OPENBLAS_LIBRARIES openblas PATHS $ENV{OPENBLASDIR} ${LIB_INSTALL_DIR})
if(OPENBLAS_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
set(OPENBLAS_LIBRARIES ${OPENBLAS_LIBRARIES} "-lpthread -lgfortran")
endif()
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(OPENBLAS DEFAULT_MSG
OPENBLAS_LIBRARIES)
mark_as_advanced(OPENBLAS_LIBRARIES)

View File

@ -1,19 +1,19 @@
aat ; "{/*1.5 A x A^T}" ; "matrix size" ; 4:3000
ata ; "{/*1.5 A^T x A}" ; "matrix size" ; 4:3000
atv ; "{/*1.5 matrix^T x vector}" ; "matrix size" ; 4:3000
aat ; "{/*1.5 A x A^T}" ; "matrix size" ; 4:5000
ata ; "{/*1.5 A^T x A}" ; "matrix size" ; 4:5000
atv ; "{/*1.5 matrix^T x vector}" ; "matrix size" ; 4:5000
axpby ; "{/*1.5 Y = alpha X + beta Y}" ; "vector size" ; 5:1000000
axpy ; "{/*1.5 Y += alpha X}" ; "vector size" ; 5:1000000
matrix_matrix ; "{/*1.5 matrix matrix product}" ; "matrix size" ; 4:3000
matrix_vector ; "{/*1.5 matrix vector product}" ; "matrix size" ; 4:3000
trmm ; "{/*1.5 triangular matrix matrix product}" ; "matrix size" ; 4:3000
trisolve_vector ; "{/*1.5 triangular solver - vector (X = inv(L) X)}" ; "size" ; 4:3000
trisolve_matrix ; "{/*1.5 triangular solver - matrix (M = inv(L) M)}" ; "size" ; 4:3000
cholesky ; "{/*1.5 Cholesky decomposition}" ; "matrix size" ; 4:3000
complete_lu_decomp ; "{/*1.5 Complete LU decomposition}" ; "matrix size" ; 4:3000
partial_lu_decomp ; "{/*1.5 Partial LU decomposition}" ; "matrix size" ; 4:3000
tridiagonalization ; "{/*1.5 Tridiagonalization}" ; "matrix size" ; 4:3000
hessenberg ; "{/*1.5 Hessenberg decomposition}" ; "matrix size" ; 4:3000
symv ; "{/*1.5 symmetric matrix vector product}" ; "matrix size" ; 4:3000
syr2 ; "{/*1.5 symmretric rank-2 update (A += u^T v + u v^T)}" ; "matrix size" ; 4:3000
ger ; "{/*1.5 general rank-1 update (A += u v^T)}" ; "matrix size" ; 4:3000
matrix_matrix ; "{/*1.5 matrix matrix product}" ; "matrix size" ; 4:5000
matrix_vector ; "{/*1.5 matrix vector product}" ; "matrix size" ; 4:5000
trmm ; "{/*1.5 triangular matrix matrix product}" ; "matrix size" ; 4:5000
trisolve_vector ; "{/*1.5 triangular solver - vector (X = inv(L) X)}" ; "size" ; 4:5000
trisolve_matrix ; "{/*1.5 triangular solver - matrix (M = inv(L) M)}" ; "size" ; 4:5000
cholesky ; "{/*1.5 Cholesky decomposition}" ; "matrix size" ; 4:5000
complete_lu_decomp ; "{/*1.5 Complete LU decomposition}" ; "matrix size" ; 4:5000
partial_lu_decomp ; "{/*1.5 Partial LU decomposition}" ; "matrix size" ; 4:5000
tridiagonalization ; "{/*1.5 Tridiagonalization}" ; "matrix size" ; 4:5000
hessenberg ; "{/*1.5 Hessenberg decomposition}" ; "matrix size" ; 4:5000
symv ; "{/*1.5 symmetric matrix vector product}" ; "matrix size" ; 4:5000
syr2 ; "{/*1.5 symmretric rank-2 update (A += u^T v + u v^T)}" ; "matrix size" ; 4:5000
ger ; "{/*1.5 general rank-1 update (A += u v^T)}" ; "matrix size" ; 4:5000
rot ; "{/*1.5 apply rotation in the plane}" ; "vector size" ; 4:1000000

View File

@ -10,7 +10,7 @@ ublas ; with lines lw 3 lt 1 lc rgbcolor "#00b7ff"
mtl4 ; with lines lw 3 lt 1 lc rgbcolor "#d18847"
blitz ; with lines lw 3 lt 1 lc rgbcolor "#ff00ff"
F77 ; with lines lw 3 lt 3 lc rgbcolor "#e6e64c"
GOTO ; with lines lw 3 lt 3 lc rgbcolor "#C05600"
GOTO2 ; with lines lw 3 lt 1 lc rgbcolor "#C05600"
OPENBLAS ; with lines lw 3 lt 1 lc rgbcolor "#C05600"
C ; with lines lw 3 lt 3 lc rgbcolor "#e6bd96"
ACML ; with lines lw 2 lt 3 lc rgbcolor "#e6e64c"
blaze ; with lines lw 3 lt 1 lc rgbcolor "#ff00ff"

View File

@ -33,7 +33,7 @@
// min matrix size for matrix vector product bench
#define MIN_MV 5
// max matrix size for matrix vector product bench
#define MAX_MV 3000
#define MAX_MV 5000
// min matrix size for matrix matrix product bench
#define MIN_MM 5
// max matrix size for matrix matrix product bench
@ -41,7 +41,7 @@
// min matrix size for LU bench
#define MIN_LU 5
// max matrix size for LU bench
#define MAX_LU 3000
#define MAX_LU 5000
// max size for tiny vector and matrix
#define TINY_MV_MAX_SIZE 16
// default nb_sample for x86 timer

View File

@ -176,7 +176,7 @@ public:
if (_config!=NULL)
{
std::vector<BtlString> config = BtlString(_config).split(" \t\n");
for (int i = 0; i<config.size(); i++)
for (unsigned int i = 0; i<config.size(); i++)
{
if (config[i].beginsWith("-a"))
{
@ -224,7 +224,7 @@ public:
return false;
BtlString name(_name);
for (int i=0; i<Instance.m_selectedActionNames.size(); ++i)
for (unsigned int i=0; i<Instance.m_selectedActionNames.size(); ++i)
if (name.contains(Instance.m_selectedActionNames[i]))
return false;

View File

@ -30,23 +30,23 @@ double simple_function(int index_i, int index_j)
return index_i+index_j;
}
double pseudo_random(int index)
double pseudo_random(int /*index*/)
{
return std::rand()/double(RAND_MAX);
}
double pseudo_random(int index_i, int index_j)
double pseudo_random(int /*index_i*/, int /*index_j*/)
{
return std::rand()/double(RAND_MAX);
}
double null_function(int index)
double null_function(int /*index*/)
{
return 0.0;
}
double null_function(int index_i, int index_j)
double null_function(int /*index_i*/, int /*index_j*/)
{
return 0.0;
}

View File

@ -23,7 +23,7 @@
#include "size_log.hh"
template<class Vector>
void size_lin_log(const int nb_point, const int size_min, const int size_max, Vector & X)
void size_lin_log(const int nb_point, const int /*size_min*/, const int size_max, Vector & X)
{
int ten=10;
int nine=9;

View File

@ -18,27 +18,14 @@ if (MKL_FOUND)
endif (MKL_FOUND)
find_package(GOTO2)
if (GOTO2_FOUND)
btl_add_bench(btl_goto2 main.cpp)
if(BUILD_btl_goto2)
target_link_libraries(btl_goto2 ${GOTO_LIBRARIES} )
set_target_properties(btl_goto2 PROPERTIES COMPILE_FLAGS "-DCBLASNAME=GOTO2")
endif(BUILD_btl_goto2)
endif (GOTO2_FOUND)
find_package(GOTO)
if (GOTO_FOUND)
if(GOTO2_FOUND)
btl_add_bench(btl_goto main.cpp OFF)
else()
btl_add_bench(btl_goto main.cpp)
endif()
if(BUILD_btl_goto)
target_link_libraries(btl_goto ${GOTO_LIBRARIES} )
set_target_properties(btl_goto PROPERTIES COMPILE_FLAGS "-DCBLASNAME=GOTO")
endif(BUILD_btl_goto)
endif (GOTO_FOUND)
find_package(OPENBLAS)
if (OPENBLAS_FOUND)
btl_add_bench(btl_openblas main.cpp)
if(BUILD_btl_openblas)
target_link_libraries(btl_openblas ${GOTO_LIBRARIES} )
set_target_properties(btl_openblas PROPERTIES COMPILE_FLAGS "-DCBLASNAME=OPENBLAS")
endif(BUILD_btl_openblas)
endif (OPENBLAS_FOUND)
find_package(ACML)
if (ACML_FOUND)

View File

@ -75,7 +75,6 @@ public :
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
int N2 = N*N;
BLAS_FUNC(copy)(&N2, X, &intone, C, &intone);
char uplo = 'L';
int info = 0;
int * ipiv = (int*)alloca(sizeof(int)*N);
BLAS_FUNC(getrf)(&N, &N, C, &N, ipiv, &info);
@ -92,7 +91,7 @@ public :
BLAS_FUNC(trsm)(&right, &lower, &notrans, &nonunit, &N, &N, &fone, L, &N, X, &N);
}
static inline void trmm(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
static inline void trmm(gene_matrix & A, gene_matrix & B, gene_matrix & /*X*/, int N){
BLAS_FUNC(trmm)(&left, &lower, &notrans,&nonunit, &N,&N,&fone,A,&N,B,&N);
}
@ -101,7 +100,6 @@ public :
static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
int N2 = N*N;
BLAS_FUNC(copy)(&N2, X, &intone, C, &intone);
char uplo = 'L';
int info = 0;
int * ipiv = (int*)alloca(sizeof(int)*N);
int * jpiv = (int*)alloca(sizeof(int)*N);
@ -134,8 +132,6 @@ public :
}
char uplo = 'U';
int info = 0;
int ilo = 1;
int ihi = N;
int bsize = 64;
int worksize = N*bsize;
SCALAR* d = new SCALAR[3*N+worksize];

View File

@ -17,12 +17,12 @@ public:
typedef real* gene_matrix;
typedef real* gene_vector;
static void free_matrix(gene_matrix & A, int N){
delete A;
static void free_matrix(gene_matrix & A, int /*N*/){
delete[] A;
}
static void free_vector(gene_vector & B){
delete B;
delete[] B;
}
static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){

View File

@ -44,9 +44,9 @@ public :
return "STL";
}
static void free_matrix(gene_matrix & A, int N){}
static void free_matrix(gene_matrix & /*A*/, int /*N*/){}
static void free_vector(gene_vector & B){}
static void free_vector(gene_vector & /*B*/){}
static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){
A = A_stl;

View File

@ -45,9 +45,9 @@ public :
return EIGEN_MAKESTRING(BTL_PREFIX);
}
static void free_matrix(gene_matrix & A, int N) {}
static void free_matrix(gene_matrix & /*A*/, int /*N*/) {}
static void free_vector(gene_vector & B) {}
static void free_vector(gene_vector & /*B*/) {}
static BTL_DONT_INLINE void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){
A.resize(A_stl[0].size(), A_stl.size());
@ -74,7 +74,7 @@ public :
}
static BTL_DONT_INLINE void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){
int N=A_stl.size();
int N=A_stl.size();
for (int j=0;j<N;j++){
A_stl[j].resize(N);
@ -84,28 +84,28 @@ public :
}
}
static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){
static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int /*N*/){
X.noalias() = A*B;
}
static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){
static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int /*N*/){
X.noalias() = A.transpose()*B.transpose();
}
// static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){
// static inline void ata_product(const gene_matrix & A, gene_matrix & X, int /*N*/){
// X.noalias() = A.transpose()*A;
// }
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int /*N*/){
X.template triangularView<Lower>().setZero();
X.template selfadjointView<Lower>().rankUpdate(A);
}
static inline void matrix_vector_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N){
static inline void matrix_vector_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int /*N*/){
X.noalias() = A*B;
}
static inline void symv(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N){
static inline void symv(const gene_matrix & A, const gene_vector & B, gene_vector & X, int /*N*/){
X.noalias() = (A.template selfadjointView<Lower>() * B);
// internal::product_selfadjoint_vector<real,0,LowerTriangularBit,false,false>(N,A.data(),N, B.data(), 1, X.data(), 1);
}
@ -155,54 +155,54 @@ public :
}
}
static EIGEN_DONT_INLINE void syr2(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){
static EIGEN_DONT_INLINE void syr2(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){
// internal::product_selfadjoint_rank2_update<real,0,LowerTriangularBit>(N,A.data(),N, X.data(), 1, Y.data(), 1, -1);
for(int j=0; j<N; ++j)
A.col(j).tail(N-j) += X[j] * Y.tail(N-j) + Y[j] * X.tail(N-j);
}
static EIGEN_DONT_INLINE void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){
static EIGEN_DONT_INLINE void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){
for(int j=0; j<N; ++j)
A.col(j) += X * Y[j];
}
static EIGEN_DONT_INLINE void rot(gene_vector & A, gene_vector & B, real c, real s, int N){
static EIGEN_DONT_INLINE void rot(gene_vector & A, gene_vector & B, real c, real s, int /*N*/){
internal::apply_rotation_in_the_plane(A, B, JacobiRotation<real>(c,s));
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int /*N*/){
X.noalias() = (A.transpose()*B);
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int /*N*/){
Y += coef * X;
}
static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){
static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int /*N*/){
Y = a*X + b*Y;
}
static EIGEN_DONT_INLINE void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
static EIGEN_DONT_INLINE void copy_matrix(const gene_matrix & source, gene_matrix & cible, int /*N*/){
cible = source;
}
static EIGEN_DONT_INLINE void copy_vector(const gene_vector & source, gene_vector & cible, int N){
static EIGEN_DONT_INLINE void copy_vector(const gene_vector & source, gene_vector & cible, int /*N*/){
cible = source;
}
static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector& X, int N){
static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector& X, int /*N*/){
X = L.template triangularView<Lower>().solve(B);
}
static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int /*N*/){
X = L.template triangularView<Upper>().solve(B);
}
static inline void trmm(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){
static inline void trmm(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int /*N*/){
X.noalias() = L.template triangularView<Lower>() * B;
}
static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){
static inline void cholesky(const gene_matrix & X, gene_matrix & C, int /*N*/){
C = X;
internal::llt_inplace<real,Lower>::blocked(C);
//C = X.llt().matrixL();
@ -211,11 +211,11 @@ public :
// Cholesky<gene_matrix>::computeInPlaceBlock(C);
}
static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int /*N*/){
C = X.fullPivLu().matrixLU();
}
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){
Matrix<DenseIndex,1,Dynamic> piv(N);
DenseIndex nb;
C = X;
@ -223,13 +223,13 @@ public :
// C = X.partialPivLu().matrixLU();
}
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){
typename Tridiagonalization<gene_matrix>::CoeffVectorType aux(N-1);
C = X;
internal::tridiagonalization_inplace(C, aux);
}
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int /*N*/){
C = HessenbergDecomposition<gene_matrix>(X).packedMatrix();
}