diff options
author | Gael Guennebaud <g.gael@free.fr> | 2008-07-09 14:04:48 +0000 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2008-07-09 14:04:48 +0000 |
commit | 28539e7597c643dbc2b8d4f49dd16bd86fb7251f (patch) | |
tree | 3bfb96ae898f0afc943a388edcc83c27d28d5b3a /bench/btl/libs/gmm | |
parent | 5f55ab524cf0aad34dd6884bcae09eaa6c43c247 (diff) |
imported a reworked version of BTL (Benchmark for Templated Libraries).
the modifications to initial code follow:
* changed build system from plain makefiles to cmake
* added eigen2 (4 versions: vec/novec and fixed/dynamic), GMM++, MTL4 interfaces
* added "transposed matrix * vector" product action
* updated blitz interface to use condensed products instead of hand coded loops
* removed some deprecated interfaces
* changed default storage order to column major for all libraries
* new generic bench timer strategy which is supposed to be more accurate
* various code clean-up
Diffstat (limited to 'bench/btl/libs/gmm')
-rw-r--r-- | bench/btl/libs/gmm/CMakeLists.txt | 3 | ||||
-rw-r--r-- | bench/btl/libs/gmm/gmm_LU_solve_interface.hh | 192 | ||||
-rw-r--r-- | bench/btl/libs/gmm/gmm_interface.hh | 119 | ||||
-rw-r--r-- | bench/btl/libs/gmm/main.cpp | 44 |
4 files changed, 358 insertions, 0 deletions
diff --git a/bench/btl/libs/gmm/CMakeLists.txt b/bench/btl/libs/gmm/CMakeLists.txt new file mode 100644 index 000000000..98759e8b6 --- /dev/null +++ b/bench/btl/libs/gmm/CMakeLists.txt @@ -0,0 +1,3 @@ + +include_directories(${GMM_INCLUDES}) +btl_add_bench(btl_gmm main.cpp) diff --git a/bench/btl/libs/gmm/gmm_LU_solve_interface.hh b/bench/btl/libs/gmm/gmm_LU_solve_interface.hh new file mode 100644 index 000000000..dcb9f567f --- /dev/null +++ b/bench/btl/libs/gmm/gmm_LU_solve_interface.hh @@ -0,0 +1,192 @@ +//===================================================== +// File : blitz_LU_solve_interface.hh +// Author : L. Plagne <laurent.plagne@edf.fr)> +// Copyright (C) EDF R&D, lun sep 30 14:23:31 CEST 2002 +//===================================================== +// +// 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 BLITZ_LU_SOLVE_INTERFACE_HH +#define BLITZ_LU_SOLVE_INTERFACE_HH + +#include "blitz/array.h" +#include <vector> + +BZ_USING_NAMESPACE(blitz) + +template<class real> +class blitz_LU_solve_interface : public blitz_interface<real> +{ + +public : + + typedef typename blitz_interface<real>::gene_matrix gene_matrix; + typedef typename blitz_interface<real>::gene_vector gene_vector; + + typedef blitz::Array<int,1> Pivot_Vector; + + inline static void new_Pivot_Vector(Pivot_Vector & pivot,int N) + { + + pivot.resize(N); + + } + + inline static void free_Pivot_Vector(Pivot_Vector & pivot) + { + + return; + + } + + + static inline real matrix_vector_product_sliced(const gene_matrix & A, gene_vector B, int row, int col_start, int col_end) + { + + real somme=0.; + + for (int j=col_start ; j<col_end+1 ; j++){ + + somme+=A(row,j)*B(j); + + } + + return somme; + + } + + + + + static inline real matrix_matrix_product_sliced(gene_matrix & A, int row, int col_start, int col_end, gene_matrix & B, int row_shift, int col ) + { + + real somme=0.; + + for (int j=col_start ; j<col_end+1 ; j++){ + + somme+=A(row,j)*B(j+row_shift,col); + + } + + return somme; + + } + + inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N) + { + + ASSERT( LU.rows()==LU.cols() ) ; + int index_max = 0 ; + real big = 0. ; + real theSum = 0. ; + real dum = 0. ; + // Get the implicit scaling information : + gene_vector ImplicitScaling( N ) ; + for( int i=0; i<N; i++ ) { + big = 0. ; + for( int j=0; j<N; j++ ) { + if( abs( LU( i, j ) )>=big ) big = abs( LU( i, j ) ) ; + } + if( big==0. ) { + INFOS( "blitz_LU_factor::Singular matrix" ) ; + exit( 0 ) ; + } + ImplicitScaling( i ) = 1./big ; + } + // Loop over columns of Crout's method : + for( int j=0; j<N; j++ ) { + for( int i=0; i<j; i++ ) { + theSum = LU( i, j ) ; + theSum -= matrix_matrix_product_sliced(LU, i, 0, i-1, LU, 0, j) ; + // theSum -= sum( LU( i, Range( fromStart, i-1 ) )*LU( Range( fromStart, i-1 ), j ) ) ; + LU( i, j ) = theSum ; + } + + // Search for the largest pivot element : + big = 0. ; + for( int i=j; i<N; i++ ) { + theSum = LU( i, j ) ; + theSum -= matrix_matrix_product_sliced(LU, i, 0, j-1, LU, 0, j) ; + // theSum -= sum( LU( i, Range( fromStart, j-1 ) )*LU( Range( fromStart, j-1 ), j ) ) ; + LU( i, j ) = theSum ; + if( (ImplicitScaling( i )*abs( theSum ))>=big ) { + dum = ImplicitScaling( i )*abs( theSum ) ; + big = dum ; + index_max = i ; + } + } + // Interchanging rows and the scale factor : + if( j!=index_max ) { + for( int k=0; k<N; k++ ) { + dum = LU( index_max, k ) ; + LU( index_max, k ) = LU( j, k ) ; + LU( j, k ) = dum ; + } + ImplicitScaling( index_max ) = ImplicitScaling( j ) ; + } + pivot( j ) = index_max ; + if ( LU( j, j )==0. ) LU( j, j ) = 1.e-20 ; + // Divide by the pivot element : + if( j<N ) { + dum = 1./LU( j, j ) ; + for( int i=j+1; i<N; i++ ) LU( i, j ) *= dum ; + } + } + + } + + inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, gene_vector &B, gene_vector X, int N) + { + + // Pour conserver le meme header, on travaille sur X, copie du second-membre B + X = B.copy() ; + ASSERT( LU.rows()==LU.cols() ) ; + firstIndex indI ; + // Forward substitution : + int ii = 0 ; + real theSum = 0. ; + for( int i=0; i<N; i++ ) { + int ip = pivot( i ) ; + theSum = X( ip ) ; + // theSum = B( ip ) ; + X( ip ) = X( i ) ; + // B( ip ) = B( i ) ; + if( ii ) { + theSum -= matrix_vector_product_sliced(LU, X, i, ii-1, i-1) ; + // theSum -= sum( LU( i, Range( ii-1, i-1 ) )*X( Range( ii-1, i-1 ) ) ) ; + // theSum -= sum( LU( i, Range( ii-1, i-1 ) )*B( Range( ii-1, i-1 ) ) ) ; + } else if( theSum ) { + ii = i+1 ; + } + X( i ) = theSum ; + // B( i ) = theSum ; + } + // Backsubstitution : + for( int i=N-1; i>=0; i-- ) { + theSum = X( i ) ; + // theSum = B( i ) ; + theSum -= matrix_vector_product_sliced(LU, X, i, i+1, N) ; + // theSum -= sum( LU( i, Range( i+1, toEnd ) )*X( Range( i+1, toEnd ) ) ) ; + // theSum -= sum( LU( i, Range( i+1, toEnd ) )*B( Range( i+1, toEnd ) ) ) ; + // Store a component of the solution vector : + X( i ) = theSum/LU( i, i ) ; + // B( i ) = theSum/LU( i, i ) ; + } + + } + +}; + +#endif diff --git a/bench/btl/libs/gmm/gmm_interface.hh b/bench/btl/libs/gmm/gmm_interface.hh new file mode 100644 index 000000000..6a81fe969 --- /dev/null +++ b/bench/btl/libs/gmm/gmm_interface.hh @@ -0,0 +1,119 @@ +//===================================================== +// 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 GMM_INTERFACE_HH +#define GMM_INTERFACE_HH + +#include <gmm/gmm.h> +#include <vector> + +using namespace gmm; + +template<class real> +class gmm_interface { + +public : + + typedef real real_type ; + + typedef std::vector<real> stl_vector; + typedef std::vector<stl_vector > stl_matrix; + + typedef gmm::dense_matrix<real> gene_matrix; + typedef stl_vector gene_vector; + + static inline std::string name( void ) + { + return "gmm"; + } + + static void free_matrix(gene_matrix & A, int N){ + return ; + } + + static void free_vector(gene_vector & B){ + return ; + } + + static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ + A.resize(A_stl[0].size(),A_stl.size()); + + for (int j=0; j<A_stl.size() ; j++){ + for (int i=0; i<A_stl[j].size() ; i++){ + A(i,j) = A_stl[j][i]; + } + } + } + + static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){ + B = B_stl; + } + + static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ + B_stl = B; + } + + static inline void matrix_to_stl(gene_matrix & A, stl_matrix & A_stl){ + int N=A_stl.size(); + + for (int j=0;j<N;j++){ + A_stl[j].resize(N); + for (int i=0;i<N;i++){ + A_stl[j][i] = A(i,j); + } + } + } + + static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ + gmm::mult(A,B, X); + } + + static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ + gmm::mult(gmm::transposed(A),gmm::transposed(B), X); + } + + static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){ + gmm::mult(gmm::transposed(A),A, X); + } + + static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){ + gmm::mult(A,gmm::transposed(A), X); + } + + static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ + gmm::mult(A,B,X); + } + + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ + gmm::mult(gmm::transposed(A),B,X); + } + + static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){ + gmm::add(gmm::scaled(X,coef), Y); + } + + static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ + gmm::copy(source,cible); + } + + static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ + gmm::copy(source,cible); + } + +}; + +#endif diff --git a/bench/btl/libs/gmm/main.cpp b/bench/btl/libs/gmm/main.cpp new file mode 100644 index 000000000..27afeedbb --- /dev/null +++ b/bench/btl/libs/gmm/main.cpp @@ -0,0 +1,44 @@ +//===================================================== +// 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. +// +#include "utilities.h" +#include "gmm_interface.hh" +#include "bench.hh" +#include "action_matrix_vector_product.hh" +#include "action_matrix_matrix_product.hh" +#include "action_axpy.hh" +#include "action_lu_solve.hh" +#include "action_ata_product.hh" +#include "action_aat_product.hh" +#include "action_atv_product.hh" + +int main() +{ + + bench<Action_axpy<gmm_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT); + bench<Action_matrix_vector_product<gmm_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_atv_product<gmm_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_matrix_matrix_product<gmm_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + bench<Action_ata_product<gmm_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + bench<Action_aat_product<gmm_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + + //bench<Action_lu_solve<blitz_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT); + + return 0; +} + + |