diff --git a/mex/sources/estimation/libmat/QRDecomposition.cc b/mex/sources/estimation/libmat/QRDecomposition.cc
index 38e3a5f5718a90d667fa6b6356be34ef85d1eb77..41f48b62c1c3750fcfd12e1954366a66900aff68 100644
--- a/mex/sources/estimation/libmat/QRDecomposition.cc
+++ b/mex/sources/estimation/libmat/QRDecomposition.cc
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -21,7 +21,7 @@
#include "QRDecomposition.hh"
-QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg) :
+QRDecomposition::QRDecomposition(ptrdiff_t rows_arg, ptrdiff_t cols_arg, ptrdiff_t cols2_arg) :
rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), cols2(cols2_arg),
lwork(rows*cols), lwork2(cols2)
{
diff --git a/mex/sources/estimation/libmat/QRDecomposition.hh b/mex/sources/estimation/libmat/QRDecomposition.hh
index d8022a7e5cc66d60fd971f9865a3345b6021081f..30e03eb80e052d2d401f5fab7459c2318df94c44 100644
--- a/mex/sources/estimation/libmat/QRDecomposition.hh
+++ b/mex/sources/estimation/libmat/QRDecomposition.hh
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010-2011 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -19,18 +19,19 @@
#include <algorithm> // For std::min()
-#include <dynlapack.h>
+#include <Eigen/Core>
+
+using namespace Eigen;
-#include "Vector.hh"
-#include "Matrix.hh"
-#include "BlasBindings.hh"
+
+#include <dynlapack.h>
class QRDecomposition
{
private:
- const size_t rows, cols, mind;
+ const ptrdiff_t rows, cols, mind;
//! Number of columns of the matrix to be left-multiplied by Q
- const size_t cols2;
+ const ptrdiff_t cols2;
lapack_int lwork, lwork2;
double *work, *work2, *tau;
public:
@@ -40,7 +41,7 @@ public:
\param[in] cols_arg Number of columns of the matrix to decompose
\param[in] cols2_arg Number of columns of the matrix to be multiplied by Q
*/
- QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg);
+ QRDecomposition(ptrdiff_t rows_arg, ptrdiff_t cols_arg, ptrdiff_t cols2_arg);
virtual ~QRDecomposition();
//! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q
/*!
@@ -48,25 +49,27 @@ public:
\param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
\param[in,out] C The matrix to be left-multiplied by Q, modified in place
*/
- template<class Mat1, class Mat2>
- void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C);
+ template<typename Derived1, typename Derived2>
+ void computeAndLeftMultByQ(PlainObjectBase<Derived1> &A, const char *trans,
+ PlainObjectBase<Derived2> &C);
};
-template<class Mat1, class Mat2>
+template<typename Derived1, typename Derived2>
void
-QRDecomposition::computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C)
+QRDecomposition::computeAndLeftMultByQ(PlainObjectBase<Derived1> &A, const char *trans,
+ PlainObjectBase<Derived2> &C)
{
- assert(A.getRows() == rows && A.getCols() == cols);
- assert(C.getRows() == rows && C.getCols() == cols2);
+ assert(A.rows() == rows && A.cols() == cols);
+ assert(C.rows() == rows && C.cols() == cols2);
- lapack_int m = rows, n = cols, lda = A.getLd();
+ lapack_int m = rows, n = cols, lda = A.outerStride();
lapack_int info;
- dgeqrf(&m, &n, A.getData(), &lda, tau, work, &lwork, &info);
+ dgeqrf(&m, &n, A.data(), &lda, tau, work, &lwork, &info);
assert(info == 0);
n = cols2;
- lapack_int k = mind, ldc = C.getLd();
- dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
+ lapack_int k = mind, ldc = C.outerStride();
+ dormqr("L", trans, &m, &n, &k, A.data(), &lda, tau, C.data(), &ldc,
work2, &lwork2, &info);
assert(info == 0);
}
diff --git a/mex/sources/estimation/libmat/tests/Makefile.am b/mex/sources/estimation/libmat/tests/Makefile.am
index 7223fe01de26281bcaf8c0f19b1454b06c18a31d..33d4a988389e990495db478fc84a649002d8494a 100644
--- a/mex/sources/estimation/libmat/tests/Makefile.am
+++ b/mex/sources/estimation/libmat/tests/Makefile.am
@@ -1,6 +1,6 @@
check_PROGRAMS = test-qr test-gsd test-lu test-repmat
-test_qr_SOURCES = ../Matrix.cc ../Vector.cc ../QRDecomposition.cc test-qr.cc
+test_qr_SOURCES = ../QRDecomposition.cc test-qr.cc
test_qr_LDADD = $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
test_qr_CPPFLAGS = -I.. -I../../../
diff --git a/mex/sources/estimation/libmat/tests/test-qr.cc b/mex/sources/estimation/libmat/tests/test-qr.cc
index f73d85488222d9aa53ff7c8e1af2574d49864f60..b560f3fc49b326777b5bb504f1355f924b30c0a1 100644
--- a/mex/sources/estimation/libmat/tests/test-qr.cc
+++ b/mex/sources/estimation/libmat/tests/test-qr.cc
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010 Dynare Team
+ * Copyright (C) 2010-2012 Dynare Team
*
* This file is part of Dynare.
*
@@ -19,14 +19,13 @@
#include <iostream>
-#include "BlasBindings.hh"
#include "QRDecomposition.hh"
int
main(int argc, char **argv)
{
size_t m = 4, n = 3;
- Matrix S(m, n), Q(m), A(m, n), B(m), S2(m, n);
+ MatrixXd S(m, n), Q(m, m), A(m, n), S2(m, n), B;
QRDecomposition QRD(m, n, m);
for (size_t i = 0; i < m; i++)
@@ -35,32 +34,27 @@ main(int argc, char **argv)
std::cout << "Matrix to be decomposed:" << std::endl << S << std::endl;
- mat::set_identity(Q);
+ Q.setIdentity();
S2 = S;
QRD.computeAndLeftMultByQ(S2, "N", Q);
std::cout << "Q =" << std::endl << Q << std::endl;
- blas::gemm("T", "N", 1.0, Q, Q, 0.0, B);
-
+ B = Q.transpose()*Q;
+
std::cout << "Q'*Q =" << std::endl << B << std::endl;
- for (size_t j = 0; j < n; j++)
- mat::col_set(S2, j, j+1, m-j-1, 0);
+ S2 = S2.triangularView<Upper>();
std::cout << "R =" << std::endl << S2 << std::endl;
- blas::gemm("N", "N", 1.0, Q, S2, 0.0, A);
+ A = Q*S2;
std::cout << "Q*R =" << std::endl << A << std::endl;
// Check values
- Matrix B2(m);
- mat::set_identity(B2);
- mat::sub(B2, B);
- assert(mat::nrminf(B2) < 1e-4);
+ assert((B - MatrixXd::Identity(m,m)).lpNorm<Infinity>() < 1e-4);
- mat::sub(A, S);
- assert(mat::nrminf(A) < 1e-4);
+ assert((A-S).lpNorm<Infinity>() < 1e-4);
}