Fixed GeneralSolve and added LeastSquares.

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"