diff options
Diffstat (limited to 'bench')
23 files changed, 340 insertions, 86 deletions
diff --git a/bench/btl/actions/action_aat_product.hh b/bench/btl/actions/action_aat_product.hh index 3d11e7d6f..92930e219 100644 --- a/bench/btl/actions/action_aat_product.hh +++ b/bench/btl/actions/action_aat_product.hh @@ -104,7 +104,7 @@ public : } void check_result( void ){ - + if (_size>128) return; // calculation check Interface::matrix_to_stl(X,resu_stl); diff --git a/bench/btl/actions/action_ata_product.hh b/bench/btl/actions/action_ata_product.hh index 47531d607..04364fe67 100644 --- a/bench/btl/actions/action_ata_product.hh +++ b/bench/btl/actions/action_ata_product.hh @@ -104,7 +104,7 @@ public : } void check_result( void ){ - + if (_size>128) return; // calculation check Interface::matrix_to_stl(X,resu_stl); diff --git a/bench/btl/actions/action_atv_product.hh b/bench/btl/actions/action_atv_product.hh index 22d25e23e..a8234514b 100644 --- a/bench/btl/actions/action_atv_product.hh +++ b/bench/btl/actions/action_atv_product.hh @@ -93,6 +93,7 @@ public : void check_result( void ) { + if (_size>128) return; Interface::vector_to_stl(X,resu_stl); STL_interface<typename Interface::real_type>::atv_product(A_stl,B_stl,X_stl,_size); diff --git a/bench/btl/actions/action_axpby.hh b/bench/btl/actions/action_axpby.hh index 7815530ac..99ade7cd8 100644 --- a/bench/btl/actions/action_axpby.hh +++ b/bench/btl/actions/action_axpby.hh @@ -91,6 +91,7 @@ public : } void check_result( void ){ + if (_size>128) return; // calculation check Interface::vector_to_stl(Y,resu_stl); diff --git a/bench/btl/actions/action_axpy.hh b/bench/btl/actions/action_axpy.hh index f127426db..e4cb3a5bd 100644 --- a/bench/btl/actions/action_axpy.hh +++ b/bench/btl/actions/action_axpy.hh @@ -102,7 +102,7 @@ public : } void check_result( void ){ - + if (_size>128) return; // calculation check Interface::vector_to_stl(Y,resu_stl); diff --git a/bench/btl/actions/action_cholesky.hh b/bench/btl/actions/action_cholesky.hh index a1a30977c..2ec83f021 100644 --- a/bench/btl/actions/action_cholesky.hh +++ b/bench/btl/actions/action_cholesky.hh @@ -38,14 +38,13 @@ public : { MESSAGE("Action_cholesky Ctor"); - // STL vector initialization - typename Interface::stl_matrix tmp; - init_matrix<pseudo_random>(tmp,_size); - init_matrix<null_function>(X_stl,_size); - STL_interface<typename Interface::real_type>::ata_product(tmp,X_stl,_size); - + // STL mat/vec initialization + init_matrix_symm<pseudo_random>(X_stl,_size); init_matrix<null_function>(C_stl,_size); - init_matrix<null_function>(resu_stl,_size); + + // make sure X is invertible + for (int i=0; i<_size; ++i) + X_stl[i][i] = X_stl[i][i] * 1e2 + 1; // generic matrix and vector initialization Interface::matrix_from_stl(X_ref,X_stl); @@ -101,8 +100,6 @@ public : void check_result( void ){ // calculation check - Interface::matrix_to_stl(C,resu_stl); - // STL_interface<typename Interface::real_type>::cholesky(X_stl,C_stl,_size); // // typename Interface::real_type error= @@ -119,7 +116,6 @@ private : typename Interface::stl_matrix X_stl; typename Interface::stl_matrix C_stl; - typename Interface::stl_matrix resu_stl; typename Interface::gene_matrix X_ref; typename Interface::gene_matrix X; diff --git a/bench/btl/actions/action_ger.hh b/bench/btl/actions/action_ger.hh new file mode 100644 index 000000000..dc766efc5 --- /dev/null +++ b/bench/btl/actions/action_ger.hh @@ -0,0 +1,128 @@ + +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License +// as published by the Free Software Foundation; either version 2 +// of the License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// +#ifndef ACTION_GER +#define ACTION_GER +#include "utilities.h" +#include "STL_interface.hh" +#include <string> +#include "init/init_function.hh" +#include "init/init_vector.hh" +#include "init/init_matrix.hh" + +using namespace std; + +template<class Interface> +class Action_ger { + +public : + + // Ctor + BTL_DONT_INLINE Action_ger( int size ):_size(size) + { + MESSAGE("Action_ger Ctor"); + + // STL matrix and vector initialization + typename Interface::stl_matrix tmp; + init_matrix<pseudo_random>(A_stl,_size); + init_vector<pseudo_random>(B_stl,_size); + init_vector<pseudo_random>(X_stl,_size); + init_vector<null_function>(resu_stl,_size); + + // generic matrix and vector initialization + Interface::matrix_from_stl(A_ref,A_stl); + Interface::matrix_from_stl(A,A_stl); + Interface::vector_from_stl(B_ref,B_stl); + Interface::vector_from_stl(B,B_stl); + Interface::vector_from_stl(X_ref,X_stl); + Interface::vector_from_stl(X,X_stl); + } + + // invalidate copy ctor + Action_ger( const Action_ger & ) + { + INFOS("illegal call to Action_ger Copy Ctor"); + exit(1); + } + + // Dtor + BTL_DONT_INLINE ~Action_ger( void ){ + MESSAGE("Action_ger Dtor"); + Interface::free_matrix(A,_size); + Interface::free_vector(B); + Interface::free_vector(X); + Interface::free_matrix(A_ref,_size); + Interface::free_vector(B_ref); + Interface::free_vector(X_ref); + + } + + // action name + static inline std::string name( void ) + { + return "ger_" + Interface::name(); + } + + double nb_op_base( void ){ + return 2.0*_size*_size; + } + + BTL_DONT_INLINE void initialize( void ){ + Interface::copy_matrix(A_ref,A,_size); + Interface::copy_vector(B_ref,B,_size); + Interface::copy_vector(X_ref,X,_size); + } + + BTL_DONT_INLINE void calculate( void ) { + BTL_ASM_COMMENT("#begin ger"); + Interface::ger(A,B,X,_size); + BTL_ASM_COMMENT("end ger"); + } + + BTL_DONT_INLINE void check_result( void ){ + // calculation check + Interface::vector_to_stl(X,resu_stl); + + STL_interface<typename Interface::real_type>::ger(A_stl,B_stl,X_stl,_size); + + typename Interface::real_type error= + STL_interface<typename Interface::real_type>::norm_diff(X_stl,resu_stl); + + if (error>1.e-3){ + INFOS("WRONG CALCULATION...residual=" << error); +// exit(0); + } + + } + +private : + + typename Interface::stl_matrix A_stl; + typename Interface::stl_vector B_stl; + typename Interface::stl_vector X_stl; + typename Interface::stl_vector resu_stl; + + typename Interface::gene_matrix A_ref; + typename Interface::gene_vector B_ref; + typename Interface::gene_vector X_ref; + + typename Interface::gene_matrix A; + typename Interface::gene_vector B; + typename Interface::gene_vector X; + + int _size; +}; + + +#endif diff --git a/bench/btl/actions/action_partial_lu.hh b/bench/btl/actions/action_partial_lu.hh new file mode 100644 index 000000000..837661027 --- /dev/null +++ b/bench/btl/actions/action_partial_lu.hh @@ -0,0 +1,125 @@ +//===================================================== +// File : action_lu_decomp.hh +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +//===================================================== +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License +// as published by the Free Software Foundation; either version 2 +// of the License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// +#ifndef ACTION_PARTIAL_LU +#define ACTION_PARTIAL_LU +#include "utilities.h" +#include "STL_interface.hh" +#include <string> +#include "init/init_function.hh" +#include "init/init_vector.hh" +#include "init/init_matrix.hh" + +using namespace std; + +template<class Interface> +class Action_partial_lu { + +public : + + // Ctor + + Action_partial_lu( int size ):_size(size) + { + MESSAGE("Action_partial_lu Ctor"); + + // STL vector initialization + init_matrix<pseudo_random>(X_stl,_size); + init_matrix<null_function>(C_stl,_size); + + // make sure X is invertible + for (int i=0; i<_size; ++i) + X_stl[i][i] = X_stl[i][i] * 1e2 + 1; + + // generic matrix and vector initialization + Interface::matrix_from_stl(X_ref,X_stl); + Interface::matrix_from_stl(X,X_stl); + Interface::matrix_from_stl(C,C_stl); + + _cost = 2.0*size*size*size/3.0 + size*size; + } + + // invalidate copy ctor + + Action_partial_lu( const Action_partial_lu & ) + { + INFOS("illegal call to Action_partial_lu Copy Ctor"); + exit(1); + } + + // Dtor + + ~Action_partial_lu( void ){ + + MESSAGE("Action_partial_lu Dtor"); + + // deallocation + Interface::free_matrix(X_ref,_size); + Interface::free_matrix(X,_size); + Interface::free_matrix(C,_size); + } + + // action name + + static inline std::string name( void ) + { + return "partial_lu_decomp_"+Interface::name(); + } + + double nb_op_base( void ){ + return _cost; + } + + inline void initialize( void ){ + Interface::copy_matrix(X_ref,X,_size); + } + + inline void calculate( void ) { + Interface::partial_lu_decomp(X,C,_size); + } + + void check_result( void ){ + // calculation check +// Interface::matrix_to_stl(C,resu_stl); + +// STL_interface<typename Interface::real_type>::lu_decomp(X_stl,C_stl,_size); +// +// typename Interface::real_type error= +// STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl); +// +// if (error>1.e-6){ +// INFOS("WRONG CALCULATION...residual=" << error); +// exit(0); +// } + + } + +private : + + typename Interface::stl_matrix X_stl; + typename Interface::stl_matrix C_stl; + + typename Interface::gene_matrix X_ref; + typename Interface::gene_matrix X; + typename Interface::gene_matrix C; + + int _size; + double _cost; +}; + +#endif diff --git a/bench/btl/actions/action_symv.hh b/bench/btl/actions/action_symv.hh index f520623a8..a32b9dfa0 100644 --- a/bench/btl/actions/action_symv.hh +++ b/bench/btl/actions/action_symv.hh @@ -40,18 +40,12 @@ public : MESSAGE("Action_symv Ctor"); // STL matrix and vector initialization - - typename Interface::stl_matrix tmp; - init_matrix<pseudo_random>(A_stl,_size); - init_matrix<pseudo_random>(tmp,_size); + init_matrix_symm<pseudo_random>(A_stl,_size); init_vector<pseudo_random>(B_stl,_size); init_vector<null_function>(X_stl,_size); init_vector<null_function>(resu_stl,_size); - STL_interface<typename Interface::real_type>::ata_product(tmp,A_stl,_size); - // generic matrix and vector initialization - Interface::matrix_from_stl(A_ref,A_stl); Interface::matrix_from_stl(A,A_stl); Interface::vector_from_stl(B_ref,B_stl); @@ -70,21 +64,13 @@ public : } // Dtor - BTL_DONT_INLINE ~Action_symv( void ){ - - MESSAGE("Action_symv Dtor"); - - // deallocation - Interface::free_matrix(A,_size); Interface::free_vector(B); Interface::free_vector(X); - Interface::free_matrix(A_ref,_size); Interface::free_vector(B_ref); Interface::free_vector(X_ref); - } // action name @@ -113,9 +99,8 @@ public : } BTL_DONT_INLINE void check_result( void ){ - + if (_size>128) return; // calculation check - Interface::vector_to_stl(X,resu_stl); STL_interface<typename Interface::real_type>::symv(A_stl,B_stl,X_stl,_size); diff --git a/bench/btl/actions/action_syr2.hh b/bench/btl/actions/action_syr2.hh index d5ef07411..7c6712b13 100644 --- a/bench/btl/actions/action_syr2.hh +++ b/bench/btl/actions/action_syr2.hh @@ -37,32 +37,23 @@ public : BTL_DONT_INLINE Action_syr2( int size ):_size(size) { - MESSAGE("Action_syr2 Ctor"); - // STL matrix and vector initialization - typename Interface::stl_matrix tmp; init_matrix<pseudo_random>(A_stl,_size); - init_matrix<pseudo_random>(tmp,_size); init_vector<pseudo_random>(B_stl,_size); init_vector<pseudo_random>(X_stl,_size); init_vector<null_function>(resu_stl,_size); - STL_interface<typename Interface::real_type>::ata_product(tmp,A_stl,_size); - // generic matrix and vector initialization - Interface::matrix_from_stl(A_ref,A_stl); Interface::matrix_from_stl(A,A_stl); Interface::vector_from_stl(B_ref,B_stl); Interface::vector_from_stl(B,B_stl); Interface::vector_from_stl(X_ref,X_stl); Interface::vector_from_stl(X,X_stl); - } // invalidate copy ctor - Action_syr2( const Action_syr2 & ) { INFOS("illegal call to Action_syr2 Copy Ctor"); @@ -70,21 +61,13 @@ public : } // Dtor - BTL_DONT_INLINE ~Action_syr2( void ){ - - MESSAGE("Action_syr2 Dtor"); - - // deallocation - Interface::free_matrix(A,_size); Interface::free_vector(B); Interface::free_vector(X); - Interface::free_matrix(A_ref,_size); Interface::free_vector(B_ref); Interface::free_vector(X_ref); - } // action name @@ -99,11 +82,9 @@ public : } BTL_DONT_INLINE void initialize( void ){ - Interface::copy_matrix(A_ref,A,_size); Interface::copy_vector(B_ref,B,_size); Interface::copy_vector(X_ref,X,_size); - } BTL_DONT_INLINE void calculate( void ) { @@ -113,9 +94,7 @@ public : } BTL_DONT_INLINE void check_result( void ){ - // calculation check -return; Interface::vector_to_stl(X,resu_stl); STL_interface<typename Interface::real_type>::syr2(A_stl,B_stl,X_stl,_size); diff --git a/bench/btl/actions/action_trisolve.hh b/bench/btl/actions/action_trisolve.hh index 0d277eb77..94e41e072 100644 --- a/bench/btl/actions/action_trisolve.hh +++ b/bench/btl/actions/action_trisolve.hh @@ -102,7 +102,8 @@ public : Interface::trisolve_lower(L,B,X,_size); } - void check_result( void ){ + void check_result(){ + if (_size>128) return; // calculation check Interface::vector_to_stl(X,resu_stl); diff --git a/bench/btl/actions/basic_actions.hh b/bench/btl/actions/basic_actions.hh index a23e58096..06a382be6 100644 --- a/bench/btl/actions/basic_actions.hh +++ b/bench/btl/actions/basic_actions.hh @@ -14,6 +14,7 @@ #include "action_symv.hh" // #include "action_symm.hh" #include "action_syr2.hh" +#include "action_ger.hh" // #include "action_lu_solve.hh" diff --git a/bench/btl/data/action_settings.txt b/bench/btl/data/action_settings.txt index ab1880e84..e2d931ddc 100644 --- a/bench/btl/data/action_settings.txt +++ b/bench/btl/data/action_settings.txt @@ -1,16 +1,17 @@ -aat ; "{/*1.5 A x A^T}" ; "matrix size" ; 4:2048 -ata ; "{/*1.5 A^T x A}" ; "matrix size" ; 4:2048 -atv ; "{/*1.5 matrix^T x vector}" ; "matrix size" ; 4:2048 +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 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:2048 -matrix_vector ; "{/*1.5 matrix vector product}" ; "matrix size" ; 4:2048 -trisolve ; "{/*1.5 triangular solver (X = inv(L) X)}" ; "size" ; 4:2048 -matrix_trisolve ; "{/*1.5 matrix triangular solver (M = inv(L) M)}" ; "size" ; 4:2048 -cholesky ; "{/*1.5 Cholesky decomposition}" ; "matrix size" ; 4:2048 -lu_decomp ; "{/*1.5 Complete LU decomposition}" ; "matrix size" ; 4:2048 -partial_lu_decomp ; "{/*1.5 Partial LU decomposition}" ; "matrix size" ; 4:2048 -tridiagonalization ; "{/*1.5 Tridiagonalization}" ; "matrix size" ; 4:2048 -hessenberg ; "{/*1.5 Hessenberg decomposition}" ; "matrix size" ; 4:2048 -symv ; "{/*1.5 symmetric matrix vector product}" ; "matrix size" ; 4:2048 -syr2 ; "{/*1.5 symmetric rank-2 update (A += u^T v + u v^T)}" ; "matrix size" ; 4:2048
\ No newline at end of file +matrix_matrix ; "{/*1.5 matrix matrix product}" ; "matrix size" ; 4:3000 +matrix_vector ; "{/*1.5 matrix vector product}" ; "matrix size" ; 4:3000 +trisolve ; "{/*1.5 triangular solver (X = inv(L) X)}" ; "size" ; 4:3000 +matrix_trisolve ; "{/*1.5 matrix triangular solver (M = inv(L) M)}" ; "size" ; 4:3000 +cholesky ; "{/*1.5 Cholesky decomposition}" ; "matrix size" ; 4:3000 +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
\ No newline at end of file diff --git a/bench/btl/data/go_mean b/bench/btl/data/go_mean index 38c50d226..5ff3a3b98 100755 --- a/bench/btl/data/go_mean +++ b/bench/btl/data/go_mean @@ -47,6 +47,7 @@ source mk_mean_script.sh tridiagonalization $1 11 100 300 1000 $mode $prefix source mk_mean_script.sh hessenberg $1 11 100 300 1000 $mode $prefix source mk_mean_script.sh symv $1 11 50 300 1000 $mode $prefix source mk_mean_script.sh syr2 $1 11 50 300 1000 $mode $prefix +source mk_mean_script.sh ger $1 11 50 300 1000 $mode $prefix fi diff --git a/bench/btl/generic_bench/bench_parameter.hh b/bench/btl/generic_bench/bench_parameter.hh index 08fea80e4..e9603e4fc 100644 --- a/bench/btl/generic_bench/bench_parameter.hh +++ b/bench/btl/generic_bench/bench_parameter.hh @@ -23,7 +23,7 @@ // minimal time for each measurement #define REAL_TYPE float // minimal time for each measurement -#define MIN_TIME 1.0 +#define MIN_TIME 0.5 // nb of point on bench curves #define NB_POINT 100 // min vector size for axpy bench @@ -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 2048 +#define MAX_MV 3000 // min matrix size for matrix matrix product bench #define MIN_MM 5 // max matrix size for matrix matrix product bench @@ -41,13 +41,13 @@ // min matrix size for LU bench #define MIN_LU 5 // max matrix size for LU bench -#define MAX_LU 2048 +#define MAX_LU 3000 // max size for tiny vector and matrix #define TINY_MV_MAX_SIZE 16 // default nb_sample for x86 timer #define DEFAULT_NB_SAMPLE 1000 // how many times we run a single bench (keep the best perf) -#define NB_TRIES 3 +#define NB_TRIES 5 #endif diff --git a/bench/btl/generic_bench/init/init_function.hh b/bench/btl/generic_bench/init/init_function.hh index c7b348957..7b3bdbafc 100644 --- a/bench/btl/generic_bench/init/init_function.hh +++ b/bench/btl/generic_bench/init/init_function.hh @@ -32,7 +32,6 @@ double simple_function(int index_i, int index_j) double pseudo_random(int index) { - // INFOS("random="<<(std::rand()/double(RAND_MAX))); return std::rand()/double(RAND_MAX); } diff --git a/bench/btl/generic_bench/init/init_matrix.hh b/bench/btl/generic_bench/init/init_matrix.hh index 6b57504c0..67cbd2073 100644 --- a/bench/btl/generic_bench/init/init_matrix.hh +++ b/bench/btl/generic_bench/init/init_matrix.hh @@ -41,13 +41,24 @@ BTL_DONT_INLINE void init_row(Vector & X, int size, int row){ // [] operator for setting rows template<double init_function(int,int),class Vector> BTL_DONT_INLINE void init_matrix(Vector & A, int size){ - A.resize(size); for (int row=0; row<A.size() ; row++){ init_row<init_function>(A[row],size,row); } +} - +template<double init_function(int,int),class Matrix> +BTL_DONT_INLINE void init_matrix_symm(Matrix& A, int size){ + A.resize(size); + for (int row=0; row<A.size() ; row++) + A[row].resize(size); + for (int row=0; row<A.size() ; row++){ + A[row][row] = init_function(row,row); + for (int col=0; col<row ; col++){ + double x = init_function(row,col); + A[row][col] = A[col][row] = x; + } + } } #endif diff --git a/bench/btl/generic_bench/timers/portable_perf_analyzer.hh b/bench/btl/generic_bench/timers/portable_perf_analyzer.hh index 67d3378fc..4298e61df 100644 --- a/bench/btl/generic_bench/timers/portable_perf_analyzer.hh +++ b/bench/btl/generic_bench/timers/portable_perf_analyzer.hh @@ -47,7 +47,6 @@ public: time_action = time_calculate(action); while (time_action < MIN_TIME) { - //Action action(size); _nb_calc *= 2; action.initialize(); time_action = time_calculate(action); @@ -65,7 +64,7 @@ public: time_action = time_action / (double(_nb_calc)); // check - if (BtlConfig::Instance.checkResults) + if (BtlConfig::Instance.checkResults && size<128) { action.initialize(); action.calculate(); diff --git a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh index c8b69883e..725f35944 100644 --- a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh +++ b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh @@ -171,6 +171,14 @@ public : #endif } + static inline void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){ + #ifdef PUREBLAS + sger_(&N,&N,&fone,X,&intone,Y,&intone,A,&N); + #else + cblas_sger(CblasColMajor,N,N,1.0,X,1,Y,1,A,N); + #endif + } + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ #ifdef PUREBLAS sgemv_(&trans,&N,&N,&fone,A,&N,B,&intone,&fzero,X,&intone); diff --git a/bench/btl/libs/C_BLAS/main.cpp b/bench/btl/libs/C_BLAS/main.cpp index bb3022d3f..6fd75812c 100644 --- a/bench/btl/libs/C_BLAS/main.cpp +++ b/bench/btl/libs/C_BLAS/main.cpp @@ -44,6 +44,8 @@ int main() bench<Action_symv<C_BLAS_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); bench<Action_syr2<C_BLAS_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_ger<C_BLAS_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_matrix_matrix_product<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); bench<Action_ata_product<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); bench<Action_aat_product<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); diff --git a/bench/btl/libs/STL/STL_interface.hh b/bench/btl/libs/STL/STL_interface.hh index 3958d4af5..0b73382f3 100644 --- a/bench/btl/libs/STL/STL_interface.hh +++ b/bench/btl/libs/STL/STL_interface.hh @@ -156,6 +156,15 @@ public : } } + static inline void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N) + { + for (int j=0; j<N; ++j) + { + for (int i=j; i<N; ++i) + A[j][i] += X[i]*Y[j]; + } + } + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) { real somme; diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh index 7a014723f..6a1bc5d61 100644 --- a/bench/btl/libs/eigen2/eigen2_interface.hh +++ b/bench/btl/libs/eigen2/eigen2_interface.hh @@ -104,12 +104,12 @@ public : } static inline void matrix_vector_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N){ - X = (A*B)/*.lazy()*/; + X = (A*B).lazy(); } static inline void symv(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N){ - //X = (A.template marked<SelfAdjoint|LowerTriangular>() * B)/*.lazy()*/; - ei_product_selfadjoint_vector<real,0,LowerTriangularBit,false,false>(N,A.data(),N, B.data(), 1, X.data(), 1); + X = (A.template selfadjointView<LowerTriangular>() * B)/*.lazy()*/; +// ei_product_selfadjoint_vector<real,0,LowerTriangularBit,false,false>(N,A.data(),N, B.data(), 1, X.data(), 1); } template<typename Dest, typename Src> static void triassign(Dest& dst, const Src& src) @@ -163,8 +163,13 @@ public : A.col(j).end(N-j) += X[j] * Y.end(N-j) + Y[j] * X.end(N-j); } + 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 inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ - X = (A.transpose()*B)/*.lazy()*/; + X = (A.transpose()*B).lazy(); } static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){ @@ -172,9 +177,7 @@ public : } static inline void axpby(real a, const gene_vector & X, real b, gene_vector & Y, int N){ - asm("#begin axpby"); Y = a*X + b*Y; - asm("#end axpby"); } static EIGEN_DONT_INLINE void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ @@ -191,20 +194,23 @@ public : static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix& X, int N){ X = L.template triangularView<LowerTriangular>().solve(B); -// -// ei_triangular_solve_matrix<real,ColMajor,ColMajor,LowerTriangular> -// ::run(L.cols(), X.cols(), L.data(), L.stride(), X.data(), X.stride()); } static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){ - C = X.llt().matrixL(); + C = X; + ei_llt_inplace<LowerTriangular>::blocked(C); + //C = X.llt().matrixL(); // C = X; // Cholesky<gene_matrix>::computeInPlace(C); // Cholesky<gene_matrix>::computeInPlaceBlock(C); } static inline void lu_decomp(const gene_matrix & X, gene_matrix & C, int N){ - C = X.lu().matrixLU(); + RowVectorXi piv(N); + int nb; + C = X; + ei_partial_lu_inplace(C,piv,nb); + //C = X.lu().matrixLU(); } static inline void partial_lu_decomp(const gene_matrix & X, gene_matrix & C, int N){ diff --git a/bench/btl/libs/eigen2/main_vecmat.cpp b/bench/btl/libs/eigen2/main_vecmat.cpp index fb00d6f79..b4434d4af 100644 --- a/bench/btl/libs/eigen2/main_vecmat.cpp +++ b/bench/btl/libs/eigen2/main_vecmat.cpp @@ -28,6 +28,7 @@ int main() bench<Action_atv_product<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); bench<Action_symv<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); bench<Action_syr2<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_ger<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); return 0; } |