devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

imported a reworked version of BTL (Benchmark for Templated Libraries).

Browse Source 
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
This commit is contained in:
Gael Guennebaud
2008-07-09 14:04:48 +00:00
parent 5f55ab524c
commit 28539e7597
115 changed files with 8825 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

118
bench/btl/libs/ATLAS/ATLAS_LU_solve_interface.hh Normal file
View File

@@ -0,0 +1,118 @@
//=====================================================
// File : ATLAS_LU_solve_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:22 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 ATLAS_LU_solve_interface_HH
#define ATLAS_LU_solve_interface_HH
#include "ATLAS_interface.hh"
extern "C"
{
#include <atlas_level1.h>
#include <atlas_level2.h>
#include <atlas_level3.h>
#include "cblas.h"
#include <atlas_lapack.h>
}
template<class real>
class ATLAS_LU_solve_interface : public ATLAS_interface<real>
{
public :
typedef typename ATLAS_interface<real>::gene_matrix gene_matrix;
typedef typename ATLAS_interface<real>::gene_vector gene_vector;
typedef int * Pivot_Vector;
inline static void new_Pivot_Vector(Pivot_Vector & pivot, int N)
{
pivot = new int[N];
}
inline static void free_Pivot_Vector(Pivot_Vector & pivot)
{
delete pivot;
}
inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N)
{
int error=ATL_dgetrf(CblasColMajor,N,N,LU,N,pivot);
}
inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, const gene_vector &B, gene_vector X, int N)
{
copy_vector(B,X,N);
ATL_dgetrs(CblasColMajor,CblasNoTrans,N,1,LU,N,pivot,X,N);
}
};
template<>
class ATLAS_LU_solve_interface<float> : public ATLAS_interface<float>
{
public :
typedef int * Pivot_Vector;
inline static void new_Pivot_Vector(Pivot_Vector & pivot, int N)
{
pivot = new int[N];
}
inline static void free_Pivot_Vector(Pivot_Vector & pivot)
{
delete pivot;
}
inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N)
{
int error=ATL_sgetrf(CblasColMajor,N,N,LU,N,pivot);
}
inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, const gene_vector &B, gene_vector X, int N)
{
copy_vector(B,X,N);
ATL_sgetrs(CblasColMajor,CblasNoTrans,N,1,LU,N,pivot,X,N);
}
};
#endif

120
bench/btl/libs/ATLAS/ATLAS_interface.hh Normal file
View File

@@ -0,0 +1,120 @@
//=====================================================
// File : ATLAS_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:21 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 ATLAS_PRODUIT_MATRICE_VECTEUR_HH
#define ATLAS_PRODUIT_MATRICE_VECTEUR_HH
#include "f77_interface_base.hh"
#include <string>
extern "C"
{
#include <atlas_level1.h>
#include <atlas_level2.h>
#include <atlas_level3.h>
#include "cblas.h"
}
template<class real>
class ATLAS_interface : public f77_interface_base<real>
{
public :
typedef typename f77_interface_base<real>::gene_matrix gene_matrix;
typedef typename f77_interface_base<real>::gene_vector gene_vector;
static inline std::string name( void )
{
return "ATLAS";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
ATL_dgemv(CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
ATL_dgemm(CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
ATL_dgemm(CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
ATL_dgemm(CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N)
{
ATL_daxpy(N,coef,X,1,Y,1);
}
};
template<>
class ATLAS_interface<float> : public f77_interface_base<float>
{
public :
static inline std::string name( void )
{
return "ATLAS";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
ATL_sgemv(CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
ATL_sgemm(CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
ATL_sgemm(CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
ATL_sgemm(CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N)
{
ATL_saxpy(N,coef,X,1,Y,1);
}
};
#endif

4
bench/btl/libs/ATLAS/CMakeLists.txt Normal file
View File

@@ -0,0 +1,4 @@
include_directories(${BLITZ_INCLUDES})
add_executable(btl_blitz main.cpp)
target_link_libraries(btl_blitz ${BLITZ_LIBRARIES})

49
bench/btl/libs/ATLAS/main.cpp Normal file
View File

@@ -0,0 +1,49 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:21 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 "ATLAS_interface.hh"
#include "ATLAS_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"
int main()
{
bench<Action_axpy<ATLAS_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<ATLAS_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<ATLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<ATLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<ATLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
// bench<Action_lu_solve<ATLAS_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT);
return 0;
}

33
bench/btl/libs/ATLAS/titi.txt Normal file
View File

@@ -0,0 +1,33 @@
Bonjour à tous,
Une dizaine de candidats Neptune se sont déjà déclarés pour la formation C++
(sur la base de 2 jours/semaine pendant 1 mois).
Il faut faire une proposition de date pour la formation. Les vacances scolaires zone C (Paris)
se terminent le 22 avril (au matin) nous pourrions commencer ce jour.
A priori il me semble que deux jours consécutifs soient préférables. Je propose les mardi et mercredi
de chaque semaine. Sachant qu'il y a 2 jeudi (1er et 8 mai) consécutifs qui sont fériés.
Dans cette hypothèse les dates de formation seraient :
Mardi 22 avril (Vincent/Marc)
Mercredi 23 avril (Marc)
Mardi 29 avril (Marc/Antoine)
Mercredi 30 avril (Antoine)
Mardi 6 mai (Antoine)
Mercredi 7 mai (Antoine)
Mardi 13 mai (Laurent)
Mercredi 14 mai (Laurent)
J'ai mis entre parenthèse les intervenants principaux (on doit choisir le deuxième formateur pour chaque session).
Qu'en pensez-vous ?
Je dois toujours présenter un programme, pouvez-vous me donner vos programmes respectifs...?
Laurent

2
bench/btl/libs/C/CMakeLists.txt Normal file
View File

@@ -0,0 +1,2 @@
include_directories(${PROJECT_SOURCE_DIR}/libs/f77)
btl_add_bench(btl_C main.cpp)

96
bench/btl/libs/C/C_interface.hh Executable file
View File

@@ -0,0 +1,96 @@
//=====================================================
// File : C_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:23 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 C_INTERFACE_HH
#define C_INTERFACE_HH
#include "f77_interface.hh"
template<class real>
class C_interface : public f77_interface_base<real> {
public :
typedef typename f77_interface_base<real>::gene_matrix gene_matrix;
typedef typename f77_interface_base<real>::gene_vector gene_vector;
static inline std::string name() { return "C"; }
static inline void matrix_vector_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N)
{
for (int i=0;i<N;i++)
{
real somme = 0.0;
for (int j=0;j<N;j++)
somme += A[j*N+i] * B[j];
X[i] = somme;
}
}
static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N)
{
real somme;
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
somme=0.0;
for (int k=0;k<N;k++){
somme += A[i+k*N] * B[k+j*N];
}
X[i+j*N] = somme;
}
}
}
static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N)
{
real somme;
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
somme=0.0;
for (int k=0;k<N;k++){
somme+=A[k+i*N]*A[k+j*N];
}
X[i+j*N]=somme;
}
}
}
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){
real somme;
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
somme=0.0;
for (int k=0;k<N;k++){
somme+=A[i+k*N]*A[j+k*N];
}
X[i+j*N] = somme;
}
}
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){
for (int i=0;i<N;i++)
Y[i]+=coef*X[i];
}
};
#endif

