mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

BTL: - added tridiagonalization and hessenberg decomposition bench

Browse Source 
- added GOTO library
This commit is contained in:
Gael Guennebaud 2008-07-28 20:48:21 +00:00
parent 44d95e0540
commit e0215ee510
21 changed files with 421 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

230
bench/btl/actions/action_hessenberg.hh Normal file
View File

@ -0,0 +1,230 @@
//=====================================================
// File : action_hessenberg.hh
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
//=====================================================
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
#ifndef ACTION_HESSENBERG
#define ACTION_HESSENBERG
#include "utilities.h"
#include "STL_interface.hh"
#include <string>
#include "init/init_function.hh"
#include "init/init_vector.hh"
#include "init/init_matrix.hh"
using namespace std;
template<class Interface>
class Action_hessenberg {
public :
// Ctor
Action_hessenberg( int size ):_size(size)
{
MESSAGE("Action_hessenberg Ctor");
// STL vector initialization
init_matrix<pseudo_random>(X_stl,_size);
init_matrix<null_function>(C_stl,_size);
init_matrix<null_function>(resu_stl,_size);
// generic matrix and vector initialization
Interface::matrix_from_stl(X_ref,X_stl);
Interface::matrix_from_stl(X,X_stl);
Interface::matrix_from_stl(C,C_stl);
_cost = 0;
for (int j=0; j<_size-2; ++j)
{
int r = std::max(0,_size-j-1);
int b = std::max(0,_size-j-2);
_cost += 6 + 3*b + r*r*4 + r*_size*4;
}
}
// invalidate copy ctor
Action_hessenberg( const Action_hessenberg & )
{
INFOS("illegal call to Action_hessenberg Copy Ctor");
exit(1);
}
// Dtor
~Action_hessenberg( void ){
MESSAGE("Action_hessenberg Dtor");
// deallocation
Interface::free_matrix(X_ref,_size);
Interface::free_matrix(X,_size);
Interface::free_matrix(C,_size);
}
// action name
static inline std::string name( void )
{
return "hessenberg_"+Interface::name();
}
double nb_op_base( void ){
return _cost;
}
inline void initialize( void ){
Interface::copy_matrix(X_ref,X,_size);
}
inline void calculate( void ) {
Interface::hessenberg(X,C,_size);
}
void check_result( void ){
// calculation check
Interface::matrix_to_stl(C,resu_stl);
// STL_interface<typename Interface::real_type>::hessenberg(X_stl,C_stl,_size);
//
// typename Interface::real_type error=
// STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
//
// if (error>1.e-6){
// INFOS("WRONG CALCULATION...residual=" << error);
// exit(0);
// }
}
private :
typename Interface::stl_matrix X_stl;
typename Interface::stl_matrix C_stl;
typename Interface::stl_matrix resu_stl;
typename Interface::gene_matrix X_ref;
typename Interface::gene_matrix X;
typename Interface::gene_matrix C;
int _size;
double _cost;
};
template<class Interface>
class Action_tridiagonalization {
public :
// Ctor
Action_tridiagonalization( int size ):_size(size)
{
MESSAGE("Action_hessenberg Ctor");
// STL vector initialization
typename Interface::stl_matrix tmp;
init_matrix<pseudo_random>(tmp,_size);
init_matrix<null_function>(X_stl,_size);
STL_interface<typename Interface::real_type>::ata_product(tmp,X_stl,_size);
init_matrix<null_function>(C_stl,_size);
init_matrix<null_function>(resu_stl,_size);
// generic matrix and vector initialization
Interface::matrix_from_stl(X_ref,X_stl);
Interface::matrix_from_stl(X,X_stl);
Interface::matrix_from_stl(C,C_stl);
_cost = 0;
for (int j=0; j<_size-2; ++j)
{
int r = std::max(0,_size-j-1);
int b = std::max(0,_size-j-2);
_cost += 6 + 3*b + r*r*8;
}
}
// invalidate copy ctor
Action_tridiagonalization( const Action_tridiagonalization & )
{
INFOS("illegal call to Action_tridiagonalization Copy Ctor");
exit(1);
}
// Dtor
~Action_tridiagonalization( void ){
MESSAGE("Action_tridiagonalization Dtor");
// deallocation
Interface::free_matrix(X_ref,_size);
Interface::free_matrix(X,_size);
Interface::free_matrix(C,_size);
}
// action name
static inline std::string name( void ) { return "tridiagonalization_"+Interface::name(); }
double nb_op_base( void ){
return _cost;
}
inline void initialize( void ){
Interface::copy_matrix(X_ref,X,_size);
}
inline void calculate( void ) {
Interface::tridiagonalization(X,C,_size);
}
void check_result( void ){
// calculation check
Interface::matrix_to_stl(C,resu_stl);
// STL_interface<typename Interface::real_type>::tridiagonalization(X_stl,C_stl,_size);
//
// typename Interface::real_type error=
// STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
//
// if (error>1.e-6){
// INFOS("WRONG CALCULATION...residual=" << error);
// exit(0);
// }
}
private :
typename Interface::stl_matrix X_stl;
typename Interface::stl_matrix C_stl;
typename Interface::stl_matrix resu_stl;
typename Interface::gene_matrix X_ref;
typename Interface::gene_matrix X;
typename Interface::gene_matrix C;
int _size;
double _cost;
};
#endif

