diff options
author | 2012-10-09 09:18:57 +0200 | |
---|---|---|
committer | 2012-10-09 09:18:57 +0200 | |
commit | 91b7fe13b24657071e06e36a5322344d02a45789 (patch) | |
tree | 75138fec7df9558bffa66044f757c763d27627bd | |
parent | Placed BLAS action under subfolder. Written first LAPACK action. Updated (diff) | |
download | auto-numerical-bench-91b7fe13b24657071e06e36a5322344d02a45789.tar.gz auto-numerical-bench-91b7fe13b24657071e06e36a5322344d02a45789.tar.bz2 auto-numerical-bench-91b7fe13b24657071e06e36a5322344d02a45789.zip |
Fixed GeneralSolve and added LeastSquares.
-rw-r--r-- | btl/NumericInterface/NI_internal/FortranDeclarations.hpp | 8 | ||||
-rw-r--r-- | btl/NumericInterface/NI_internal/FortranInterface.hpp | 13 | ||||
-rw-r--r-- | btl/actions/LAPACK/action_GeneralSolve.hpp (renamed from btl/actions/LAPACK/GeneralSolve.hpp) | 3 | ||||
-rw-r--r-- | btl/actions/LAPACK/action_LeastSquaresSolve.hpp | 84 | ||||
-rw-r--r-- | btl/actions/actionsLAPACK.hpp | 19 |
5 files changed, 122 insertions, 5 deletions
diff --git a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp index 2378cf2..2ae59e5 100644 --- a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp +++ b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp @@ -86,13 +86,13 @@ extern "C" { * LAPACK SOLVERS * ******************/ - void sgesv(const int&, const int&, float*, const int&, int*, float*, const int&, int&); - void dgesv(const int&, const int&, double*, const int&, int*, double*, const int&, int&); + void sgesv_(const int*, const int*, float*, const int*, int*, float*, const int*, int*); + void dgesv_(const int*, const int*, double*, const int*, int*, double*, const int*, int*); + void sgels_(const char*, const int*, const int*, const int*, float*, const int*, float*, const int*, float*, const int*, int*); + void dgels_(const char*, const int*, const int*, const int*, double*, const int*, double*, const int*, double*, const int*, int*); - FORTFUNC(gesv)(&N, &ONE, A, &N, ipiv, b, &N, &info); - } diff --git a/btl/NumericInterface/NI_internal/FortranInterface.hpp b/btl/NumericInterface/NI_internal/FortranInterface.hpp index 2073e80..c6cb3b2 100644 --- a/btl/NumericInterface/NI_internal/FortranInterface.hpp +++ b/btl/NumericInterface/NI_internal/FortranInterface.hpp @@ -28,6 +28,7 @@ # define NI_NAME CAT(CAT(FortranInterface,$),NI_SCALAR) #endif +#include <vector> template<> class NumericInterface<NI_SCALAR> @@ -173,6 +174,18 @@ public: int info; FORTFUNC(gesv)(&N, &ONE, A, &N, ipiv, b, &N, &info); } + + static void LeastSquaresSolve(const int& N, Scalar *A, Scalar *b) + { + int lw, info; + const int mONE = -1; + Scalar lworks; + + FORTFUNC(gels)("N", &N, &N, &ONE, A, &N, b, &N, &lworks, &mONE, &info); + lw = static_cast<int>(lworks); + std::vector<Scalar> work(lw); + FORTFUNC(gels)("N", &N, &N, &ONE, A, &N, b, &N, &work[0], &lw, &info); + } }; const int NumericInterface<NI_SCALAR>::ZERO = 0; diff --git a/btl/actions/LAPACK/GeneralSolve.hpp b/btl/actions/LAPACK/action_GeneralSolve.hpp index 7fb98ba..c6842b1 100644 --- a/btl/actions/LAPACK/GeneralSolve.hpp +++ b/btl/actions/LAPACK/action_GeneralSolve.hpp @@ -68,7 +68,8 @@ public: std::copy(b.begin(), b.end(), b_res.begin()); Interface::MatrixVector(false, _size, _size, -1., &A[0], &x_work[0], 1., &b_res[0]); - return Interface::norm(_size, &b_res[0]); + const Scalar Anorm = Interface::norm(_size*_size, &A[0]); + return Interface::norm(_size, &b_res[0])/Anorm; } private: diff --git a/btl/actions/LAPACK/action_LeastSquaresSolve.hpp b/btl/actions/LAPACK/action_LeastSquaresSolve.hpp new file mode 100644 index 0000000..4bd7da4 --- /dev/null +++ b/btl/actions/LAPACK/action_LeastSquaresSolve.hpp @@ -0,0 +1,84 @@ +//===================================================== +// 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. +// +#ifndef ACTION_LEASTSQUARESSOLVE +#define ACTION_LEASTSQUARESSOLVE + +#include "LinearCongruential.hpp" +#include <vector> +#include <algorithm> + +template<class Interface> +class Action_LeastSquaresSolve { + + typedef typename Interface::Scalar Scalar; + typedef std::vector<Scalar> vector_t; + +private: + // Invalidate copy constructor + Action_LeastSquaresSolve(const Action_LeastSquaresSolve&); + +public: + + // Constructor + Action_LeastSquaresSolve(int size) + : _size(size), lc(10), + A(lc.fillVector<Scalar>(size*size)), b(lc.fillVector<Scalar>(size)), + A_work(size*size), x_work(size), b_res(size) + { + MESSAGE("Action_LeastSquaresSolve Constructor"); + } + + // Action name + static std::string name() + { + return "LeastSquaresSolve_" + Interface::name(); + } + + double fpo() { + return double(_size)*double(_size)*double(_size); + } + + inline void initialize(){ + std::copy(A.begin(), A.end(), A_work.begin()); + std::copy(b.begin(), b.end(), x_work.begin()); + } + + inline void calculate() { + Interface::LeastSquaresSolve(_size, &A_work[0], &x_work[0]); + } + + Scalar getResidual() { + initialize(); + calculate(); + std::copy(b.begin(), b.end(), b_res.begin()); + Interface::MatrixVector(false, _size, _size, -1., &A[0], &x_work[0], + 1., &b_res[0]); + const Scalar Anorm = Interface::norm(_size*_size, &A[0]); + return Interface::norm(_size, &b_res[0])/Anorm; + } + +private: + const int _size; + LinearCongruential<> lc; + + const vector_t A, b; + vector_t A_work, x_work, b_res; + +}; + +#endif // ACTION_LEASTSQUARESSOLVE diff --git a/btl/actions/actionsLAPACK.hpp b/btl/actions/actionsLAPACK.hpp new file mode 100644 index 0000000..3353ec4 --- /dev/null +++ b/btl/actions/actionsLAPACK.hpp @@ -0,0 +1,19 @@ +//===================================================== +// 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 "LAPACK/action_GeneralSolve.hpp" +#include "LAPACK/action_LeastSquaresSolve.hpp" |