diff options
author | 2012-10-15 23:14:41 +0200 | |
---|---|---|
committer | 2012-10-15 23:14:41 +0200 | |
commit | 93de8d4c28ecd6808ebe798b0a8860239d05a023 (patch) | |
tree | f2cbb3261e3c014b66fc9a840419af376d05aa90 | |
parent | Fixed GeneralSolve and added LeastSquares. (diff) | |
download | auto-numerical-bench-93de8d4c28ecd6808ebe798b0a8860239d05a023.tar.gz auto-numerical-bench-93de8d4c28ecd6808ebe798b0a8860239d05a023.tar.bz2 auto-numerical-bench-93de8d4c28ecd6808ebe798b0a8860239d05a023.zip |
Added LU and Cholesky decompositions (errors not working yet)
-rw-r--r-- | btl/NumericInterface/NI_internal/FortranDeclarations.hpp | 13 | ||||
-rw-r--r-- | btl/NumericInterface/NI_internal/FortranInterface.hpp | 20 | ||||
-rw-r--r-- | btl/actions/LAPACK/action_Choleskydecomp.hpp | 81 | ||||
-rw-r--r-- | btl/actions/LAPACK/action_LUdecomp.hpp | 102 | ||||
-rw-r--r-- | btl/actions/actionsLAPACK.hpp | 3 |
5 files changed, 219 insertions, 0 deletions
diff --git a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp index 2ae59e5..30309bf 100644 --- a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp +++ b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp @@ -93,6 +93,19 @@ extern "C" { void dgels_(const char*, const int*, const int*, const int*, double*, const int*, double*, const int*, double*, const int*, int*); + + + /************************* + * LAPACK FACTORIZATIONS * + *************************/ + + void sgetrf_(const int*, const int*, float*, const int*, int*, int*); + void dgetrf_(const int*, const int*, double*, const int*, int*, int*); + + void spotrf_(const char*, const int*, float*, const int*, int*); + void dpotrf_(const char*, const int*, double*, const int*, int*); + + } diff --git a/btl/NumericInterface/NI_internal/FortranInterface.hpp b/btl/NumericInterface/NI_internal/FortranInterface.hpp index c6cb3b2..d9ce272 100644 --- a/btl/NumericInterface/NI_internal/FortranInterface.hpp +++ b/btl/NumericInterface/NI_internal/FortranInterface.hpp @@ -186,6 +186,26 @@ public: std::vector<Scalar> work(lw); FORTFUNC(gels)("N", &N, &N, &ONE, A, &N, b, &N, &work[0], &lw, &info); } + + + + /****************** + * LAPACK SOLVERS * + ******************/ + + static void LUdecomp(const int& N, Scalar* A, int* ipiv) + { + int info; + FORTFUNC(getrf)(&N, &N, A, &N, ipiv, &info); + if (info != 0) std::cerr << "Error: " << info << std::endl; + } + + static void Choleskydecomp(const char& uplo, const int& N, Scalar* A) + { + int info; + FORTFUNC(potrf)(&uplo, &N, A, &N, &info); + if (info != 0) std::cerr << "Error: " << info << std::endl; + } }; const int NumericInterface<NI_SCALAR>::ZERO = 0; diff --git a/btl/actions/LAPACK/action_Choleskydecomp.hpp b/btl/actions/LAPACK/action_Choleskydecomp.hpp new file mode 100644 index 0000000..debb2f1 --- /dev/null +++ b/btl/actions/LAPACK/action_Choleskydecomp.hpp @@ -0,0 +1,81 @@ +//===================================================== +// 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_CHOLESKYDECOMP +#define ACTION_CHOLESKYDECOMP + +#include "LinearCongruential.hpp" +#include <vector> +#include <algorithm> + +template<class Interface> +class Action_Choleskydecomp{ + + typedef typename Interface::Scalar Scalar; + typedef std::vector<Scalar> vector_t; + +private: + // Invalidate copy constructor + Action_Choleskydecomp(const Action_Choleskydecomp&); + +public: + + // Constructor + Action_Choleskydecomp(int size) + : _size(size), lc(10), + A(lc.fillVector<Scalar>(size*size)), A_work(size*size) + { + MESSAGE("Action_Choleskydecomp Constructor"); + + // Make A positive definite + for (int i = 0; i < size; ++i) + A[i+i*size] += 1.2*size; + } + + // Action name + static std::string name() + { + return "Choleskydecomp_" + Interface::name(); + } + + double fpo() { + return double(_size)*double(_size)*double(_size); + } + + inline void initialize(){ + std::copy(A.begin(), A.end(), A_work.begin()); + } + + inline void calculate() { + Interface::Choleskydecomp('U', _size, &A_work[0]); + } + + Scalar getResidual() { + } + +private: + const int _size; + LinearCongruential<> lc; + + const vector_t A; + vector_t A_work, eye_work;; + std::vector<int> ipiv; + +}; + +#endif // ACTION_CHOLESKYDECOMP + diff --git a/btl/actions/LAPACK/action_LUdecomp.hpp b/btl/actions/LAPACK/action_LUdecomp.hpp new file mode 100644 index 0000000..69324de --- /dev/null +++ b/btl/actions/LAPACK/action_LUdecomp.hpp @@ -0,0 +1,102 @@ +//===================================================== +// 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_LUDECOMP +#define ACTION_LUDECOMP + +#include "LinearCongruential.hpp" +#include <vector> +#include <algorithm> + +template<class Interface> +class Action_LUdecomp{ + + typedef typename Interface::Scalar Scalar; + typedef std::vector<Scalar> vector_t; + +private: + // Invalidate copy constructor + Action_LUdecomp(const Action_LUdecomp&); + +public: + + // Constructor + Action_LUdecomp(int size) + : _size(size), lc(10), + A(lc.fillVector<Scalar>(size*size)), A_work(size*size), + eye_work(size*size), ipiv(size) + { + MESSAGE("Action_LUdecomp Constructor"); + + // Make sure A is invertible + for (int i = 0; i < size; ++i) + A[i+i*size] += 1.2*size; + } + + // Action name + static std::string name() + { + return "LUdecomp_" + Interface::name(); + } + + double fpo() { + return double(_size)*double(_size)*double(_size); + } + + inline void initialize(){ + std::copy(A.begin(), A.end(), A_work.begin()); + } + + inline void calculate() { + Interface::LUdecomp(_size, &A_work[0], &ipiv[0]); + } + + Scalar getResidual() { + initialize(); + calculate(); + + // Initialize identity + for (int r = 0; r < _size; ++r) + for (int c = 0; c < _size; ++c) + eye_work[r+_size*c] = (r == c); + + Interface::TriMatrixMatrix('u', _size, _size, &A_work[0], &eye_work[0]); + + // FIXME: hard-coded unitary diagonal + for (int r = 0; r < _size; ++r) + A_work[r+_size*r] = 1.; + + Interface::TriMatrixMatrix('l', _size, _size, &A_work[0], &eye_work[0]); + + Interface::axpy(_size*_size, -1., &A[0], &eye_work[0]); + Scalar n = Interface::norm(_size*_size, &eye_work[0]); + //n /= Interface::norm(_size*_size, &A[0]); + return n; + } + +private: + const int _size; + LinearCongruential<> lc; + + const vector_t A; + vector_t A_work, eye_work;; + std::vector<int> ipiv; + +}; + +#endif // ACTION_LUDECOMP + diff --git a/btl/actions/actionsLAPACK.hpp b/btl/actions/actionsLAPACK.hpp index 3353ec4..11eb379 100644 --- a/btl/actions/actionsLAPACK.hpp +++ b/btl/actions/actionsLAPACK.hpp @@ -17,3 +17,6 @@ // #include "LAPACK/action_GeneralSolve.hpp" #include "LAPACK/action_LeastSquaresSolve.hpp" + +#include "LAPACK/action_LUdecomp.hpp" +#include "LAPACK/action_Choleskydecomp.hpp" |