7
bench/btl/cmake/FindATLAS.cmake
View File

@ -18,8 +18,11 @@ find_library(ATLAS_LIB atlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_CBLAS libcblas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_CBLAS cblas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_LAPACK liblapack_atlas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_LAPACK lapack_atlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
# find_file(ATLAS_LAPACK liblapack_atlas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
# find_library(ATLAS_LAPACK lapack_atlas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_LAPACK liblapack.so.3 PATHS /usr/lib/atlas $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
# find_library(ATLAS_LAPACK lapack PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_file(ATLAS_F77BLAS libf77blas.so.3 PATHS /usr/lib $ENV{ATLASDIR} ${LIB_INSTALL_DIR})
find_library(ATLAS_F77BLAS f77blas PATHS $ENV{ATLASDIR} ${LIB_INSTALL_DIR})

22
bench/btl/cmake/FindGOTO.cmake Normal file
View File

@ -0,0 +1,22 @@
# include(FindLibraryWithDebug)
if (GOTO_INCLUDES AND GOTO_LIBRARIES)
set(GOTO_FIND_QUIETLY TRUE)
endif (GOTO_INCLUDES AND GOTO_LIBRARIES)
find_path(GOTO_INCLUDES
NAMES
cblas.h
PATHS
$ENV{GOTODIR}/include
${INCLUDE_INSTALL_DIR}
)
find_file(GOTO_LIBRARIES libgotoblas.so PATHS /usr/lib $ENV{GOTODIR} ${LIB_INSTALL_DIR})
find_library(GOTO_LIBRARIES gotoblas PATHS $ENV{GOTODIR} ${LIB_INSTALL_DIR})
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(GOTO DEFAULT_MSG
GOTO_INCLUDES GOTO_LIBRARIES)
mark_as_advanced(GOTO_INCLUDES GOTO_LIBRARIES)

4
bench/btl/data/CMakeLists.txt
View File

@ -1,7 +1,9 @@
ADD_CUSTOM_TARGET(copy_scripts)
SET(script_files go_mean aat.hh ata.hh axpy.hh axpby.hh cholesky.hh trisolve.hh order_lib mk_mean_script.sh mk_new_gnuplot.sh mk_gnuplot_script.sh matrix_matrix.hh matrix_vector.hh atv.hh perlib_plot_settings.txt gnuplot_common_settings.hh )
SET(script_files go_mean aat.hh ata.hh axpy.hh axpby.hh cholesky.hh trisolve.hh hessenberg.hh tridiagonalization.hh
order_lib mk_mean_script.sh mk_new_gnuplot.sh mk_gnuplot_script.sh
matrix_matrix.hh matrix_vector.hh atv.hh perlib_plot_settings.txt gnuplot_common_settings.hh )
FOREACH(script_file ${script_files})
ADD_CUSTOM_COMMAND(

2
bench/btl/data/axpby.hh
View File

@ -1,3 +1,3 @@
set title "{/*1.5 Y = alpha * X + beta * Y}" 0.000000,0.000000
set xlabel "vector size" 0.000000,0.000000
set xrange [1:1000000]
set xrange [5:1000000]

2
bench/btl/data/axpy.hh
View File

@ -1,3 +1,3 @@
set title "{/*1.5 Y += alpha * X}" 0.000000,0.000000
set xlabel "vector size" 0.000000,0.000000
set xrange [1:1000000]
set xrange [5:1000000]

5
bench/btl/data/go_mean
View File

@ -1,5 +1,5 @@
#! /bin/bash
mkdir $1
mkdir -p $1
##cp ../libs/*/*.dat $1
EIGENDIR=`cat eigen_root_dir.txt`
@ -17,6 +17,7 @@ echo '<ul>'\
'</p>' >> $webpagefilename
source mk_mean_script.sh axpy $1 11 2500 100000 250000 > $1/axpy.html
source mk_mean_script.sh axpby $1 11 2500 100000 250000 > $1/axpby.html
source mk_mean_script.sh matrix_vector $1 11 50 300 1000 > $1/matrix_vector.html
source mk_mean_script.sh atv $1 11 50 300 1000 > $1/atv.html
source mk_mean_script.sh matrix_matrix $1 11 100 300 1000 > $1/matrix_matrix.html
@ -24,6 +25,8 @@ source mk_mean_script.sh aat $1 11 100 300 1000 > $1/aat.html
source mk_mean_script.sh ata $1 11 100 300 1000 > $1/ata.html
source mk_mean_script.sh trisolve $1 11 100 300 1000 > $1/trisolve.html
source mk_mean_script.sh cholesky $1 11 100 300 1000 > $1/cholesky.html
source mk_mean_script.sh hessenberg $1 11 100 300 1000 > $1/hessenberg.html
source mk_mean_script.sh tridiagonalization $1 11 100 300 1000 > $1/tridiagonalization.html
## compile the web page ##

3
bench/btl/data/hessenberg.hh Normal file
View File

@ -0,0 +1,3 @@
set title "{/*1.5 Hessenberg decomposition}" 0.000000,0.000000
set xlabel "matrix size" 0.000000,0.000000
set xrange [6:1000]

3
bench/btl/data/mk_mean_script.sh
View File

@ -10,7 +10,8 @@ WORK_DIR=tmp
mkdir $WORK_DIR
DATA_FILE=`find $DIR -name "*.dat" | grep _${WHAT}`
echo
echo ""
echo $"1..."
for FILE in $DATA_FILE
do
##echo hello world

1
bench/btl/data/perlib_plot_settings.txt
View File

@ -7,4 +7,5 @@ ublas ; with lines lw 3 lt 1 lc rgbcolor "#ff0000"
mtl4 ; with lines lw 3 lt 1 lc rgbcolor "#d18847"
blitz ; with lines lw 3 lt 1 lc rgbcolor "#ff00ff"
F77 ; with lines lw 3 lt 3 lc rgbcolor "#e6e64c"
GOTO ; with lines lw 3 lt 3 lc rgbcolor "#e6bd96"
C ; with lines lw 3 lt 3 lc rgbcolor "#e6bd96"

3
bench/btl/data/tridiagonalization.hh Normal file
View File

@ -0,0 +1,3 @@
set title "{/*1.5 Tridiagonalization}" 0.000000,0.000000
set xlabel "matrix size" 0.000000,0.000000
set xrange [6:1000]

10
bench/btl/libs/C_BLAS/CMakeLists.txt
View File

@ -18,3 +18,13 @@ if (MKL_FOUND)
set_target_properties(btl_mkl PROPERTIES COMPILE_FLAGS "-DCBLASNAME=INTEL_MKL -DHAS_LAPACK=1")
endif(BUILD_btl_mkl)
endif (MKL_FOUND)
find_package(GOTO)
if (GOTO_FOUND)
include_directories(${GOTO_INCLUDES} ${PROJECT_SOURCE_DIR}/libs/f77)
btl_add_bench(btl_goto main.cpp)
if(BUILD_btl_goto)
target_link_libraries(btl_goto ${GOTO_LIBRARIES} )
set_target_properties(btl_goto PROPERTIES COMPILE_FLAGS "-DCBLASNAME=GOTO -DPUREBLAS")
endif(BUILD_btl_goto)
endif (GOTO_FOUND)

108
bench/btl/libs/C_BLAS/C_BLAS_interface.hh
View File

@ -25,10 +25,32 @@
extern "C"
{
#include "cblas.h"
#ifdef HAS_LAPACK
void sgemm_(const char *transa, const char *transb, const int *m, const int *n, const int *k,
const float *alpha, const float *a, const int *lda, const float *b, const int *ldb,
const float *beta, float *c, const int *ldc);
void sgemv_(const char *trans, const int *m, const int *n, const float *alpha,
const float *a, const int *lda, const float *x, const int *incx,
const float *beta, float *y, const int *incy);
void sscal_(const int *n, const float *alpha, const float *x, const int *incx);
void saxpy_(const int *n, const float *alpha, const float *x, const int *incx,
float *y, const int *incy);
void strsv_(const char *uplo, const char *trans, const char *diag, const int *n,
const float *a, const int *lda, float *x, const int *incx);
void scopy_(const int *n, const float *x, const int *incx, float *y, const int *incy);
// Cholesky Factorization
// #include "mkl_lapack.h"
void spotrf_(const char* uplo, const int* n, float *a, const int* ld, int* info);
void dpotrf_(const char* uplo, const int* n, double *a, const int* ld, int* info);
void ssytrd_(char *uplo, const int *n, float *a, const int *lda, float *d, float *e, float *tau, float *work, int *lwork, int *info );
void sgehrd_( const int *n, int *ilo, int *ihi, float *a, const int *lda, float *tau, float *work, int *lwork, int *info );
#ifdef HAS_LAPACK
#endif
}
@ -85,9 +107,18 @@ public :
};
static const float fone = 1;
static const float fzero = 0;
static const char notrans = 'N';
static const char trans = 'T';
static const char nonunit = 'N';
static const char lower = 'L';
static const int intone = 1;
template<>
class C_BLAS_interface<float> : public f77_interface_base<float>
{
public :
static inline std::string name( void )
@ -96,59 +127,114 @@ public :
}
static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
#ifdef PUREBLAS
sgemv_(&notrans,&N,&N,&fone,A,&N,B,&intone,&fzero,X,&intone);
#else
cblas_sgemv(CblasColMajor,CblasNoTrans,N,N,1.0,A,N,B,1,0.0,X,1);
#endif
}
static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
#ifdef PUREBLAS
sgemv_(&trans,&N,&N,&fone,A,&N,B,&intone,&fzero,X,&intone);
#else
cblas_sgemv(CblasColMajor,CblasTrans,N,N,1.0,A,N,B,1,0.0,X,1);
#endif
}
static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
#ifdef PUREBLAS
sgemm_(&notrans,&notrans,&N,&N,&N,&fone,A,&N,B,&N,&fzero,X,&N);
#else
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
#endif
}
static inline void transposed_matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
#ifdef PUREBLAS
sgemm_(&notrans,&notrans,&N,&N,&N,&fone,A,&N,B,&N,&fzero,X,&N);
#else
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,1.0,A,N,B,N,0.0,X,N);
#endif
}
static inline void ata_product(gene_matrix & A, gene_matrix & X, int N){
#ifdef PUREBLAS
sgemm_(&trans,&notrans,&N,&N,&N,&fone,A,&N,A,&N,&fzero,X,&N);
#else
cblas_sgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
#endif
}
static inline void aat_product(gene_matrix & A, gene_matrix & X, int N){
#ifdef PUREBLAS
sgemm_(&notrans,&trans,&N,&N,&N,&fone,A,&N,A,&N,&fzero,X,&N);
#else
cblas_sgemm(CblasColMajor,CblasNoTrans,CblasTrans,N,N,N,1.0,A,N,A,N,0.0,X,N);
#endif
}
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N){
#ifdef PUREBLAS
saxpy_(&N,&coef,X,&intone,Y,&intone);
#else
cblas_saxpy(N,coef,X,1,Y,1);
#endif
}
static inline void axpby(float a, const gene_vector & X, float b, gene_vector & Y, int N){
#ifdef PUREBLAS
sscal_(&N,&b,Y,&intone);
saxpy_(&N,&a,X,&intone,Y,&intone);
#else
cblas_sscal(N,b,Y,1);
cblas_saxpy(N,a,X,1,Y,1);
#endif
}
#ifdef HAS_LAPACK
static inline void cholesky(const gene_vector & X, gene_vector & C, int N){
cblas_scopy(N*N, X, 1, C, 1);
char uplo = 'L';
int info = 0;
spotrf_(&uplo, &N, C, &N, &info);
// float tmp[N];
// for (int j = 1; j < N; ++j)
// {
// int endSize = N-j-1;
// if (endSize>0) {
//cblas_sgemv(CblasColMajor, CblasNoTrans, N-j-1, j, ATL_rnone, A+j+1, lda, A+j, lda, ATL_rone, Ac+j+1, 1);
// cblas_sgemv(CblasColMajor, CblasNoTrans,endSize,j, 1.0, &(C[j+1]),N, &(C[j]),N, 0.0, &(C[j+1+N*j]),1);
// }
// }
}
#ifdef HAS_LAPACK
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
cblas_scopy(N*N, X, 1, C, 1);
int info = 0;
int ilo = 1;
int ihi = N;
int bsize = 64;
int worksize = N*bsize;
float* d = new float[N+worksize];
sgehrd_(&N, &ilo, &ihi, C, &N, d, d+N, &worksize, &info);
delete[] d;
}
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){
cblas_scopy(N*N, X, 1, C, 1);
char uplo = 'U';
int info = 0;
int ilo = 1;
int ihi = N;
int bsize = 64;
int worksize = N*bsize;
float* d = new float[3*N+worksize];
ssytrd_(&uplo, &N, C, &N, d, d+N, d+2*N, d+3*N, &worksize, &info);
delete[] d;
}
#endif
static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector & X, int N){
#ifdef PUREBLAS
scopy_(&N, B, &intone, X, &intone);
strsv_(&lower, &notrans, &nonunit, &N, L, &N, X, &intone);
#else
cblas_scopy(N, B, 1, X, 1);
cblas_strsv(CblasColMajor, CblasLower, CblasNoTrans, CblasNonUnit, N, L, N, X, 1);
#endif
}
};

