Accuracy: ready for single and double precision.accuracy

3 files changed, 92 insertions, 22 deletions

diff --git a/accuracy/TestAccuracy.hpp b/accuracy/TestAccuracy.hpp
index a79b17a..300d4fb 100644
--- a/accuracy/TestAccuracy.hpp
+++ b/accuracy/TestAccuracy.hpp
@@ -78,8 +78,7 @@ void testAccuracy(
         t.start();
         do {
             Action<value_t> action(size, N);
-            action.compute();
-            e = action.getResidual();
+            e = action.compute();
             emean += e;
             e2mean += e*e;
         } while(++N < 100 && t.elapsed() < 1. || N < 4);
diff --git a/accuracy/actions/ActionGESV.hpp b/accuracy/actions/ActionGESV.hpp
index cd73164..aa9326c 100644
--- a/accuracy/actions/ActionGESV.hpp
+++ b/accuracy/actions/ActionGESV.hpp
@@ -4,9 +4,22 @@
 #include <vector>
 #include <iostream>
 #include <string>
-#include "LinearCongruential.hpp"
+#include "utilities/LinearCongruential.hpp"
+
+/////////////////////////////////////
+// Declaration of fortran routines //
+/////////////////////////////////////
 
 extern "C" {
+    void sgesv_(const int*, const int*, float*, const int*, float*,
+                float*, const int*, int*);
+
+    void sgemv_(const char*, const int*, const int*, const float*,
+                const float*, const int*, const float*, const int*,
+                const float*, float*, const int*);
+
+    float snrm2_(const int*, const float*, const int*);
+
     void dgesv_(const int*, const int*, double*, const int*, double*,
                 double*, const int*, int*);
 
@@ -17,6 +30,11 @@ extern "C" {
     double dnrm2_(const int*, const double*, const int*);
 }
 
+
+/*
+ * Action class
+ */
+
 template<typename value_t = double>
 class ActionGESV {
     typedef std::vector<value_t> storage_t;
@@ -24,26 +42,11 @@ class ActionGESV {
 
 public:
     ActionGESV(const int& size, const unsigned& seed=0)
-        : size(size), rg(seed), A(rg.fillVector<>(size*size, 0.)), Acopy(A),
-          x(rg.fillVector(size, 0.)), b(x), bcopy(b)
+        : size(size), rg(seed), A(rg.fillVector<value_t>(size*size, 0.)),
+          Acopy(A), x(rg.fillVector<float>(size, 0.)), b(x)
     { }
 
-    void compute() {
-        const int ONE = 1;
-        int info;
-        std::vector<value_t> ipiv(size);
-        dgesv_(&size, &ONE, &A[0], &size, &ipiv[0], &x[0], &size, &info);
-        if (info != 0)
-            std::cerr << "Info: " << info << "\n";
-    }
-
-    double getResidual() {
-        const double alpha = -1., beta = 1.;
-        const int ONE = 1;
-        dgemv_("N", &size, &size, &alpha, &Acopy[0], &size, &x[0], &ONE, &beta,
-                &b[0], &ONE);
-        return dnrm2_(&size, &b[0], &ONE)/dnrm2_(&size, &bcopy[0], &ONE);
-    }
+    inline value_t compute();
 
     static std::string fileName() {
         return "accuracy_general_solve.dat";
@@ -52,7 +55,58 @@ public:
 private:
     const int size;
     rangen_t rg;
-    storage_t A, Acopy, x, b, bcopy;
+    storage_t A, Acopy, x, b;
 };
 
+
+
+//////////////////////////////////////
+// compute() method implementations //
+//////////////////////////////////////
+
+template<>
+float ActionGESV<float>::compute()
+{
+    const int ONE = 1;
+    const float alpha = -1., beta = 1.;
+    int info;
+    std::vector<float> ipiv(size);
+
+    // Get input (b) norm
+    const float bnorm = snrm2_(&size, &b[0], &ONE);
+
+    // Perform computation
+    sgesv_(&size, &ONE, &A[0], &size, &ipiv[0], &x[0], &size, &info);
+    if (info != 0)
+        std::cerr << "Info: " << info << "\n";
+
+    // Compute residual
+    sgemv_("N", &size, &size, &alpha, &Acopy[0], &size, &x[0], &ONE, &beta,
+            &b[0], &ONE);
+    return snrm2_(&size, &b[0], &ONE)/bnorm;
+}
+
+
+template<>
+double ActionGESV<double>::compute()
+{
+    const int ONE = 1;
+    const double alpha = -1., beta = 1.;
+    int info;
+    std::vector<double> ipiv(size);
+
+    // Get input (b) norm
+    const double bnorm = dnrm2_(&size, &b[0], &ONE);
+
+    // Perform computation
+    dgesv_(&size, &ONE, &A[0], &size, &ipiv[0], &x[0], &size, &info);
+    if (info != 0)
+        std::cerr << "Info: " << info << "\n";
+
+    // Compute residual
+    dgemv_("N", &size, &size, &alpha, &Acopy[0], &size, &x[0], &ONE, &beta,
+            &b[0], &ONE);
+    return dnrm2_(&size, &b[0], &ONE)/bnorm;
+}
+
 #endif //  ACTIONGESV_HPP
diff --git a/accuracy/libraries/LAPACK/main.cpp b/accuracy/libraries/LAPACK/main.cpp
index ae55cd9..481d5a1 100644
--- a/accuracy/libraries/LAPACK/main.cpp
+++ b/accuracy/libraries/LAPACK/main.cpp
@@ -1,3 +1,20 @@
+//=====================================================
+// 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.
+//
 #include <iostream>
 #include <string>