44
bench/btl/libs/C/main.cpp Normal file
View File

@@ -0,0 +1,44 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:23 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 "bench.hh"
#include "C_interface.hh"
#include "action_matrix_vector_product.hh"
#include "action_matrix_matrix_product.hh"
#include "action_axpy.hh"
#include "action_ata_product.hh"
#include "action_aat_product.hh"
//#include "action_lu_solve.hh"
#include "timers/mixed_perf_analyzer.hh"
int main()
{
bench<Action_matrix_vector_product<C_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<C_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<C_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<C_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_axpy<C_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
return 0;
}

4
bench/btl/libs/C_BLAS/CMakeLists.txt Normal file
View File

@@ -0,0 +1,4 @@
include_directories(${CBLAS_INCLUDES} ${PROJECT_SOURCE_DIR}/libs/f77)
btl_add_bench(btl_cblas main.cpp)
target_link_libraries(btl_cblas ${CBLAS_LIBRARIES})

131
bench/btl/libs/C_BLAS/C_BLAS_interface.hh Normal file
View File

@@ -0,0 +1,131 @@
//=====================================================
// File : C_BLAS_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:28 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 C_BLAS_PRODUIT_MATRICE_VECTEUR_HH
#define C_BLAS_PRODUIT_MATRICE_VECTEUR_HH
#include "f77_interface.hh"
extern "C"
{
#include "cblas.h"
}
template<class real>
class C_BLAS_interface : public f77_interface_base<real>
{
public :
typedef typename f77_interface_base<real>::gene_matrix gene_matrix;
typedef typename f77_interface_base<real>::gene_vector gene_vector;
static inline std::string name( void )
{
return "C_BLAS";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
cblas_dgemv(CblasColMajor,CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
cblas_dgemv(CblasColMajor,CblasTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void transposed_matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
cblas_dgemm(CblasColMajor,CblasTrans,CblasTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N)
{
cblas_daxpy(N,coef,X,1,Y,1);
}
};
template<>
class C_BLAS_interface<float> : public f77_interface_base<float>
{
public :
static inline std::string name( void )
{
return "C_BLAS";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
cblas_sgemv(CblasColMajor,CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
cblas_sgemv(CblasColMajor,CblasTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void transposed_matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
cblas_sgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N)
{
cblas_saxpy(N,coef,X,1,Y,1);
}
};
#endif

49
bench/btl/libs/C_BLAS/main.cpp Normal file
View File

@@ -0,0 +1,49 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:28 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 "C_BLAS_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"
int main()
{
bench<Action_axpy<C_BLAS_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<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);
//bench<Action_lu_solve<C_BLAS_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT);
return 0;
}

127
bench/btl/libs/INTEL_BLAS/INTEL_BLAS_LU_solve_interface.hh Normal file
View File

@@ -0,0 +1,127 @@
//=====================================================
// File : INTEL_BLAS_LU_solve_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:29 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 INTEL_BLAS_LU_solve_interface_HH
#define INTEL_BLAS_LU_solve_interface_HH
#include "INTEL_BLAS_interface.hh"
extern "C"
{
// void dgetrf_(int *M, int *N, double *A, int *LDA, int *IPIV, int *INFO);
// void dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO);
// void sgetrf_(int *M, int *N, float *A, int *LDA, int *IPIV, int *INFO);
// void sgetrs_(char *TRANS, int *N, int *NRHS, float *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO);
#include "mkl_lapack.h"
}
template<class real>
class INTEL_BLAS_LU_solve_interface : public INTEL_BLAS_interface<real>
{
public :
typedef typename INTEL_BLAS_interface<real>::gene_matrix gene_matrix;
typedef typename INTEL_BLAS_interface<real>::gene_vector gene_vector;
typedef int * Pivot_Vector;
inline static void new_Pivot_Vector(Pivot_Vector & pivot, int N)
{
pivot = new int[N];
}
inline static void free_Pivot_Vector(Pivot_Vector & pivot)
{
delete pivot;
}
inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N)
{
int info;
DGETRF(&N,&N,LU,&N,pivot,&info);
}
inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, const gene_vector &B, gene_vector X, int N)
{
int info;
int one=1;
char * transpose="N";
copy_vector(B,X,N);
DGETRS(transpose,&N,&one,LU,&N,pivot,X,&N,&info);
}
};
template<>
class INTEL_BLAS_LU_solve_interface<float> : public INTEL_BLAS_interface<float>
{
public :
typedef int * Pivot_Vector;
inline static void new_Pivot_Vector(Pivot_Vector & pivot, int N)
{
pivot = new int[N];
}
inline static void free_Pivot_Vector(Pivot_Vector & pivot)
{
delete pivot;
}
inline static void LU_factor(gene_matrix & LU, Pivot_Vector & pivot, int N)
{
int info;
SGETRF(&N,&N,LU,&N,pivot,&info);
}
inline static void LU_solve(const gene_matrix & LU, const Pivot_Vector pivot, const gene_vector &B, gene_vector X, int N)
{
char * transpose="N";
int info;
int one=1;
copy_vector(B,X,N);
SGETRS(transpose,&N,&one,LU,&N,pivot,X,&N,&info);
}
};
#endif

95
bench/btl/libs/INTEL_BLAS/INTEL_BLAS_interface.hh Normal file
View File

@@ -0,0 +1,95 @@
//=====================================================
// File : INTEL_BLAS_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:29 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 INTEL_BLAS_PRODUIT_MATRICE_VECTEUR_HH
#define INTEL_BLAS_PRODUIT_MATRICE_VECTEUR_HH
#include "f77_interface.hh"
extern "C"
{
#include "mkl_cblas.h"
}
template<class real>
class INTEL_BLAS_interface : public f77_interface_base<real>
{
public :
typedef typename f77_interface_base<real>::gene_matrix gene_matrix;
typedef typename f77_interface_base<real>::gene_vector gene_vector;
static inline std::string name( void ) { return "INTEL_BLAS"; }
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) {
cblas_dgemv(CblasColMajor,CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N) {
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N) {
cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N) {
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N) {
cblas_daxpy(N,coef,X,1,Y,1);
}
};
template<>
class INTEL_BLAS_interface<float> : public f77_interface_base<float>
{
public :
static inline std::string name() { return "INTEL_BLAS"; }
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N) {
// cblas_sgemv(CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
cblas_sgemv(CblasColMajor,CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N) {
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N) {
cblas_sgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N) {
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
}
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N) {
cblas_saxpy(N,coef,X,1,Y,1);
}
};
#endif

2
bench/btl/libs/INTEL_BLAS/config.sh Executable file
View File

@@ -0,0 +1,2 @@
#! /bin/bash
export LD_LIBRARY_PATH=/opt/intel/mkl/lib/32:${LD_LIBRARY_PATH}

49
bench/btl/libs/INTEL_BLAS/main.cpp Normal file
View File

@@ -0,0 +1,49 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:29 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 "INTEL_BLAS_interface.hh"
#include "INTEL_BLAS_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"
int main()
{
bench<Action_axpy<INTEL_BLAS_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<INTEL_BLAS_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<INTEL_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<INTEL_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<INTEL_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
// bench<Action_lu_solve<INTEL_BLAS_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT);
return 0;
}

2
bench/btl/libs/STL/CMakeLists.txt Normal file
View File

@@ -0,0 +1,2 @@
btl_add_bench(btl_STL main.cpp)

190
bench/btl/libs/STL/STL_interface.hh Normal file
View File

@@ -0,0 +1,190 @@
//=====================================================
// File : STL_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:24 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 STL_INTERFACE_HH
#define STL_INTERFACE_HH
#include <string>
#include <vector>
#include "utilities.h"
using namespace std;
template<class real>
class STL_interface{
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector > stl_matrix;
typedef stl_matrix gene_matrix;
typedef stl_vector gene_vector;
static inline std::string name( void )
{
return "STL";
}
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 = A_stl;
}
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){
A_stl = A ;
}
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 copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
for (int i=0;i<N;i++)
for (int j=0;j<N;j++)
cible[i][j]=source[i][j];
}
static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N)
{
real somme;
for (int j=0;j<N;j++){
for (int i=0;i<N;i++){
somme=0.0;
for (int k=0;k<N;k++)
somme += A[i][k]*A[j][k];
X[j][i]=somme;
}
}
}
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N)
{
real somme;
for (int j=0;j<N;j++){
for (int i=0;i<N;i++){
somme=0.0;
for (int k=0;k<N;k++){
somme+=A[k][i]*A[k][j];
}
X[j][i]=somme;
}
}
}
static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N)
{
real somme;
for (int j=0;j<N;j++){
for (int i=0;i<N;i++){
somme=0.0;
for (int k=0;k<N;k++)
somme+=A[k][i]*B[j][k];
X[j][i]=somme;
}
}
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
real somme;
for (int i=0;i<N;i++){
somme=0.0;
for (int j=0;j<N;j++)
somme+=A[j][i]*B[j];
X[i]=somme;
}
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
real somme;
for (int i=0;i<N;i++){
somme = 0.0;
for (int j=0;j<N;j++)
somme += A[i][j]*B[j];
X[i] = somme;
}
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N){
for (int i=0;i<N;i++)
Y[i]+=coef*X[i];
}
static inline real norm_diff(const stl_vector & A, const stl_vector & B)
{
int N=A.size();
real somme=0.0;
real somme2=0.0;
for (int i=0;i<N;i++){
real diff=A[i]-B[i];
somme+=diff*diff;
somme2+=A[i]*A[i];
}
return somme/somme2;
}
static inline real norm_diff(const stl_matrix & A, const stl_matrix & B)
{
int N=A[0].size();
real somme=0.0;
real somme2=0.0;
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
real diff=A[i][j] - B[i][j];
somme += diff*diff;
somme2 += A[i][j]*A[i][j];
}
}
return somme/somme2;
}
static inline void display_vector(const stl_vector & A)
{
int N=A.size();
for (int i=0;i<N;i++){
INFOS("A["<<i<<"]="<<A[i]<<endl);
}
}
};
#endif