9
bench/btl/libs/C_BLAS/main.cpp
View File

@ -22,8 +22,10 @@
#include "bench.hh"
#include "basic_actions.hh"
#ifdef HAS_LAPACK
#include "action_cholesky.hh"
#ifdef HAS_LAPACK
#include "action_hessenberg.hh"
#endif
BTL_MAIN;
@ -43,8 +45,11 @@ int main()
bench<Action_trisolve<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
// bench<Action_cholesky<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
#ifdef HAS_LAPACK
bench<Action_cholesky<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_hessenberg<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_tridiagonalization<C_BLAS_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
#endif
//bench<Action_lu_solve<C_BLAS_LU_solve_interface<REAL_TYPE> > >(MIN_LU,MAX_LU,NB_POINT);

11
bench/btl/libs/eigen2/eigen2_interface.hh
View File

@ -20,6 +20,7 @@
#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/QR>
#include <vector>
#include "btl.hh"
@ -136,11 +137,17 @@ public :
}
static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){
// C = X;
// Cholesky<gene_matrix>::computeInPlace(C);
C = X.cholesky().matrixL();
}
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
C = HessenbergDecomposition<gene_matrix>(X).packedMatrix();
}
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){
C = Tridiagonalization<gene_matrix>(X).packedMatrix();
}
};
#endif

