BTL-based lapack test begin. general_solve (xgesv) and least_squares

(xgels) implemented. Written main.cpp (not tested).

6 files changed, 364 insertions, 11 deletions

diff --git a/btl/actions/action_general_solve.hh b/btl/actions/action_general_solve.hh
new file mode 100644
index 0000000..e18795c
--- /dev/null
+++ b/btl/actions/action_general_solve.hh
@@ -0,0 +1,120 @@
+#ifndef ACTION_GENERAL_SOLVE
+#define ACTION_GENERAL_SOLVE
+
+#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_general_solve {
+
+public:
+
+	// Ctor
+	Action_general_solve( int size ) : _size(size)
+	{
+	    MESSAGE("Action_general_solve Ctor");
+
+	    // STL matrix and vector initialization
+
+	    init_matrix<pseudo_random>(A_stl,_size);
+	    init_vector<pseudo_random>(B_stl,_size);
+	    init_vector<null_function>(X_stl,_size);
+	    init_vector<null_function>(resu_stl,_size);
+
+	    // generic matrix and vector initialization
+
+	    Interface::matrix_from_stl(A_ref,A_stl);
+	    Interface::matrix_from_stl(A,A_stl);
+	    Interface::vector_from_stl(B_ref,B_stl);
+	    Interface::vector_from_stl(B,B_stl);
+	    Interface::vector_from_stl(X_ref,X_stl);
+	    Interface::vector_from_stl(X,X_stl);
+
+
+	}
+
+	// invalidate copy ctor
+
+	Action_general_solve( const  Action_general_solve & )
+	{
+	    INFOS("illegal call to Action_general_solve Copy Ctor");
+	    exit(0);
+	}
+
+	// Dtor
+
+	BTL_DONT_INLINE ~Action_general_solve( void ){
+
+	    MESSAGE("Action_general_solve Dtor");
+
+	    // deallocation
+
+	    Interface::free_matrix(A,_size);
+	    Interface::free_vector(B);
+	    Interface::free_vector(X);
+
+	    Interface::free_matrix(A_ref,_size);
+	    Interface::free_vector(B_ref);
+	    Interface::free_vector(X_ref);
+
+	}
+
+	  // action name
+
+	  static inline std::string name( void )
+	  {
+	    return "general_solve_" + Interface::name();
+	  }
+
+	  double nb_op_base( void ){
+	    return 2.0*_size*_size*_size;
+	  }
+
+	  BTL_DONT_INLINE  void initialize( void ){
+
+	    Interface::copy_matrix(A_ref,A,_size);
+	    Interface::copy_vector(B_ref,B,_size);
+	    Interface::copy_vector(X_ref,X,_size);
+
+	  }
+
+	  BTL_DONT_INLINE void calculate( void ) {
+	      BTL_ASM_COMMENT("#begin general_solve");
+	      Interface::general_solve(A,B,X,_size);
+	      BTL_ASM_COMMENT("end general_solve");
+	  }
+
+	  BTL_DONT_INLINE void check_result() {
+		  // TODO: check result
+	  }
+
+
+private:
+
+	  typename Interface::stl_matrix A_stl;
+	  typename Interface::stl_vector B_stl;
+	  typename Interface::stl_vector X_stl;
+	  typename Interface::stl_vector resu_stl;
+
+	  typename Interface::gene_matrix A_ref;
+	  typename Interface::gene_vector B_ref;
+	  typename Interface::gene_vector X_ref;
+
+	  typename Interface::gene_matrix A;
+	  typename Interface::gene_vector B;
+	  typename Interface::gene_vector X;
+
+
+	  int _size;
+
+
+};
+
+
+#endif// ACTION_GENERAL_SOLVE
diff --git a/btl/actions/action_least_squares.hh b/btl/actions/action_least_squares.hh
new file mode 100644
index 0000000..a37ac66
--- /dev/null
+++ b/btl/actions/action_least_squares.hh
@@ -0,0 +1,120 @@
+#ifndef ACTION_LEAST_SQUARES
+#define ACTION_LEAST_SQUARES
+
+#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_least_squares {
+
+public:
+
+	// Ctor
+	Action_least_squares( int size ) : _size(size)
+	{
+	    MESSAGE("Action_least_squares Ctor");
+
+	    // STL matrix and vector initialization
+
+	    init_matrix<pseudo_random>(A_stl,_size);
+	    init_vector<pseudo_random>(B_stl,_size);
+	    init_vector<null_function>(X_stl,_size);
+	    init_vector<null_function>(resu_stl,_size);
+
+	    // generic matrix and vector initialization
+
+	    Interface::matrix_from_stl(A_ref,A_stl);
+	    Interface::matrix_from_stl(A,A_stl);
+	    Interface::vector_from_stl(B_ref,B_stl);
+	    Interface::vector_from_stl(B,B_stl);
+	    Interface::vector_from_stl(X_ref,X_stl);
+	    Interface::vector_from_stl(X,X_stl);
+
+
+	}
+
+	// invalidate copy ctor
+
+	Action_least_squares( const  Action_least_squares & )
+	{
+	    INFOS("illegal call to Action_least_squares Copy Ctor");
+	    exit(0);
+	}
+
+	// Dtor
+
+	BTL_DONT_INLINE ~Action_least_squares( void ){
+
+	    MESSAGE("Action_least_squares Dtor");
+
+	    // deallocation
+
+	    Interface::free_matrix(A,_size);
+	    Interface::free_vector(B);
+	    Interface::free_vector(X);
+
+	    Interface::free_matrix(A_ref,_size);
+	    Interface::free_vector(B_ref);
+	    Interface::free_vector(X_ref);
+
+	}
+
+	  // action name
+
+	  static inline std::string name( void )
+	  {
+	    return "least_squares_" + Interface::name();
+	  }
+
+	  double nb_op_base( void ){
+	    return 2.0*_size*_size*_size;
+	  }
+
+	  BTL_DONT_INLINE  void initialize( void ){
+
+	    Interface::copy_matrix(A_ref,A,_size);
+	    Interface::copy_vector(B_ref,B,_size);
+	    Interface::copy_vector(X_ref,X,_size);
+
+	  }
+
+	  BTL_DONT_INLINE void calculate( void ) {
+	      BTL_ASM_COMMENT("#begin least_squares");
+	      Interface::least_squares(A,B,X,_size);
+	      BTL_ASM_COMMENT("end least_squares");
+	  }
+
+	  BTL_DONT_INLINE void check_result() {
+		  // TODO: check result
+	  }
+
+
+private:
+
+	  typename Interface::stl_matrix A_stl;
+	  typename Interface::stl_vector B_stl;
+	  typename Interface::stl_vector X_stl;
+	  typename Interface::stl_vector resu_stl;
+
+	  typename Interface::gene_matrix A_ref;
+	  typename Interface::gene_vector B_ref;
+	  typename Interface::gene_vector X_ref;
+
+	  typename Interface::gene_matrix A;
+	  typename Interface::gene_vector B;
+	  typename Interface::gene_vector X;
+
+
+	  int _size;
+
+
+};
+
+
+#endif// ACTION_LEAST_SQUARES
diff --git a/btl/libs/BLAS/cblas_interface_impl.hh b/btl/libs/BLAS/cblas_interface_impl.hh
index e481c57..61caa6b 100644
--- a/btl/libs/BLAS/cblas_interface_impl.hh
+++ b/btl/libs/BLAS/cblas_interface_impl.hh
@@ -15,19 +15,19 @@ public :
   }
 
   static inline void matrix_vector_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