43
bench/btl/libs/STL/main.cpp Normal file
View File

@@ -0,0 +1,43 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:23 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 "STL_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<STL_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<STL_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_atv_product<STL_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<STL_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<STL_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<STL_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}

2
bench/btl/libs/STL_algo/CMakeLists.txt Normal file
View File

@@ -0,0 +1,2 @@
btl_add_bench(btl_STL_algo main.cpp)

138
bench/btl/libs/STL_algo/STL_algo_interface.hh Normal file
View File

@@ -0,0 +1,138 @@
//=====================================================
// File : STL_algo_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:24 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 STL_ALGO_INTERFACE_HH
#define STL_ALGO_INTERFACE_HH
#include <string>
#include <vector>
#include <numeric>
#include <algorithm>
#include "utilities.h"
template<class real>
class STL_algo_interface{
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector > stl_matrix;
typedef stl_matrix gene_matrix;
typedef stl_vector gene_vector;
static inline std::string name( void )
{
return "STL_algo";
}
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=A_stl ;
}
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){
A_stl=A ;
}
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 copy_matrix(const gene_matrix & source, gene_matrix & cible, int N)
{
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
cible[i][j]=source[i][j];
}
}
}
class somme {
public:
somme(real coef):_coef(coef){};
real operator()(const real & val1, const real & val2)
{
return _coef * val1 + val2;
}
private:
real _coef;
};
class vector_generator {
public:
vector_generator(const gene_matrix & a_matrix, const gene_vector & a_vector):
_matrice(a_matrix),
_vecteur(a_vector),
_index(0)
{};
real operator()( void )
{
const gene_vector & ai=_matrice[_index];
int N=ai.size();
_index++;
return std::inner_product(&ai[0],&ai[N],&_vecteur[0],0.0);
}
private:
int _index;
const gene_matrix & _matrice;
const gene_vector & _vecteur;
};
static inline void atv_product(const gene_matrix & A, const gene_vector & B, gene_vector & X, int N)
{
std::generate(&X[0],&X[N],vector_generator(A,B));
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N)
{
std::transform(&X[0],&X[N],&Y[0],&Y[0],somme(coef));
}
};
#endif

