summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--accuracy/TestAccuracy.hpp3
-rw-r--r--accuracy/actions/ActionGESV.hpp94
-rw-r--r--accuracy/libraries/LAPACK/main.cpp17
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>