diff options
Diffstat (limited to 'eigen/bench/btl/libs/blitz')
-rw-r--r-- | eigen/bench/btl/libs/blitz/CMakeLists.txt | 17 | ||||
-rw-r--r-- | eigen/bench/btl/libs/blitz/blitz_LU_solve_interface.hh | 192 | ||||
-rw-r--r-- | eigen/bench/btl/libs/blitz/blitz_interface.hh | 147 | ||||
-rw-r--r-- | eigen/bench/btl/libs/blitz/btl_blitz.cpp | 51 | ||||
-rw-r--r-- | eigen/bench/btl/libs/blitz/btl_tiny_blitz.cpp | 38 | ||||
-rw-r--r-- | eigen/bench/btl/libs/blitz/tiny_blitz_interface.hh | 106 |
6 files changed, 551 insertions, 0 deletions
diff --git a/eigen/bench/btl/libs/blitz/CMakeLists.txt b/eigen/bench/btl/libs/blitz/CMakeLists.txt new file mode 100644 index 0000000..880ab73 --- /dev/null +++ b/eigen/bench/btl/libs/blitz/CMakeLists.txt @@ -0,0 +1,17 @@ + +find_package(Blitz) + +if (BLITZ_FOUND) + include_directories(${BLITZ_INCLUDES}) + + btl_add_bench(btl_blitz btl_blitz.cpp) + if (BUILD_btl_blitz) + target_link_libraries(btl_blitz ${BLITZ_LIBRARIES}) + endif (BUILD_btl_blitz) + + btl_add_bench(btl_tiny_blitz btl_tiny_blitz.cpp OFF) + if (BUILD_btl_tiny_blitz) + target_link_libraries(btl_tiny_blitz ${BLITZ_LIBRARIES}) + endif (BUILD_btl_tiny_blitz) + +endif (BLITZ_FOUND) diff --git a/eigen/bench/btl/libs/blitz/blitz_LU_solve_interface.hh b/eigen/bench/btl/libs/blitz/blitz_LU_solve_interface.hh new file mode 100644 index 0000000..dcb9f56 --- /dev/null +++ b/eigen/bench/btl/libs/blitz/blitz_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/eigen/bench/btl/libs/blitz/blitz_interface.hh b/eigen/bench/btl/libs/blitz/blitz_interface.hh new file mode 100644 index 0000000..a67c47c --- /dev/null +++ b/eigen/bench/btl/libs/blitz/blitz_interface.hh @@ -0,0 +1,147 @@ +//===================================================== +// File : blitz_interface.hh +// Author : L. Plagne <laurent.plagne@edf.fr)> +// Copyright (C) EDF R&D, lun sep 30 14:23:30 CEST 2002 +// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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 BLITZ_INTERFACE_HH +#define BLITZ_INTERFACE_HH + +#include <blitz/blitz.h> +#include <blitz/array.h> +#include <blitz/vector-et.h> +#include <blitz/vecwhere.h> +#include <blitz/matrix.h> +#include <vector> + +BZ_USING_NAMESPACE(blitz) + +template<class real> +class blitz_interface{ + +public : + + typedef real real_type ; + + typedef std::vector<real> stl_vector; + typedef std::vector<stl_vector > stl_matrix; + + typedef blitz::Array<real, 2> gene_matrix; + typedef blitz::Array<real, 1> gene_vector; +// typedef blitz::Matrix<real, blitz::ColumnMajor> gene_matrix; +// typedef blitz::Vector<real> gene_vector; + + static inline std::string name() { return "blitz"; } + + static void free_matrix(gene_matrix & A, int N){} + + static void free_vector(gene_vector & B){} + + 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.resize(B_stl.size()); + for (int i=0; i<B_stl.size() ; i++){ + B(i)=B_stl[i]; + } + } + + static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ + for (int i=0; i<B_stl.size() ; i++){ + B_stl[i]=B(i); + } + } + + 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) + { + firstIndex i; + secondIndex j; + thirdIndex k; + X = sum(A(i,k) * B(k,j), k); + } + + static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N) + { + firstIndex i; + secondIndex j; + thirdIndex k; + X = sum(A(k,i) * A(k,j), k); + } + + static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N) + { + firstIndex i; + secondIndex j; + thirdIndex k; + X = sum(A(i,k) * A(j,k), k); + } + + static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) + { + firstIndex i; + secondIndex j; + X = sum(A(i,j)*B(j),j); + } + + static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) + { + firstIndex i; + secondIndex j; + X = sum(A(j,i) * B(j),j); + } + + static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N) + { + firstIndex i; + Y = Y(i) + coef * X(i); + //Y += coef * X; + } + + static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ + cible = source; + //cible.template operator=<gene_matrix>(source); +// for (int i=0;i<N;i++){ +// for (int j=0;j<N;j++){ +// cible(i,j)=source(i,j); +// } +// } + } + + static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ + //cible.template operator=<gene_vector>(source); + cible = source; + } + +}; + +#endif diff --git a/eigen/bench/btl/libs/blitz/btl_blitz.cpp b/eigen/bench/btl/libs/blitz/btl_blitz.cpp new file mode 100644 index 0000000..16d2b59 --- /dev/null +++ b/eigen/bench/btl/libs/blitz/btl_blitz.cpp @@ -0,0 +1,51 @@ +//===================================================== +// File : main.cpp +// Author : L. Plagne <laurent.plagne@edf.fr)> +// Copyright (C) EDF R&D, lun sep 30 14:23:30 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. +// +#include "utilities.h" +#include "blitz_interface.hh" +#include "blitz_LU_solve_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" + +BTL_MAIN; + +int main() +{ + + bench<Action_matrix_vector_product<blitz_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + bench<Action_atv_product<blitz_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT); + + bench<Action_matrix_matrix_product<blitz_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + bench<Action_ata_product<blitz_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + bench<Action_aat_product<blitz_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT); + + bench<Action_axpy<blitz_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT); + + //bench<Action_lu_solve<blitz_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT); + + return 0; +} + + diff --git a/eigen/bench/btl/libs/blitz/btl_tiny_blitz.cpp b/eigen/bench/btl/libs/blitz/btl_tiny_blitz.cpp new file mode 100644 index 0000000..9fddde7 --- /dev/null +++ b/eigen/bench/btl/libs/blitz/btl_tiny_blitz.cpp @@ -0,0 +1,38 @@ +//===================================================== +// File : main.cpp +// Author : L. Plagne <laurent.plagne@edf.fr)> +// Copyright (C) EDF R&D, lun sep 30 14:23:30 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. +// +#include "utilities.h" +#include "tiny_blitz_interface.hh" +#include "static/bench_static.hh" +#include "action_matrix_vector_product.hh" +#include "action_matrix_matrix_product.hh" +#include "action_axpy.hh" + +BTL_MAIN; + +int main() +{ + bench_static<Action_axpy,tiny_blitz_interface>(); + bench_static<Action_matrix_matrix_product,tiny_blitz_interface>(); + bench_static<Action_matrix_vector_product,tiny_blitz_interface>(); + + return 0; +} + + diff --git a/eigen/bench/btl/libs/blitz/tiny_blitz_interface.hh b/eigen/bench/btl/libs/blitz/tiny_blitz_interface.hh new file mode 100644 index 0000000..6b26db7 --- /dev/null +++ b/eigen/bench/btl/libs/blitz/tiny_blitz_interface.hh @@ -0,0 +1,106 @@ +//===================================================== +// File : tiny_blitz_interface.hh +// Author : L. Plagne <laurent.plagne@edf.fr)> +// Copyright (C) EDF R&D, lun sep 30 14:23:30 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 TINY_BLITZ_INTERFACE_HH +#define TINY_BLITZ_INTERFACE_HH + +#include "blitz/array.h" +#include "blitz/tiny.h" +#include "blitz/tinymat.h" +#include "blitz/tinyvec.h" +#include <blitz/tinyvec-et.h> + +#include <vector> + +BZ_USING_NAMESPACE(blitz) + +template<class real, int SIZE> +class tiny_blitz_interface +{ + +public : + + typedef real real_type ; + + typedef std::vector<real> stl_vector; + typedef std::vector<stl_vector > stl_matrix; + + typedef TinyVector<real,SIZE> gene_vector; + typedef TinyMatrix<real,SIZE,SIZE> gene_matrix; + + static inline std::string name() { return "tiny_blitz"; } + + static void free_matrix(gene_matrix & A, int N){} + + static void free_vector(gene_vector & B){} + + static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ + 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){ + for (int i=0; i<B_stl.size() ; i++) + B(i) = B_stl[i]; + } + + static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){ + for (int i=0; i<B_stl.size() ; i++) + B_stl[i] = B(i); + } + + 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 copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){ + for (int j=0;j<N;j++) + for (int i=0;i<N;i++) + cible(i,j) = source(i,j); + } + + static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){ + for (int i=0;i<N;i++){ + cible(i) = source(i); + } + } + + static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){ + X = product(A,B); + } + + static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){ + X = product(A,B); + } + + static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){ + Y += coef * X; + } + +}; + + +#endif |