-    BLAS_FUNC(gemv)(CblasRowMajor,CblasNoTrans,N,N,fone,A,N,B,intone,fzero,X,intone);
+    BLAS_FUNC(gemv)(CblasColMajor,CblasNoTrans,N,N,fone,A,N,B,intone,fzero,X,intone);
   }
 
   static inline void symv(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
-    BLAS_FUNC(symv)(CblasRowMajor,CblasLower,N,fone,A,N,B,intone,fzero,X,intone);
+    BLAS_FUNC(symv)(CblasColMajor,CblasLower,N,fone,A,N,B,intone,fzero,X,intone);
   }
 
   static inline void syr2(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
-    BLAS_FUNC(syr2)(CblasRowMajor,CblasLower,N,fone,B,intone,X,intone,A,N);
+    BLAS_FUNC(syr2)(CblasColMajor,CblasLower,N,fone,B,intone,X,intone,A,N);
   }
 
   static inline void ger(gene_matrix & A, gene_vector & X, gene_vector & Y, int N){
-    BLAS_FUNC(ger)(CblasRowMajor,N,N,fone,X,intone,Y,intone,A,N);
+    BLAS_FUNC(ger)(CblasColMajor,N,N,fone,X,intone,Y,intone,A,N);
   }
 
   static inline void rot(gene_vector & A,  gene_vector & B, SCALAR c, SCALAR s, int N){
@@ -35,19 +35,19 @@ public :
   }
 
   static inline void atv_product(gene_matrix & A, gene_vector & B, gene_vector & X, int N){
-    BLAS_FUNC(gemv)(CblasRowMajor,CblasTrans,N,N,fone,A,N,B,intone,fzero,X,intone);
+    BLAS_FUNC(gemv)(CblasColMajor,CblasTrans,N,N,fone,A,N,B,intone,fzero,X,intone);
   }
 
   static inline void matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
-    BLAS_FUNC(gemm)(CblasRowMajor,CblasNoTrans,CblasNoTrans,N,N,N,fone,A,N,B,N,fzero,X,N);
+    BLAS_FUNC(gemm)(CblasColMajor,CblasNoTrans,CblasNoTrans,N,N,N,fone,A,N,B,N,fzero,X,N);
   }
 
   static inline void transposed_matrix_matrix_product(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
-    BLAS_FUNC(gemm)(CblasRowMajor,CblasTrans,CblasNoTrans,N,N,N,fone,A,N,B,N,fzero,X,N);
+    BLAS_FUNC(gemm)(CblasColMajor,CblasTrans,CblasNoTrans,N,N,N,fone,A,N,B,N,fzero,X,N);
   }
 
   static inline void aat_product(gene_matrix & A, gene_matrix & X, int N){
-    BLAS_FUNC(syrk)(CblasRowMajor,CblasLower,CblasNoTrans,N,N,fone,A,N,fzero,X,N);
+    BLAS_FUNC(syrk)(CblasColMajor,CblasLower,CblasNoTrans,N,N,fone,A,N,fzero,X,N);
   }
 
   static inline void axpy(SCALAR coef, const gene_vector & X, gene_vector & Y, int N){
@@ -61,16 +61,16 @@ public :
 
   static inline void trisolve_lower(const gene_matrix & L, const gene_vector& B, gene_vector & X, int N){
     BLAS_FUNC(copy)(N, B, intone, X, intone);
-    BLAS_FUNC(trsv)(CblasRowMajor,CblasLower, CblasNoTrans, CblasNonUnit, N, L, N, X, intone);
+    BLAS_FUNC(trsv)(CblasColMajor,CblasLower, CblasNoTrans, CblasNonUnit, N, L, N, X, intone);
   }
 
   static inline void trisolve_lower_matrix(const gene_matrix & L, const gene_matrix& B, gene_matrix & X, int N){
     BLAS_FUNC(copy)(N, B, intone, X, intone);
-    BLAS_FUNC(trsm)(CblasRowMajor,CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, N, N, fone, L, N, X, N);
+    BLAS_FUNC(trsm)(CblasColMajor,CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, N, N, fone, L, N, X, N);
   }
 
   static inline void trmm(gene_matrix & A, gene_matrix & B, gene_matrix & X, int N){
-    BLAS_FUNC(trmm)(CblasRowMajor,CblasLeft, CblasLower, CblasNoTrans,CblasNonUnit, N,N,fone,A,N,B,N);
+    BLAS_FUNC(trmm)(CblasColMajor,CblasLeft, CblasLower, CblasNoTrans,CblasNonUnit, N,N,fone,A,N,B,N);
   }
 
 };
