New BTL module lapacke.

5 files changed, 160 insertions, 66 deletions

diff --git a/btl/libs/LAPACK/lapack_.hh b/btl/libs/LAPACK/lapack_.hh
index e7c8c3f..bc56dda 100644
--- a/btl/libs/LAPACK/lapack_.hh
+++ b/btl/libs/LAPACK/lapack_.hh
@@ -32,6 +32,8 @@ void LAPACKFUNC(gesvd)(const char*, const char*, const int*, const int*, SCALAR*
 void LAPACKFUNC(syev)(const char*, const char*, const int*, SCALAR*, const int*, SCALAR*, SCALAR*, const int*, int*);
 void LAPACKFUNC(stev)(const char*, const int*, SCALAR*, SCALAR*, SCALAR*, const int*, SCALAR*, int*);
 
+void LAPACKFUNC(gels)(char*, int*, int*, int*, SCALAR*, int*, SCALAR*, int*, SCALAR*, int*, int*);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/btl/libs/LAPACK/lapack_interface.hh b/btl/libs/LAPACK/lapack_interface.hh
index a0ec52b..87b3cbb 100644
--- a/btl/libs/LAPACK/lapack_interface.hh
+++ b/btl/libs/LAPACK/lapack_interface.hh
@@ -20,26 +20,7 @@
 
 #include <../BLAS/c_interface_base.h>
 #include <complex>
-#include "lapack.hh"
-
-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*);
-
-//void sgetrf_(int*, int*, float *, int*, int*, int*);
-//void dgetrf_(int*, int*, double*, int*, int*, int*);
-
-//void spotrf_(char*, int*, float *, int*, int*);
-//void dpotrf_(char*, int*, double*, int*, int*);
-
-//void ssyev_(char*, char*, int*, float *, int*, float *, float *, int*, int*);
-//void dsyev_(char*, char*, int*, double*, int*, double*, double*, int*, int*);
-}
+//#include "lapack.hh"
 
 
 #define MAKE_STRING2(S) #S
@@ -65,13 +46,21 @@ static int zeroint = 0;
 
 #define SCALAR        float
 #define SCALAR_PREFIX s
-#include "lapack_interface_impl.hh"
+#ifdef LAPACKE_INTERFACE
+#  include "lapacke_interface_impl.hh"
+#else
+#  include "lapack_interface_impl.hh"
+#endif
 #undef SCALAR
 #undef SCALAR_PREFIX
 
 #define SCALAR        double
 #define SCALAR_PREFIX d
-#include "lapack_interface_impl.hh"
+#ifdef LAPACKE_INTERFACE
+#  include "lapacke_interface_impl.hh"
+#else
+#  include "lapack_interface_impl.hh"
+#endif
 #undef SCALAR
 #undef SCALAR_PREFIX
 
diff --git a/btl/libs/LAPACK/lapack_interface_impl.hh b/btl/libs/LAPACK/lapack_interface_impl.hh
index a95bc5b..523e53a 100644
--- a/btl/libs/LAPACK/lapack_interface_impl.hh
+++ b/btl/libs/LAPACK/lapack_interface_impl.hh
@@ -17,6 +17,9 @@
 //
 #define LAPACKFUNC(NAME) CAT(CAT(SCALAR_PREFIX,NAME),_)
 