37
bench/btl/libs/STL_algo/main.cpp Normal file
View File

@@ -0,0 +1,37 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:23 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 "STL_algo_interface.hh"
#include "bench.hh"
#include "action_atv_product.hh"
#include "action_axpy.hh"
int main()
{
bench<Action_atv_product<STL_algo_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_axpy<STL_algo_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
return 0;
}

4
bench/btl/libs/blitz/CMakeLists.txt Normal file
View File

@@ -0,0 +1,4 @@
include_directories(${BLITZ_INCLUDES})
btl_add_bench(btl_blitz main.cpp)
target_link_libraries(btl_blitz ${BLITZ_LIBRARIES})

192
bench/btl/libs/blitz/blitz_LU_solve_interface.hh Normal file
View File

@@ -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

147
bench/btl/libs/blitz/blitz_interface.hh Normal file
View File

@@ -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 <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 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

49
bench/btl/libs/blitz/main.cpp Normal file
View File

@@ -0,0 +1,49 @@
//=====================================================
// 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"
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;
}

8
bench/btl/libs/eigen2/CMakeLists.txt Normal file
View File

@@ -0,0 +1,8 @@
include_directories(${EIGEN2_INCLUDE_DIR})
btl_add_bench(btl_eigen2 main.cpp)
IF(NOT BTL_NOVEC)
btl_add_bench(btl_eigen2_novec main.cpp)
set_target_properties(btl_eigen2_novec PROPERTIES COMPILE_FLAGS "-DEIGEN_DONT_VECTORIZE")
ENDIF(NOT BTL_NOVEC)

192
bench/btl/libs/eigen2/eigen2_LU_solve_interface.hh Normal file
View File

@@ -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

128
bench/btl/libs/eigen2/eigen2_interface.hh Normal file
View File

@@ -0,0 +1,128 @@
//=====================================================
// 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 EIGEN2_INTERFACE_HH
#define EIGEN2_INTERFACE_HH
#include <Eigen/Core>
#include <vector>
using namespace Eigen;
template<class real, int SIZE=Dynamic>
class eigen2_interface
{
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector> stl_matrix;
typedef Eigen::Matrix<real,SIZE,SIZE> gene_matrix;
typedef Eigen::Matrix<real,SIZE,1> gene_vector;
static inline std::string name( void )
{
#if defined(EIGEN_VECTORIZE_SSE)
if (SIZE==Dynamic) return "eigen2_SSE"; else return "tiny_eigen2_SSE";
#elif defined(EIGEN_VECTORIZE_ALTIVEC)
if (SIZE==Dynamic) return "eigen2_AltiVec"; else return "tiny_eigen2_AltiVec";
#else
if (SIZE==Dynamic) return "eigen2"; else return "tiny_eigen2";
#endif
}
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.coeffRef(i,j) = A_stl[j][i];
}
}
}
static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){
B.resize(B_stl.size(),1);
for (int i=0; i<B_stl.size() ; i++){
B.coeffRef(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.coeff(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.coeff(i,j);
}
}
}
static inline void matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){
X = (A*B).lazy();
}
static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){
X = (A.transpose()*B.transpose()).lazy();
}
static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){
X = (A.transpose()*A).lazy();
}
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){
X = (A*A.transpose()).lazy();
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X = (A*B).lazy();
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X = (A*B).lazy();
}
static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){
Y += coef * X;
}
static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
cible = source;
}
static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){
cible = source;
}
};
#endif