diff --git a/btl/libs/LAPACK/lapack_interface.hh b/btl/libs/LAPACK/lapack_interface.hh
new file mode 100644
index 0000000..73aab38
--- /dev/null
+++ b/btl/libs/LAPACK/lapack_interface.hh
@@ -0,0 +1,44 @@
+#include <../BLAS/c_interface_base.h>
+#include <complex>
+
+extern "C" {
+#include "../BLAS/blas.h"
+
+//void sgesv_(int*, int*, float *, int*, int*, float *, int*, int*);
+//void dgesv_(int*, int*, double*, int*, int*, double*, int*, int*);
+
+void sgels_(char*, int*, int*, int*, float *, int*, float *, int*, float *, int*, int*);
+void dgels_(char*, int*, int*, int*, double*, int*, double*, int*, double*, int*, int*);
+}
+
+
+#define MAKE_STRING2(S) #S
+#define MAKE_STRING(S) MAKE_STRING2(S)
+
+#define CAT2(A,B) A##B
+#define CAT(A,B) CAT2(A,B)
+
+template <typename real> class lapack_interface;
+
+
+static char notrans = 'N';
+static char trans = 'T';
+static char nonunit = 'N';
+static char lower = 'L';
+static char right = 'R';
+static char left = 'L';
+static int intone = 1;
+static int zeroint = 0;
+
+
+#define SCALAR        float
+#define SCALAR_PREFIX s
+#include "lapack_interface_impl.hh"
+#undef SCALAR
+#undef SCALAR_PREFIX
+
+#define SCALAR        double
+#define SCALAR_PREFIX d
+#include "lapack_interface_impl.hh"
+#undef SCALAR
+#undef SCALAR_PREFIX
diff --git a/btl/libs/LAPACK/lapack_interface_impl.hh b/btl/libs/LAPACK/lapack_interface_impl.hh
new file mode 100644
index 0000000..38e3639
--- /dev/null
+++ b/btl/libs/LAPACK/lapack_interface_impl.hh
@@ -0,0 +1,37 @@
+#define LPF(NAME) CAT(CAT(SCALAR_PREFIX,NAME),_)
+
+template<> class lapack_interface<SCALAR> : public c_interface_base<SCALAR>
+{
+public:
+
+	static inline std::string name()
+	{
+		return MAKE_STRING(LAPACKNAME);
+	}
+
+	static inline void general_solve(gene_matrix& A, gene_vector& b, gene_vector& x, int N)
+	{
+		int *ipiv = new int[N];
+		int info;
+		LPF(copy)(&N, b, &intone, x, &intone);
+		LPF(gesv)(&N, &intone, A, &N, ipiv, x, &N, &info);
+		delete[] ipiv;
+	}
+
+	static inline void least_squares(gene_matrix& A, gene_vector& b, gene_vector& x, int N)
+	{
+		int *ipiv = new int[N];
+		int info;
+		LPF(copy)(&N, b, &intone, x, &intone);
+		SCALAR work1;
+		int MONE = -1;
+		LPF(gels)(&notrans, &N, &N, &intone, A, &N, x, &N, &work1, &MONE, &info);
+		int lwork = (int)work1;
+		SCALAR *work2 = new SCALAR[lwork];
+		LPF(gels)(&notrans, &N, &N, &intone, A, &N, x, &N, work2, &lwork, &info);
+		delete[] work2;
+		delete[] ipiv;
+	}
+
+
+};
diff --git a/btl/libs/LAPACK/main.cpp b/btl/libs/LAPACK/main.cpp
new file mode 100644
index 0000000..caa076c
--- /dev/null
+++ b/btl/libs/LAPACK/main.cpp
@@ -0,0 +1,32 @@
+#include <btl.hh>
+#include <bench.hh>
+#include <action_general_solve.hh>
+#include <action_least_squares.hh>
+#include <lapack_interface.hh>
+
+#include <string>
+
+BTL_MAIN;
+
+int main(int argc, char **argv)
+{
+	bool
+	general_solve=false, least_squares=false
+	;
+
+
+	for (int i = 1; i < argv; ++i) {
+		std::string arg = argc[i];
+		if (arg == "general_solve") general_solve = true;
+		else if (arg == "least_squares") least_squares = true;
+	}
+
+
+	if (general_solve)
+	bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(0,2000,NB_POINT);
+
+	if (least_squares)
+	bench<Action_least_squares<lapack_interface<REAL_TYPE> > >(0,2000,NB_POINT);
+
+	return 0;
+}