+#include "lapack_.hh"
+#include "blas.h"
+
 template<> class lapack_interface<SCALAR> : public c_interface_base<SCALAR>
 {
 public:
@@ -114,7 +117,6 @@ public:
 	static inline void stev(const gene_vector& D, const gene_vector& E, gene_vector& W, gene_matrix& V, const int& N)
 	{
 		int N0 = N;
-		int N1 = N-1;
 		LAPACKFUNC(copy)(&N0, D, &intone, W, &intone);
 		stl_vector E_(E, E+N-1), work(max(1, 2*N-2));
 
diff --git a/btl/libs/LAPACK/lapacke_interface_impl.hh b/btl/libs/LAPACK/lapacke_interface_impl.hh
new file mode 100644
index 0000000..eef3638
--- /dev/null
+++ b/btl/libs/LAPACK/lapacke_interface_impl.hh
@@ -0,0 +1,101 @@
+//=====================================================
+// Copyright (C) 2012 Andrea Arteaga <andyspiros@gmail.com>
+//=====================================================
+//
+// 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.
+//
+#define LAPACKEFUNC(NAME) CAT(LAPACKE_,CAT(SCALAR_PREFIX,NAME))
+#define BLASFUNC(NAME) CAT(CAT(SCALAR_PREFIX,NAME),_)
+
+#include "lapacke.h"
+
+
+// Define BLAS functions used by LAPACKE interface
+extern "C" {
+    int BLASFUNC(copy) (int *, SCALAR*, int *, SCALAR*, int *);
+}
+
+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)
+	{
+		std::vector<int> ipiv(N);
+		BLASFUNC(copy)(&N, b, &intone, x, &intone);
+		LAPACKEFUNC(gesv)(LAPACK_COL_MAJOR, N, 1, A, N, &ipiv[0], x, N);
+	}
+
+	static inline void least_squares(gene_matrix& A, gene_vector& b, gene_vector& x, int N)
+	{
+		BLASFUNC(copy)(&N, b, &intone, x, &intone);
+		LAPACKEFUNC(gels)(LAPACK_COL_MAJOR, 'N', N, N, 1, A, N, x, N);
+	}
+
+	static inline void lu_decomp(const gene_matrix& X, gene_matrix& C, int N)
+	{
+		int N2 = N*N;
+		std::vector<int> ipiv(N);
+		BLASFUNC(copy)(&N2, X, &intone, C, &intone);
+		LAPACKEFUNC(getrf)(LAPACK_COL_MAJOR, N, N, C, N, &ipiv[0]);
+	}
+
+	static inline void cholesky(const gene_matrix& X, gene_matrix& C, int N)
+	{
+		int N2 = N*N;
+		BLASFUNC(copy)(&N2, X, &intone, C, &intone);
+		LAPACKEFUNC(potrf)(LAPACK_COL_MAJOR, 'L', N, C, N);
+	}
+	
+	static inline void qr_decomp(const gene_matrix& X, gene_matrix& QR, gene_vector& tau, const int& N)
+	{
+		int N2 = N*N;
+		BLASFUNC(copy)(&N2, X, &intone, QR, &intone);
+		LAPACKEFUNC(geqrf)(LAPACK_COL_MAJOR, N, N, QR, N, tau);
+	}
+
+	static inline void svd_decomp(const gene_matrix& X, gene_matrix& U, gene_vector& S, gene_matrix& VT, const int& N)
+	{
+		int N2 = N*N;
+		stl_vector Xcopy(N2), superb(N-1);
+		BLASFUNC(copy)(&N2, X, &intone, &Xcopy[0], &intone);
+		LAPACKEFUNC(gesvd)(LAPACK_COL_MAJOR, 'A', 'A', N, N, &Xcopy[0], N, S, U, N, VT, N, &superb[0]);
+	}
+
+	static inline void syev(const gene_matrix& X, gene_matrix& V, gene_vector& W, const int& N)
+	{
+		int N2 = N*N;
+		BLASFUNC(copy)(&N2, X, &intone, V, &intone);
+		LAPACKEFUNC(syev)(LAPACK_COL_MAJOR, 'V', 'U', N, V, N, W);
+	}
+
+	/* Size of D, W: N; size of E: N-1, size of V: NxN */
+	static inline void stev(const gene_vector& D, const gene_vector& E, gene_vector& W, gene_matrix& V, int N)
+	{
+		stl_vector E_(E, E+N-1);
+		BLASFUNC(copy)(&N, D, &intone, W, &intone);
+		LAPACKEFUNC(stev)(LAPACK_COL_MAJOR, 'V', N, W, &E_[0], V, N);
+	}
+
+	static inline void symm_ev(const gene_matrix& X, gene_vector& W, int N)
+	{
+		LAPACKEFUNC(syev)(LAPACK_COL_MAJOR, 'N', 'L', N, X, N, W);
+	}
+
+};
diff --git a/btl/libs/LAPACK/main.cpp b/btl/libs/LAPACK/main.cpp
index e70c161..78f314a 100644
--- a/btl/libs/LAPACK/main.cpp
+++ b/btl/libs/LAPACK/main.cpp
@@ -35,61 +35,61 @@ BTL_MAIN;
 
 int main(int argc, char **argv)
 {
-	bool
-	general_solve=false, least_squares=false, lu_decomp=false, cholesky=false, qr_decomp=false, svd_decomp=false,
-	syev=false, stev=false, symm_ev=false
-	;
-	int N = 100;
+    bool
+    general_solve=false, least_squares=false, lu_decomp=false, cholesky=false,
+        qr_decomp=false, svd_decomp=false, syev=false, stev=false,
+        symm_ev=false;
+    int N = 100;
 
 
-	for (int i = 1; i < argc; ++i) {
-		std::string arg = argv[i];
-		if (arg == "general_solve") general_solve = true;
-		else if (arg == "least_squares") least_squares = true;
-		else if (arg == "lu_decomp") lu_decomp = true;
-		else if (arg == "cholesky") cholesky = true;
-		else if (arg == "qr_decomp") qr_decomp = true;
-		else if (arg == "svd_decomp") svd_decomp = true;
-		else if (arg == "syev") syev = true;
-		else if (arg == "stev") stev = true;
-		else if (arg == "symm_ev") symm_ev = true;
+    for (int i = 1; i < argc; ++i) {
+        std::string arg = argv[i];
+        if (arg == "general_solve") general_solve = true;
+        else if (arg == "least_squares") least_squares = true;
+        else if (arg == "lu_decomp") lu_decomp = true;
+        else if (arg == "cholesky") cholesky = true;
+        else if (arg == "qr_decomp") qr_decomp = true;
+        else if (arg == "svd_decomp") svd_decomp = true;
+        else if (arg == "syev") syev = true;
+        else if (arg == "stev") stev = true;
+        else if (arg == "symm_ev") symm_ev = true;
 
-		// Check switch -N
-		else if (arg[0] == '-' && arg[1] == 'N') {
-		    if (arg[2] != '\0')
-		        N = atoi(arg.c_str()+2);
-		    else
-		        N = atoi(argv[++i]);
-		}
-	}
+        // Check switch -N
+        else if (arg[0] == '-' && arg[1] == 'N') {
+            if (arg[2] != '\0')
+                N = atoi(arg.c_str()+2);
+            else
+                N = atoi(argv[++i]);
+        }
+    }
 
 
-	if (general_solve)
-	bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (general_solve)
+    bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	if (least_squares)
-	bench<Action_least_squares<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (least_squares)
+    bench<Action_least_squares<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	if (lu_decomp)
-	bench<Action_lu_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (lu_decomp)
+    bench<Action_lu_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	if (cholesky)
-	bench<Action_cholesky<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (cholesky)
+    bench<Action_cholesky<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	if (qr_decomp)
-	bench<Action_qr_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (qr_decomp)
+    bench<Action_qr_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	if (svd_decomp)
-	bench<Action_svd_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,750,N);
+    if (svd_decomp)
+    bench<Action_svd_decomp<lapack_interface<REAL_TYPE> > >(MIN_MM,750,N);
 
-	if (syev)
-	bench<Action_syev<lapack_interface<REAL_TYPE> > >(MIN_MM,750,N);
+    if (syev)
+    bench<Action_syev<lapack_interface<REAL_TYPE> > >(MIN_MM,750,N);
 
-	if (stev)
-	bench<Action_stev<lapack_interface<REAL_TYPE> > >(MIN_MM,1000,N);
+    if (stev)
+    bench<Action_stev<lapack_interface<REAL_TYPE> > >(MIN_MM,1000,N);
 
-	if (symm_ev)
-	bench<Action_symm_ev<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
+    if (symm_ev)
+    bench<Action_symm_ev<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,N);
 
-	return 0;
+    return 0;
 }