BTL: - added tridiagonalization and hessenberg decomposition bench

     - added GOTO library

21 files changed, 421 insertions, 28 deletions

diff --git a/bench/btl/actions/action_hessenberg.hh b/bench/btl/actions/action_hessenberg.hh
new file mode 100644
index 000000000..3cb5034e4
--- /dev/null
+++ b/bench/btl/actions/action_hessenberg.hh
@@ -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
diff --git a/bench/btl/cmake/FindATLAS.cmake b/bench/btl/cmake/FindATLAS.cmake
index 1a7d7e85a..91a27a0bd 100644
--- a/bench/btl/cmake/FindATLAS.cmake
+++ b/bench/btl/cmake/FindATLAS.cmake
@@ -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})
diff --git a/bench/btl/cmake/FindGOTO.cmake b/bench/btl/cmake/FindGOTO.cmake
new file mode 100644
index 000000000..b2b648b14
--- /dev/null
+++ b/bench/btl/cmake/FindGOTO.cmake
@@ -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)
diff --git a/bench/btl/data/CMakeLists.txt b/bench/btl/data/CMakeLists.txt
index 9b7e627d3..345e44fad 100644
--- a/bench/btl/data/CMakeLists.txt
+++ b/bench/btl/data/CMakeLists.txt
@@ -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(
diff --git a/bench/btl/data/axpby.hh b/bench/btl/data/axpby.hh
index 744e3063f..84897d8e2 100644
--- a/bench/btl/data/axpby.hh
+++ b/bench/btl/data/axpby.hh
@@ -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]
diff --git a/bench/btl/data/axpy.hh b/bench/btl/data/axpy.hh
index 47a652964..d3fe8df69 100644
--- a/bench/btl/data/axpy.hh
+++ b/bench/btl/data/axpy.hh
@@ -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]
diff --git a/bench/btl/data/go_mean b/bench/btl/data/go_mean
index 05835a8c4..7ca94e091 100755
--- a/bench/btl/data/go_mean
+++ b/bench/btl/data/go_mean
@@ -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 ##
 
diff --git a/bench/btl/data/hessenberg.hh b/bench/btl/data/hessenberg.hh
new file mode 100644
index 000000000..24d889f68
--- /dev/null
+++ b/bench/btl/data/hessenberg.hh
@@ -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]
diff --git a/bench/btl/data/mk_mean_script.sh b/bench/btl/data/mk_mean_script.sh
index 77337d86f..3d9fa1d58 100644
--- a/bench/btl/data/mk_mean_script.sh
+++ b/bench/btl/data/mk_mean_script.sh
@@ -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
diff --git a/bench/btl/data/perlib_plot_settings.txt b/bench/btl/data/perlib_plot_settings.txt
index 536a0ae7b..497914be1 100644
--- a/bench/btl/data/perlib_plot_settings.txt
+++ b/bench/btl/data/perlib_plot_settings.txt
@@ -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"
diff --git a/bench/btl/data/tridiagonalization.hh b/bench/btl/data/tridiagonalization.hh
new file mode 100644
index 000000000..9124398f0
--- /dev/null
+++ b/bench/btl/data/tridiagonalization.hh
@@ -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]
diff --git a/bench/btl/libs/C_BLAS/CMakeLists.txt b/bench/btl/libs/C_BLAS/CMakeLists.txt
index f72e9caea..d567a00d0 100644
--- a/bench/btl/libs/C_BLAS/CMakeLists.txt
+++ b/bench/btl/libs/C_BLAS/CMakeLists.txt
@@ -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)
diff --git a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
index 49ae90f3d..a3ea73740 100644
--- a/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
+++ b/bench/btl/libs/C_BLAS/C_BLAS_interface.hh
@@ -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
   }
 
 };
diff --git a/bench/btl/libs/C_BLAS/main.cpp b/bench/btl/libs/C_BLAS/main.cpp
index 3e38d6985..81b350325 100644
--- a/bench/btl/libs/C_BLAS/main.cpp
+++ b/bench/btl/libs/C_BLAS/main.cpp
@@ -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);
diff --git a/bench/btl/libs/eigen2/eigen2_interface.hh b/bench/btl/libs/eigen2/eigen2_interface.hh
index 13c7058ed..5a85be54d 100644
--- a/bench/btl/libs/eigen2/eigen2_interface.hh
+++ b/bench/btl/libs/eigen2/eigen2_interface.hh
@@ -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
diff --git a/bench/btl/libs/eigen2/main_adv.cpp b/bench/btl/libs/eigen2/main_adv.cpp
index db38b67f4..243edb98f 100644
--- a/bench/btl/libs/eigen2/main_adv.cpp
+++ b/bench/btl/libs/eigen2/main_adv.cpp
@@ -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;
 }
 
diff --git a/bench/btl/libs/f77/daxpy.f b/bench/btl/libs/f77/daxpy.f
index 94bb846de..29514a222 100644
--- a/bench/btl/libs/f77/daxpy.f
+++ b/bench/btl/libs/f77/daxpy.f
@@ -1,4 +1,4 @@
-      SUBROUTINE DAXPY(N,A,X,Y)
+      SUBROUTINE DAXPYF(N,A,X,Y)
 **  ***************************************
 **    CALCULE  Y = Y + A*X
 **  ***************************************
diff --git a/bench/btl/libs/f77/f77_interface.hh b/bench/btl/libs/f77/f77_interface.hh
index c8c8ee3a1..8267dbd42 100644
--- a/bench/btl/libs/f77/f77_interface.hh
+++ b/bench/btl/libs/f77/f77_interface.hh
@@ -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);
   }
 
 };
diff --git a/bench/btl/libs/f77/saxpy.f b/bench/btl/libs/f77/saxpy.f
index 5eda4c2d2..d0f74fd70 100644
--- a/bench/btl/libs/f77/saxpy.f
+++ b/bench/btl/libs/f77/saxpy.f
@@ -1,4 +1,4 @@
-      SUBROUTINE SAXPY(N,A,X,Y)
+      SUBROUTINE SAXPYF(N,A,X,Y)
 **  ***************************************
 **    CALCULE  Y = Y + A*X
 **  ***************************************
diff --git a/bench/btl/libs/gmm/gmm_interface.hh b/bench/btl/libs/gmm/gmm_interface.hh
index 7ef41c0cf..93f827f9b 100644
--- a/bench/btl/libs/gmm/gmm_interface.hh
+++ b/bench/btl/libs/gmm/gmm_interface.hh
@@ -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);
+  }
 
 };
 
diff --git a/bench/btl/libs/gmm/main.cpp b/bench/btl/libs/gmm/main.cpp
index 26fdf76bd..b7c36fc9a 100644
--- a/bench/btl/libs/gmm/main.cpp
+++ b/bench/btl/libs/gmm/main.cpp
@@ -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;
 }