* [gentoo-commits] proj/auto-numerical-bench:newinterfaces commit in: btl/actions/, btl/NumericInterface/NI_internal/, btl/actions/LAPACK/
@ 2012-10-09 7:19 Andrea Arteaga
0 siblings, 0 replies; 3+ messages in thread
From: Andrea Arteaga @ 2012-10-09 7:19 UTC (permalink / raw
To: gentoo-commits
commit: 91b7fe13b24657071e06e36a5322344d02a45789
Author: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
AuthorDate: Tue Oct 9 07:18:57 2012 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Tue Oct 9 07:18:57 2012 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=91b7fe13
Fixed GeneralSolve and added LeastSquares.
---
.../NI_internal/FortranDeclarations.hpp | 8 +++---
.../NI_internal/FortranInterface.hpp | 13 ++++++++++
.../{GeneralSolve.hpp => action_GeneralSolve.hpp} | 3 +-
...neralSolve.hpp => action_LeastSquaresSolve.hpp} | 24 ++++++++++----------
btl/actions/actionsLAPACK.hpp | 19 +++++++++++++++
5 files changed, 50 insertions(+), 17 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
similarity index 95%
copy from btl/actions/LAPACK/GeneralSolve.hpp
copy to 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/GeneralSolve.hpp b/btl/actions/LAPACK/action_LeastSquaresSolve.hpp
similarity index 77%
rename from btl/actions/LAPACK/GeneralSolve.hpp
rename to btl/actions/LAPACK/action_LeastSquaresSolve.hpp
index 7fb98ba..4bd7da4 100644
--- a/btl/actions/LAPACK/GeneralSolve.hpp
+++ b/btl/actions/LAPACK/action_LeastSquaresSolve.hpp
@@ -15,38 +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_GENERALSOLVE
-#define ACTION_GENERALSOLVE
+#ifndef ACTION_LEASTSQUARESSOLVE
+#define ACTION_LEASTSQUARESSOLVE
#include "LinearCongruential.hpp"
#include <vector>
#include <algorithm>
template<class Interface>
-class Action_GeneralSolve {
+class Action_LeastSquaresSolve {
typedef typename Interface::Scalar Scalar;
typedef std::vector<Scalar> vector_t;
private:
// Invalidate copy constructor
- Action_GeneralSolve(const Action_GeneralSolve&);
+ Action_LeastSquaresSolve(const Action_LeastSquaresSolve&);
public:
// Constructor
- Action_GeneralSolve(int size)
+ 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), ipiv(size)
+ A_work(size*size), x_work(size), b_res(size)
{
- MESSAGE("Action_GeneralSolve Constructor");
+ MESSAGE("Action_LeastSquaresSolve Constructor");
}
// Action name
static std::string name()
{
- return "GeneralSolve_" + Interface::name();
+ return "LeastSquaresSolve_" + Interface::name();
}
double fpo() {
@@ -59,7 +59,7 @@ public:
}
inline void calculate() {
- Interface::GeneralSolve(_size, &A_work[0], &x_work[0], &ipiv[0]);
+ Interface::LeastSquaresSolve(_size, &A_work[0], &x_work[0]);
}
Scalar getResidual() {
@@ -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:
@@ -77,8 +78,7 @@ private:
const vector_t A, b;
vector_t A_work, x_work, b_res;
- std::vector<int> ipiv;
};
-#endif // ACTION_GENERALSOLVE
+#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"
^ permalink raw reply related [flat|nested] 3+ messages in thread
* [gentoo-commits] proj/auto-numerical-bench:newinterfaces commit in: btl/actions/, btl/NumericInterface/NI_internal/, btl/actions/LAPACK/
@ 2012-10-15 21:14 Andrea Arteaga
0 siblings, 0 replies; 3+ messages in thread
From: Andrea Arteaga @ 2012-10-15 21:14 UTC (permalink / raw
To: gentoo-commits
commit: 93de8d4c28ecd6808ebe798b0a8860239d05a023
Author: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
AuthorDate: Mon Oct 15 21:14:41 2012 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Mon Oct 15 21:14:41 2012 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=93de8d4c
Added LU and Cholesky decompositions (errors not working yet)
---
.../NI_internal/FortranDeclarations.hpp | 13 +++
.../NI_internal/FortranInterface.hpp | 20 ++++
btl/actions/LAPACK/action_Choleskydecomp.hpp | 81 ++++++++++++++++
btl/actions/LAPACK/action_LUdecomp.hpp | 102 ++++++++++++++++++++
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"
^ permalink raw reply related [flat|nested] 3+ messages in thread
* [gentoo-commits] proj/auto-numerical-bench:newinterfaces commit in: btl/actions/, btl/NumericInterface/NI_internal/, btl/actions/LAPACK/
@ 2012-10-15 22:17 Andrea Arteaga
0 siblings, 0 replies; 3+ messages in thread
From: Andrea Arteaga @ 2012-10-15 22:17 UTC (permalink / raw
To: gentoo-commits
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"
^ permalink raw reply related [flat|nested] 3+ messages in thread
end of thread, other threads:[~2012-10-15 22:18 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2012-10-15 21:14 [gentoo-commits] proj/auto-numerical-bench:newinterfaces commit in: btl/actions/, btl/NumericInterface/NI_internal/, btl/actions/LAPACK/ Andrea Arteaga
-- strict thread matches above, loose matches on Subject: below --
2012-10-15 22:17 Andrea Arteaga
2012-10-09 7:19 Andrea Arteaga
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox