Added LU and Cholesky decompositions (errors not working yet)

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"