diff options
| -rw-r--r-- | accuracy/TestAccuracy.hpp | 3 | ||||
| -rw-r--r-- | accuracy/actions/ActionGESV.hpp | 94 | ||||
| -rw-r--r-- | accuracy/libraries/LAPACK/main.cpp | 17 |
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> |