Convert QRDecomposition to Eigen

mex/sources/estimation/libmat/QRDecomposition.cc

 /*
- * 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)
 {

mex/sources/estimation/libmat/QRDecomposition.hh

 /*
- * 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);
 }

mex/sources/estimation/libmat/tests/Makefile.am

 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../../../
 

mex/sources/estimation/libmat/tests/test-qr.cc

 /*
- * 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);
 }