49
bench/btl/libs/eigen2/main.cpp Normal file
View File

@@ -0,0 +1,49 @@
//=====================================================
// 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 "eigen2_interface.hh"
#include "bench.hh"
#include "static/bench_static.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_matrix_vector_product<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_atv_product<eigen2_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_axpy<eigen2_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_matrix_product<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<eigen2_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);
bench_static<Action_axpy,eigen2_interface>();
bench_static<Action_matrix_matrix_product,eigen2_interface>();
bench_static<Action_matrix_vector_product,eigen2_interface>();
return 0;
}

3
bench/btl/libs/f77/CMakeLists.txt Normal file
View File

@@ -0,0 +1,3 @@
enable_language(Fortran)
btl_add_bench(btl_f77 main.cpp dmxv.f smxv.f dmxm.f smxm.f daxpy.f saxpy.f data.f sata.f daat.f saat.f)

14
bench/btl/libs/f77/daat.f Normal file
View File

@@ -0,0 +1,14 @@
SUBROUTINE DAAT(A,X,N)
**
** X = AT * A
REAL*8 A(N,N),X(N,N),R
DO 20 I=1,N
DO 20 J=1,N
R=0.
DO 10 K=1,N
R=R+A(I,K)*A(J,K)
10 CONTINUE
X(I,J)=R
20 CONTINUE
RETURN
END

14
bench/btl/libs/f77/data.f Normal file
View File

@@ -0,0 +1,14 @@
SUBROUTINE DATA(A,X,N)
**
** X = AT * A
REAL*8 A(N,N),X(N,N),R
DO 20 I=1,N
DO 20 J=1,N
R=0.
DO 10 K=1,N
R=R+A(K,I)*A(K,J)
10 CONTINUE
X(I,J)=R
20 CONTINUE
RETURN
END

18
bench/btl/libs/f77/daxpy.f Normal file
View File

@@ -0,0 +1,18 @@
SUBROUTINE DAXPY(N,A,X,Y)
** ***************************************
** CALCULE Y = Y + A*X
** ***************************************
*>N NOMBRE D'OPERATIONS A FAIRE
*>A CONSTANTE MULTIPLICATIVE
*>X TABLEAU
*=Y TABLEAU DES RESULTATS
*A R. SANCHEZ ( EARLY WINTER 1987 )
*V M.F. ROBEAU
REAL*8 X(1),Y(1)
REAL*8 A
DO 10 I=1,N
Y(I)=Y(I)+A*X(I)
10 CONTINUE
RETURN
END

32
bench/btl/libs/f77/dmxm.f Normal file
View File

@@ -0,0 +1,32 @@
SUBROUTINE DMXM(A,N,B,M,C,L)
**
** C = A * B
** A ET B MATRICES A(N,M) B(M,L) ==> C(N,L)
**
*>A PREMIERE MATRICE
*>N PREMIERE DIMENSION DE A ET DE C
*>B DEUXIEME MATRICE
*>M DEUXIEME DIMENSION DE A ET PERMIERE DE B
*<C MATRICE PRODUIT DE A ET DE B
*>L DEUXIEME DIMENSION DE B ET DE C
*A R. SANCHEZ ( EARLY WINTER 1987 )
*V M.F. ROBEAU
*M AM BAUDRON - AVRIL 94
*: ERREUR DANS L'APPEL A L'UTILITAIRE SGEMM
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR HP
*M AM BAUDRON - NOVEMBRE 1991
*: ERREUR ( SOMME SUR LES TERMES PAS FAITE )
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR RISC
*M AM BAUDRON - MAI 1993
*: CHANGEMENT DES %IF LOCAL SUN MIPS SUITE A INTRODUCTION VERSION IBM
REAL*8 A(N,M),B(M,L),C(N,L),R
DO 20 I=1,N
DO 20 J=1,L
R=0.
DO 10 K=1,M
R=R+A(I,K)*B(K,J)
10 CONTINUE
C(I,J)=R
20 CONTINUE
RETURN
END

36
bench/btl/libs/f77/dmxm.f.mfr Normal file
View File

@@ -0,0 +1,36 @@
SUBROUTINE DMXM(A,N,B,M,C,L)
**
** C = A * B
** A ET B MATRICES A(N,M) B(M,L) ==> C(N,L)
**
*>A PREMIERE MATRICE
*>N PREMIERE DIMENSION DE A ET DE C
*>B DEUXIEME MATRICE
*>M DEUXIEME DIMENSION DE A ET PERMIERE DE B
*<C MATRICE PRODUIT DE A ET DE B
*>L DEUXIEME DIMENSION DE B ET DE C
*A R. SANCHEZ ( EARLY WINTER 1987 )
*V M.F. ROBEAU
*M AM BAUDRON - AVRIL 94
*: ERREUR DANS L'APPEL A L'UTILITAIRE SGEMM
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR HP
*M AM BAUDRON - NOVEMBRE 1991
*: ERREUR ( SOMME SUR LES TERMES PAS FAITE )
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR RISC
*M AM BAUDRON - MAI 1993
*: CHANGEMENT DES %IF LOCAL SUN MIPS SUITE A INTRODUCTION VERSION IBM
REAL*8 A(N,M),B(M,L),C(N,L),R
DO 5 J=1,L
DO 5 I=1,N
C(I,J)=0.
5 CONTINUE
DO 10 K=1,M
DO 20 J=1,L
R=B(K,J)
DO 20 I=1,N
C(I,J)=C(I,J)+A(I,K)*R
20 CONTINUE
10 CONTINUE
RETURN
END

39
bench/btl/libs/f77/dmxv.f Normal file
View File

