Added eigensolver tests (syev, stev) to lapack_accuracy

5 files changed, 103 insertions, 8 deletions

diff --git a/accuracy/lapack/lapack_STEV.hh b/accuracy/lapack/lapack_STEV.hh
new file mode 100644
index 0000000..0c9dea3
--- /dev/null
+++ b/accuracy/lapack/lapack_STEV.hh
@@ -0,0 +1,35 @@
+#ifndef LAPACK_STEV_HH
+#define LAPACK_STEV_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_STEV(const int& N, const unsigned& seed = 0)
+{
+  LinearCongruential lc(seed);
+  vector<double> A(N*N), D(N), E(N-1), Z(N*N), DD(N*N), work(2*N-2), tmp(N*N);
+
+  /* Fill D, E and A */
+  for (int i = 0; i < N-1; ++i) {
+      D[i] = lc.get_01(); A[i+i*N] = D[i];
+      E[i] = lc.get_01(); A[(i+1)+i*N] = E[i]; A[i+(i+1)*N] = E[i];
+  }
+  D[N-1] = lc.get_01(); A[N*N-1] = D[N-1];
+
+  /* Perform computation */
+  int info;
+  dstev_("V", &N, &D[0], &E[0], &Z[0], &N, &work[0], &info);
+
+  /* Construct DD */
+  for (int i = 0; i < N; ++i)
+    DD[i+i*N] = D[i];
+
+  /* A = V*D*V' */
+  const double dONE = 1., dZERO = 0.;
+  dgemm_("N", "N", &N, &N, &N, &dONE,   &Z[0], &N, &DD[0], &N, &dZERO, &tmp[0], &N);
+  dgemm_("N", "T", &N, &N, &N, &dONE, &tmp[0], &N,  &Z[0], &N, &dZERO,  &DD[0], &N);
+
+  return diff(A, DD);
+}
+
+#endif /* LAPACK_STEV_HH */
diff --git a/accuracy/lapack/lapack_SYEV.hh b/accuracy/lapack/lapack_SYEV.hh
new file mode 100644
index 0000000..a0acad0
--- /dev/null
+++ b/accuracy/lapack/lapack_SYEV.hh
@@ -0,0 +1,43 @@
+#ifndef LAPACK_SYEV_HH
+#define LAPACK_SYEV_HH
+
+#include "LinearCongruential.hh"
+#include "diff.hh"
+
+double test_SYEV(const int& N, const unsigned& seed = 0)
+{
+  LinearCongruential lc(seed);
+  vector<double> A(N*N), Acopy, W(N), D(N*N), work(1), tmp(N*N);
+
+  /* Fill A (symmetric) and Acopy */
+  for (int r = 0; r < N; ++r) {
+      A[r+r*N] = lc.get_01();
+      for (int c = r+1; c < N; ++c) {
+          A[r+c*N] = lc.get_01();
+          A[c+r*N] = A[r+c*N];
+      }
+  }
+  Acopy = A;
+
+  /* Retrieve lwork */
+  int lwork = -1, info;
+  dsyev_("V", "U", &N, &A[0], &N, &W[0], &work[0], &lwork, &info);
+  lwork = work[0];
+  work.resize(lwork);
+
+  /* Perform computation */
+  dsyev_("V", "U", &N, &A[0], &N, &W[0], &work[0], &lwork, &info);
+
+  /* Construct D */
+  for (int i = 0; i < N; ++i)
+    D[i+i*N] = W[i];
+
+  /* A = V*D*V' */
+  const double dONE = 1., dZERO = 0.;
+  dgemm_("N", "N", &N, &N, &N, &dONE,   &A[0], &N, &D[0], &N, &dZERO, &tmp[0], &N);
+  dgemm_("N", "T", &N, &N, &N, &dONE, &tmp[0], &N, &A[0], &N, &dZERO,   &D[0], &N);
+
+  return diff(Acopy, D);
+}
+
+#endif /* LAPACK_SYEV_HH */
diff --git a/accuracy/lapack/main_lapack.cpp b/accuracy/lapack/main_lapack.cpp
index 2426cc6..18b916c 100644
--- a/accuracy/lapack/main_lapack.cpp
+++ b/accuracy/lapack/main_lapack.cpp
@@ -26,6 +26,8 @@ extern "C" {
 #include "lapack_cholesky.hh"
 #include "lapack_SVD.hh"
 #include "lapack_QR.hh"
+#include "lapack_SYEV.hh"
+#include "lapack_STEV.hh"
 
 template<typename exec_t>
 void test(exec_t exec, const std::string& testname, const int& max = 3000, const int& N = 100)
@@ -59,7 +61,8 @@ void test(exec_t exec, const std::string& testname, const int& max = 3000, const
 int main(int argc, char **argv)
 {
   bool
-  lu_decomp=false, cholesky = false, svd_decomp=false, qr_decomp=false
+  lu_decomp=false, cholesky = false, svd_decomp=false, qr_decomp=false,
+  syev=false, stev=false
   ;
 
   for (int i = 1; i < argc; ++i) {
@@ -68,6 +71,8 @@ int main(int argc, char **argv)
           else if (arg == "cholesky") cholesky = true;
           else if (arg == "svd_decomp") svd_decomp = true;
           else if (arg == "qr_decomp") qr_decomp = true;
+          else if (arg == "syev") syev = true;
+          else if (arg == "stev") stev = true;
   }
 
   if (lu_decomp) {
@@ -79,11 +84,19 @@ int main(int argc, char **argv)
   }
 
   if (svd_decomp) {
-      test(test_SVD, "svd_decomp", 3000, 100);
+      test(test_SVD, "svd_decomp", 1000, 90);
   }
 
   if (qr_decomp) {
-      test(test_QR, "qr_decomp", 3000, 100);
+      test(test_QR, "qr_decomp", 500, 70);
+  }
+
+  if (syev) {
+      test(test_SYEV, "syev", 1000, 90);
+  }
+
+  if (stev) {
+      test(test_STEV, "stev", 1000, 90);
   }
 
   return 0;
diff --git a/benchconfig.py b/benchconfig.py
index 6b96b33..c684685 100644
--- a/benchconfig.py
+++ b/benchconfig.py
@@ -17,7 +17,10 @@ if needsinitialization:
     # Script directories
     curdir = os.path.abspath('.')
     scriptdir = os.path.dirname(os.path.realpath(__file__))
-    btldir = '/usr/include/btl'
+    if os.environ.has_key('BTLDIR'):
+        btldir = os.environ['BTLDIR']
+    else:
+        btldir = '/usr/include/btl'
     
     # Library directory (lib32 vs. lib64)
     libdir = sp.Popen \
diff --git a/lapack_accuracy.py b/lapack_accuracy.py
index 1f062fa..ef0c504 100644
--- a/lapack_accuracy.py
+++ b/lapack_accuracy.py
@@ -12,7 +12,8 @@ class Module(basemodule.BaseModule):
     
     def _initialize(self):
         self.libname = 'lapack'
-        self.avail=['lu_decomp', 'cholesky', 'svd_decomp', 'qr_decomp']
+        self.avail=['lu_decomp', 'cholesky', 'svd_decomp', 'qr_decomp', \
+                    'syev', 'stev']
     
     def _parse_args(self, args):     
         # Parse arguments
@@ -125,7 +126,7 @@ class LAPACK_accuracyTest(basemodule.BaseTest):
         logfile.write(' '.join([n+'='+v for n,v in self.runenv.items()]) + ' ')
         logfile.write(' '.join(args) + '\n')
         logfile.write(80*'-' + '\n')
-        proc = sp.Popen(args, bufsize=1, stdout=sp.PIPE, stderr=sp.PIPE, 
+        proc = sp.Popen(args, bufsize=1, stdout=sp.PIPE, stderr=logfile, 
           env=self.runenv, cwd=self.testdir)
         
         # Interpret output
@@ -143,7 +144,7 @@ class LAPACK_accuracyTest(basemodule.BaseTest):
                 Print.down()
             else:
                 Print(line[1:-1])
-        Print.up()          
-        logfile.close()      
+        Print.up()
+        logfile.close()
         proc.wait()
         return proc.returncode