[gentoo-commits] proj/auto-numerical-bench:newinterfaces commit in: btl/actions/, btl/NumericInterface/NI_internal/, btl/actions/LAPACK/

["Andrea Arteaga" <andyspiros@gmail.com>]

commit:     1ee19ae464e67b8aa63ce93a9d7cec49ce8d9c52
Author:     Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
AuthorDate: Mon Oct 15 21:31:07 2012 +0000
Commit:     Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Mon Oct 15 21:31:07 2012 +0000
URL:        http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=1ee19ae4

Added QRdecomp. Solved issue with constant A in LU and Cholesky. Solved
getResidual in LU.

---
 .../NI_internal/FortranDeclarations.hpp            |    3 ++
 .../NI_internal/FortranInterface.hpp               |   14 ++++++++
 btl/actions/LAPACK/action_Choleskydecomp.hpp       |    2 +-
 btl/actions/LAPACK/action_LUdecomp.hpp             |    4 +-
 ...tion_Choleskydecomp.hpp => action_QRdecomp.hpp} |   32 ++++++++++----------
 btl/actions/actionsLAPACK.hpp                      |    1 +
 6 files changed, 37 insertions(+), 19 deletions(-)

diff --git a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp
index 30309bf..60f6afc 100644
--- a/btl/NumericInterface/NI_internal/FortranDeclarations.hpp
+++ b/btl/NumericInterface/NI_internal/FortranDeclarations.hpp
@@ -105,6 +105,9 @@ extern "C" {
     void spotrf_(const char*, const int*, float*, const int*, int*);
     void dpotrf_(const char*, const int*, double*, const int*, int*);
 
+    void sgeqp3_(const int*, const int*, float*, const int*, int*, float*, float*, const int*, int*);
+    void dgeqp3_(const int*, const int*, double*, const int*, int*, double*, double*, const int*, int*);
+
 
 }
 

diff --git a/btl/NumericInterface/NI_internal/FortranInterface.hpp b/btl/NumericInterface/NI_internal/FortranInterface.hpp
index d9ce272..57e0d1d 100644
--- a/btl/NumericInterface/NI_internal/FortranInterface.hpp
+++ b/btl/NumericInterface/NI_internal/FortranInterface.hpp
@@ -206,6 +206,20 @@ public:
         FORTFUNC(potrf)(&uplo, &N, A, &N, &info);
         if (info != 0) std::cerr << "Error: " << info << std::endl;
     }
+
+    static void QRdecomp(const int& N, Scalar* A, int* jpiv, Scalar* tau)
+    {
+        int lw, info;
+        const int mONE = -1;
+        Scalar lworks;
+
+        FORTFUNC(geqp3)(&N, &N, A, &N, jpiv, tau, &lworks, &mONE, &info);
+        lw = static_cast<int>(lworks);
+        std::vector<Scalar> work(lw);
+        FORTFUNC(geqp3)(&N, &N, A, &N, jpiv, tau, &work[0], &lw, &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
index debb2f1..b60bac1 100644
--- a/btl/actions/LAPACK/action_Choleskydecomp.hpp
+++ b/btl/actions/LAPACK/action_Choleskydecomp.hpp
@@ -71,7 +71,7 @@ private:
     const int _size;
     LinearCongruential<> lc;
 
-    const vector_t A;
+    vector_t A;
     vector_t A_work, eye_work;;
     std::vector<int> ipiv;
 

diff --git a/btl/actions/LAPACK/action_LUdecomp.hpp b/btl/actions/LAPACK/action_LUdecomp.hpp
index 69324de..09a161e 100644
--- a/btl/actions/LAPACK/action_LUdecomp.hpp
+++ b/btl/actions/LAPACK/action_LUdecomp.hpp
@@ -84,7 +84,7 @@ public:
 
         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]);
+        n /= Interface::norm(_size*_size, &A[0]);
         return n;
     }
 
@@ -92,7 +92,7 @@ private:
     const int _size;
     LinearCongruential<> lc;
 
-    const vector_t A;
+    vector_t A;
     vector_t A_work, eye_work;;
     std::vector<int> ipiv;
 

diff --git a/btl/actions/LAPACK/action_Choleskydecomp.hpp b/btl/actions/LAPACK/action_QRdecomp.hpp
similarity index 70%
copy from btl/actions/LAPACK/action_Choleskydecomp.hpp
copy to btl/actions/LAPACK/action_QRdecomp.hpp
index debb2f1..a079407 100644
--- a/btl/actions/LAPACK/action_Choleskydecomp.hpp
+++ b/btl/actions/LAPACK/action_QRdecomp.hpp
@@ -15,41 +15,38 @@
 // 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
+#ifndef ACTION_QRDECOMP
+#define ACTION_QRDECOMP
 
 #include "LinearCongruential.hpp"
 #include <vector>
 #include <algorithm>
 
 template<class Interface>
-class Action_Choleskydecomp{
+class Action_QRdecomp{
 
     typedef typename Interface::Scalar Scalar;
     typedef std::vector<Scalar> vector_t;
 
 private:
     // Invalidate copy constructor
-    Action_Choleskydecomp(const Action_Choleskydecomp&);
+    Action_QRdecomp(const Action_QRdecomp&);
 
 public:
 
     // Constructor
-    Action_Choleskydecomp(int size)
+    Action_QRdecomp(int size)
     : _size(size), lc(10),
-      A(lc.fillVector<Scalar>(size*size)), A_work(size*size)
+      A(lc.fillVector<Scalar>(size*size)), A_work(size*size),
+      tau_work(size), jpiv(size, 0), jpiv_work(size)
     {
-        MESSAGE("Action_Choleskydecomp Constructor");
-
-        // Make A positive definite
-        for (int i = 0; i < size; ++i)
-            A[i+i*size] += 1.2*size;
+        MESSAGE("Action_QRdecomp Constructor");
     }
 
     // Action name
     static std::string name()
     {
-        return "Choleskydecomp_" + Interface::name();
+        return "QRdecomp_" + Interface::name();
     }
 
     double fpo() {
@@ -61,7 +58,10 @@ public:
     }
 
     inline void calculate() {
-        Interface::Choleskydecomp('U', _size, &A_work[0]);
+        // It is crucial that jpiv_work is filled with zero
+        std::copy(jpiv.begin(), jpiv.end(), jpiv_work.begin());
+
+        Interface::QRdecomp(_size, &A_work[0], &jpiv_work[0], &tau_work[0]);
     }
 
     Scalar getResidual() {
@@ -72,10 +72,10 @@ private:
     LinearCongruential<> lc;
 
     const vector_t A;
-    vector_t A_work, eye_work;;
-    std::vector<int> ipiv;
+    vector_t A_work, tau_work;
+    std::vector<int> jpiv, jpiv_work;
 
 };
 
-#endif // ACTION_CHOLESKYDECOMP
+#endif // ACTION_QRDECOMP
 

diff --git a/btl/actions/actionsLAPACK.hpp b/btl/actions/actionsLAPACK.hpp
index 11eb379..12c4180 100644
--- a/btl/actions/actionsLAPACK.hpp
+++ b/btl/actions/actionsLAPACK.hpp
@@ -20,3 +20,4 @@
 
 #include "LAPACK/action_LUdecomp.hpp"
 #include "LAPACK/action_Choleskydecomp.hpp"
+#include "LAPACK/action_QRdecomp.hpp"