@@ -0,0 +1,39 @@
SUBROUTINE DMXV(A,N,X,M,R)
C
**
** VERSION DOUBLE PRECISION DE MXV
** R = A * X
** A MATRICE A(N,M)
** R ET X VECTEURS
**
*>A PREMIERE MATRICE
*>N PREMIERE DIMENSION DE A
*>X VECTEUR
*>M DEUXIEME DIMENSION DE A
*<R VECTEUR PRODUIT DE A ET DE X
**
*A M. COSTE
*V M.F. ROBEAU
*M
*
REAL*8 X(1),R(1),A(N,M)
REAL*8 S
C DO 20 I=1,N
C S=0.
C DO 10 J=1,M
C S=S+A(I,J)*X(J)
C 10 CONTINUE
C R(I)=S
C 20 CONTINUE
DO 5 I=1,N
R(I)=0
5 CONTINUE
DO 10 J=1,M
S=X(J)
DO 20 I=1,N
R(I)=R(I)+A(I,J)*S
20 CONTINUE
10 CONTINUE
RETURN
END

145
bench/btl/libs/f77/f77_interface.hh Normal file
View File

@@ -0,0 +1,145 @@
//=====================================================
// File : f77_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:24 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 F77_INTERFACE_HH
#define F77_INTERFACE_HH
#include "f77_interface_base.hh"
#include <string>
extern "C"
{
void dmxv_(double * A, int * N, double * X, int * M, double *R);
void smxv_(float * A, int * N, float * X, int * M, float *R);
void dmxm_(double * A, int * N, double * B, int * M, double *C, int * K);
void smxm_(float * A, int * N, float * B, int * M, float *C, int * K);
void data_(double * A, double *X, int * N);
void sata_(float * A, float *X, int * N);
void daat_(double * A, double *X, int * N);
void saat_(float * A, float *X, int * N);
void saxpy_(int * N, float * coef, float * X, float *Y);
void daxpy_(int * N, double * coef, double * X, double *Y);
}
template<class real>
class f77_interface : public f77_interface_base<real>
{
public :
typedef typename f77_interface_base<real>::gene_matrix gene_matrix;
typedef typename f77_interface_base<real>::gene_vector gene_vector;
static inline std::string name( void )
{
return "f77";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
dmxv_(A,&N,B,&N,X);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
dmxm_(A,&N,B,&N,X,&N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
data_(A,X,&N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
daat_(A,X,&N);
}
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N)
{
int one=1;
daxpy_(&N,&coef,X,Y);
}
};
template<>
class f77_interface<float> : public f77_interface_base<float>
{
public :
static inline std::string name( void )
{
return "F77";
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N)
{
smxv_(A,&N,B,&N,X);
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N)
{
smxm_(A,&N,B,&N,X,&N);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N)
{
sata_(A,X,&N);
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N)
{
saat_(A,X,&N);
}
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N)
{
saxpy_(&N,&coef,X,Y);
}
};
#endif

91
bench/btl/libs/f77/f77_interface_base.hh Normal file
View File

@@ -0,0 +1,91 @@
//=====================================================
// File : f77_interface_base.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:25 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 F77_INTERFACE_BASE_HH
#define F77_INTERFACE_BASE_HH
#include "utilities.h"
#include <vector>
template<class real>
class f77_interface_base{
public:
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector > stl_matrix;
typedef real * gene_matrix;
typedef real * gene_vector;
static void free_matrix(gene_matrix & A, int N){
delete A;
}
static void free_vector(gene_vector & B){
delete B;
}
static inline void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){
int N = A_stl.size();
A = new real[N*N];
for (int j=0;j<N;j++)
for (int i=0;i<N;i++)
A[i+N*j] = A_stl[j][i];
}
static inline void vector_from_stl(gene_vector & B, stl_vector & B_stl){
int N = B_stl.size();
B = new real[N];
for (int i=0;i<N;i++)
B[i] = B_stl[i];
}
static inline void vector_to_stl(gene_vector & B, stl_vector & B_stl){
int N = B_stl.size();
for (int i=0;i<N;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+N*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 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+N*j] = source[i+N*j];
}
}
}
};
#endif

46
bench/btl/libs/f77/main.cpp Normal file
View File

@@ -0,0 +1,46 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:25 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 "f77_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"
int main()
{
bench<Action_axpy<f77_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<f77_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<f77_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<f77_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<f77_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}

14
bench/btl/libs/f77/saat.f Normal file
View File

@@ -0,0 +1,14 @@
SUBROUTINE SAAT(A,X,N)
**
** X = AT * A
REAL*4 A(N,N),X(N,N)
DO 20 I=1,N
DO 20 J=1,N
R=0.
DO 10 K=1,N
R=R+A(I,K)*A(J,K)
10 CONTINUE
X(I,J)=R
20 CONTINUE
RETURN
END

14
bench/btl/libs/f77/sata.f Normal file
View File

@@ -0,0 +1,14 @@
SUBROUTINE SATA(A,X,N)
**
** X = AT * A
REAL*4 A(N,N),X(N,N)
DO 20 I=1,N
DO 20 J=1,N
R=0.
DO 10 K=1,N
R=R+A(K,I)*A(K,J)
10 CONTINUE
X(I,J)=R
20 CONTINUE
RETURN
END

16
bench/btl/libs/f77/saxpy.f Normal file
View File

@@ -0,0 +1,16 @@
SUBROUTINE SAXPY(N,A,X,Y)
** ***************************************
** CALCULE Y = Y + A*X
** ***************************************
*>N NOMBRE D'OPERATIONS A FAIRE
*>A CONSTANTE MULTIPLICATIVE
*>X TABLEAU
*=Y TABLEAU DES RESULTATS
*A R. SANCHEZ ( EARLY WINTER 1987 )
*V M.F. ROBEAU
DIMENSION X(1),Y(1)
DO 10 I=1,N
Y(I)=Y(I)+A*X(I)
10 CONTINUE
RETURN
END