4
bench/btl/libs/eigen2/main_adv.cpp
View File

@ -20,6 +20,7 @@
#include "bench.hh"
#include "action_trisolve.hh"
#include "action_cholesky.hh"
#include "action_hessenberg.hh"
BTL_MAIN;
@ -28,6 +29,9 @@ int main()
bench<Action_trisolve<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_cholesky<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_hessenberg<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_tridiagonalization<eigen2_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}

2
bench/btl/libs/f77/daxpy.f
View File

@ -1,4 +1,4 @@
SUBROUTINE DAXPY(N,A,X,Y)
SUBROUTINE DAXPYF(N,A,X,Y)
** ***************************************
** CALCULE Y = Y + A*X
** ***************************************

8
bench/btl/libs/f77/f77_interface.hh
View File

@ -36,8 +36,8 @@ extern "C"
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);
void saxpyf_(int * N, float * coef, float * X, float *Y);
void daxpyf_(int * N, double * coef, double * X, double *Y);
}
template<class real>
@ -76,7 +76,7 @@ public :
static inline void axpy(real coef, const gene_vector & X, gene_vector & Y, int N)
{
int one=1;
daxpy_(&N,&coef,X,Y);
daxpyf_(&N,&coef,X,Y);
}
@ -117,7 +117,7 @@ public :
static inline void axpy(float coef, const gene_vector & X, gene_vector & Y, int N)
{
saxpy_(&N,&coef,X,Y);
saxpyf_(&N,&coef,X,Y);
}
};

2
bench/btl/libs/f77/saxpy.f
View File

@ -1,4 +1,4 @@
SUBROUTINE SAXPY(N,A,X,Y)
SUBROUTINE SAXPYF(N,A,X,Y)
** ***************************************
** CALCULE Y = Y + A*X
** ***************************************

9
bench/btl/libs/gmm/gmm_interface.hh
View File

@ -123,6 +123,15 @@ public :
gmm::lower_tri_solve(L, X, false);
}
static inline void hessenberg(const gene_matrix & X, gene_matrix & C, int N){
gmm::copy(X,C);
gmm::Hessenberg_reduction(C,X,false);
}
static inline void tridiagonalization(const gene_matrix & X, gene_matrix & C, int N){
gmm::copy(X,C);
gmm::Householder_tridiagonalization(C,X,false);
}
};

4
bench/btl/libs/gmm/main.cpp
View File

@ -19,6 +19,7 @@
#include "gmm_interface.hh"
#include "bench.hh"
#include "basic_actions.hh"
#include "action_hessenberg.hh"
BTL_MAIN;
@ -38,6 +39,9 @@ int main()
bench<Action_trisolve<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);
bench<Action_hessenberg<gmm_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
bench<Action_tridiagonalization<gmm_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT);
return 0;
}