32
bench/btl/libs/f77/smxm.f Normal file
View File

@@ -0,0 +1,32 @@
SUBROUTINE SMXM(A,N,B,M,C,L)
**
** C = A * B
** A ET B MATRICES A(N,M) B(M,L) ==> C(N,L)
**
*>A PREMIERE MATRICE
*>N PREMIERE DIMENSION DE A ET DE C
*>B DEUXIEME MATRICE
*>M DEUXIEME DIMENSION DE A ET PERMIERE DE B
*<C MATRICE PRODUIT DE A ET DE B
*>L DEUXIEME DIMENSION DE B ET DE C
*A R. SANCHEZ ( EARLY WINTER 1987 )
*V M.F. ROBEAU
*M AM BAUDRON - AVRIL 94
*: ERREUR DANS L'APPEL A L'UTILITAIRE SGEMM
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR HP
*M AM BAUDRON - NOVEMBRE 1991
*: ERREUR ( SOMME SUR LES TERMES PAS FAITE )
*: APPEL A L'UTILITAIRE SGEMM DE LA LIBRAIRIE BLAS SUR RISC
*M AM BAUDRON - MAI 1993
*: CHANGEMENT DES %IF LOCAL SUN MIPS SUITE A INTRODUCTION VERSION IBM
DIMENSION A(N,M),B(M,L),C(N,L)
DO 20 I=1,N
DO 20 J=1,L
R=0.
DO 10 K=1,M
R=R+A(I,K)*B(K,J)
10 CONTINUE
C(I,J)=R
20 CONTINUE
RETURN
END

38
bench/btl/libs/f77/smxv.f Normal file
View File

@@ -0,0 +1,38 @@
SUBROUTINE SMXV(A,N,X,M,R)
C
**
** VERSION DOUBLE PRECISION DE MXV
** R = A * X
** A MATRICE A(N,M)
** R ET X VECTEURS
**
*>A PREMIERE MATRICE
*>N PREMIERE DIMENSION DE A
*>X VECTEUR
*>M DEUXIEME DIMENSION DE A
*<R VECTEUR PRODUIT DE A ET DE X
**
*A M. COSTE
*V M.F. ROBEAU
*M
*
REAL*4 X(1),R(1),A(N,M)
REAL*4 S
C DO 20 I=1,N
C S=0.
C DO 10 J=1,M
C S=S+A(I,J)*X(J)
C 10 CONTINUE
C R(I)=S
C 20 CONTINUE
DO 5 I=1,N
R(I)=0
5 CONTINUE
DO 10 J=1,M
S=X(J)
DO 20 I=1,N
R(I)=R(I)+A(I,J)*S
20 CONTINUE
10 CONTINUE
RETURN
END

36
bench/btl/libs/f77/test_interface.hh Normal file
View File

@@ -0,0 +1,36 @@
//=====================================================
// File : test_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:25 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 TEST_INTERFACE_HH
#define TEST_INTERFACE_HH
template< class Interface >
void test_interface( void ){
Interface::interface_name();
typename Interface::gene_matrix A;
}
#endif

3
bench/btl/libs/gmm/CMakeLists.txt Normal file
View File

@@ -0,0 +1,3 @@
include_directories(${GMM_INCLUDES})
btl_add_bench(btl_gmm main.cpp)

192
bench/btl/libs/gmm/gmm_LU_solve_interface.hh Normal file
View File

@@ -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

119
bench/btl/libs/gmm/gmm_interface.hh Normal file
View File

@@ -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

44
bench/btl/libs/gmm/main.cpp Normal file
View File

@@ -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;
}

12
bench/btl/libs/mtl4/.kdbgrc.main Normal file
View File

@@ -0,0 +1,12 @@
[General]
DebuggerCmdStr=
DriverName=GDB
FileVersion=1
OptionsSelected=
ProgramArgs=
TTYLevel=7
WorkingDirectory=
[Memory]
ColumnWidths=80,0
NumExprs=0

42
bench/btl/libs/mtl4/main.cpp Normal file
View File

@@ -0,0 +1,42 @@
//=====================================================
// 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 "mtl4_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<mtl4_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<mtl4_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_atv_product<mtl4_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<mtl4_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<mtl4_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<mtl4_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}

192
bench/btl/libs/mtl4/mtl4_LU_solve_interface.hh Normal file
View File

@@ -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

120
bench/btl/libs/mtl4/mtl4_interface.hh Normal file
View File

@@ -0,0 +1,120 @@
//=====================================================
// 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 MTL4_INTERFACE_HH
#define MTL4_INTERFACE_HH
#include <boost/numeric/mtl/mtl.hpp>
#include <vector>
using namespace mtl;
template<class real>
class mtl4_interface {
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector > stl_matrix;
typedef mtl::dense2D<real, mtl::matrix::parameters<mtl::tag::col_major> > gene_matrix;
typedef mtl::dense_vector<real> gene_vector;
static inline std::string name() { return "mtl4"; }
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.change_dim(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.change_dim(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){
X = (A*B);
}
static inline void transposed_matrix_matrix_product(const gene_matrix & A, const gene_matrix & B, gene_matrix & X, int N){
X = (trans(A)*trans(B));
}
static inline void ata_product(const gene_matrix & A, gene_matrix & X, int N){
X = (trans(A)*A);
}
static inline void aat_product(const gene_matrix & A, gene_matrix & X, int N){
X = (A*trans(A));
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X = (A*B);
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X = (trans(A)*B);
}
static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){
Y += coef * X;
}
static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
cible = source;
}
static inline void copy_vector(const gene_vector & source, gene_vector & cible, int N){
cible = source;
}
};
#endif

4
bench/btl/libs/tiny_blitz/CMakeLists.txt Normal file
View File

@@ -0,0 +1,4 @@
include_directories(${BLITZ_INCLUDES})
btl_add_bench(btl_tiny_blitz main.cpp)
target_link_libraries(btl_tiny_blitz ${BLITZ_LIBRARIES})

37
bench/btl/libs/tiny_blitz/main.cpp Normal file
View File

@@ -0,0 +1,37 @@
//=====================================================
// 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"
#include "timers/x86_perf_analyzer.hh"
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;
}

106
bench/btl/libs/tiny_blitz/tiny_blitz_interface.hh Normal file
View File

@@ -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

4
bench/btl/libs/tvmet/CMakeLists.txt Normal file
View File

@@ -0,0 +1,4 @@
include_directories(${BLITZ_INCLUDES})
add_executable(btl_blitz main.cpp)
target_link_libraries(btl_blitz ${BLITZ_LIBRARIES})

36
bench/btl/libs/tvmet/main.cpp Normal file
View File

@@ -0,0 +1,36 @@
//=====================================================
// 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 "tvmet_interface.hh"
#include "static/bench_static.hh"
#include "action_matrix_vector_product.hh"
#include "action_matrix_matrix_product.hh"
#include "action_axpy.hh"
int main()
{
bench_static<Action_axpy,tvmet_interface>();
bench_static<Action_matrix_matrix_product,tvmet_interface>();
bench_static<Action_matrix_vector_product,tvmet_interface>();
return 0;
}

141
bench/btl/libs/tvmet/tvmet_interface.hh Normal file
View File

@@ -0,0 +1,141 @@
//=====================================================
// File : tvmet_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 TVMET_INTERFACE_HH
#define TVMET_INTERFACE_HH
#include <tvmet/Vector.h>
#include <tvmet/Matrix.h>
#include <vector>
using namespace tvmet;
template<class real, int SIZE>
class tvmet_interface{
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector > stl_matrix;
typedef Vector<real,SIZE> gene_vector;
typedef Matrix<real,SIZE,SIZE> gene_matrix;
static inline std::string name( void )
{
return "tvmet";
}
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){
for (int i=0; i<A_stl.size() ; i++){
for (int j=0; j<A_stl[i].size() ; j++){
A(i,j)=A_stl[i][j];
}
}
}
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 i=0;i<N;i++){
A_stl[i].resize(N);
for (int j=0;j<N;j++){
A_stl[i][j]=A(i,j);
}
}
}
static inline void copy_matrix(const gene_matrix & source, gene_matrix & cible, int N)
{
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)
{
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

3
bench/btl/libs/ublas/CMakeLists.txt Normal file
View File

@@ -0,0 +1,3 @@
include_directories(${Boost_INCLUDES})
btl_add_bench(btl_ublas main.cpp)

44
bench/btl/libs/ublas/main.cpp Normal file
View File

@@ -0,0 +1,44 @@
//=====================================================
// File : main.cpp
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:27 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 "ublas_interface.hh"
#include "bench.hh"
#include "action_matrix_vector_product.hh"
#include "action_matrix_matrix_product.hh"
#include "action_axpy.hh"
#include "action_ata_product.hh"
#include "action_aat_product.hh"
#include "action_atv_product.hh"
int main()
{
bench<Action_axpy<ublas_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT);
bench<Action_matrix_vector_product<ublas_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_atv_product<ublas_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,NB_POINT);
bench<Action_matrix_matrix_product<ublas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_ata_product<ublas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_aat_product<ublas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}

131
bench/btl/libs/ublas/ublas_interface.hh Normal file
View File

@@ -0,0 +1,131 @@
//=====================================================
// File : ublas_interface.hh
// Author : L. Plagne <laurent.plagne@edf.fr)>
// Copyright (C) EDF R&D, lun sep 30 14:23:27 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 UBLAS_INTERFACE_HH
#define UBLAS_INTERFACE_HH
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
template <class real>
class ublas_interface{
public :
typedef real real_type ;
typedef std::vector<real> stl_vector;
typedef std::vector<stl_vector> stl_matrix;
typedef typename boost::numeric::ublas::matrix<real,boost::numeric::ublas::column_major> gene_matrix;
typedef typename boost::numeric::ublas::vector<real> gene_vector;
static inline std::string name( void ) { return "ublas"; }
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.size(),A_stl[0].size());
for (int i=0; i<A_stl.size() ; i++)
for (int j=0; j<A_stl[i].size() ; j++)
A(i,j)=A_stl[i][j];
}
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 i=0;i<N;i++)
{
A_stl[i].resize(N);
for (int j=0;j<N;j++)
A_stl[i][j]=A(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 copy_matrix(const gene_matrix & source, gene_matrix & cible, int N){
for (int i=0;i<N;i++){
for (int j=0;j<N;j++){
cible(i,j) = source(i,j);
}
}
}
static inline void matrix_vector_product_slow(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X = prod(A,B);
}
static inline void matrix_matrix_product_slow(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
X = prod(A,B);
}
static inline void axpy_slow(const real coef, const gene_vector & X, gene_vector & Y, int N){
Y+=coef*X;
}
// alias free assignements
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X.assign(prod(A,B));
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
X.assign(prod(trans(A),B));
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
X.assign(prod(A,B));
}
static inline void axpy(const real coef, const gene_vector & X, gene_vector & Y, int N){
Y.plus_assign(coef*X);
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N){
// X = prod(trans(A),A);
X.assign(prod(trans(A),A));
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N){
// X = prod(A,trans(A));
X.assign(prod(A,trans(A)));
}
};
#endif