Estimation DLL, matrix library: fixed problem in QRDecomposition (two workspaces are needed of different sizes) by imposing in the constructor the size of the matrix to be left-multiplied by Q

mex/sources/estimation/DecisionRules.cc

@@ -40,7 +40,7 @@ DecisionRules::DecisionRules(size_t n_arg, size_t p_arg,
   D(n_fwrd + n_back + 2*n_mixed),
   E(n_fwrd + n_back + 2*n_mixed),
   Z(n_fwrd + n_back + 2*n_mixed),
-  QR(n, n_static),
+  QR(n, n_static, n_back_mixed + n + n_fwrd_mixed),
   GSD(n_fwrd + n_back + 2*n_mixed, qz_criterium),
   LU1(n_fwrd_mixed),
   LU2(n_back_mixed),
@@ -95,7 +95,7 @@ DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (
     mat::col_copy(jacobian, n_back_mixed+zeta_static[i], S, i);
 
   A = MatrixConstView(jacobian, 0, 0, n, n_back_mixed + n + n_fwrd_mixed);
-  QR.computeAndMultByQ(S, "L", "N", A);
+  QR.computeAndLeftMultByQ(S, "N", A);
 
   // Construct matrix D
   for (size_t j = 0; j < n_back_mixed; j++)

mex/sources/estimation/libmat/QRDecomposition.cc

@@ -21,11 +21,13 @@
 
 #include "QRDecomposition.hh"
 
-QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg) :
-  rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), lwork(rows*cols),
-  H(rows), Q2(rows), v(rows)
+QRDecomposition::QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg) :
+  rows(rows_arg), cols(cols_arg), mind(std::min(rows, cols)), cols2(cols2_arg),
+  lwork(rows*cols), lwork2(cols2), H(rows), Q2(rows), v(rows)
 {
   work = new double[lwork];
+  work2 = new double[lwork2];
+
   tau = new double[mind];
 }
 

mex/sources/estimation/libmat/QRDecomposition.hh

@@ -29,42 +29,46 @@ class QRDecomposition
 {
 private:
   const size_t rows, cols, mind;
-  lapack_int lwork;
-  double *work, *tau;
+  //! Number of columns of the matrix to be left-multiplied by Q
+  const size_t cols2;
+  lapack_int lwork, lwork2;
+  double *work, *work2, *tau;
   Matrix H, Q2;
   Vector v;
 public:
   /*!
     \todo Replace heuristic choice for workspace size by a query to determine the optimal size
+    \param[in] rows_arg Number of rows of the matrix to decompose
+    \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);
+  QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg);
   ~QRDecomposition();
-  //! Performs the QR decomposition of a matrix, and multiplies another matrix by Q
+  //! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q
   /*!
     \param[in,out] A On input, the matrix to be decomposed. On output, equals to the output of dgeqrf
-    \param[in] side Specifies on which side is Q in the multiplication, either "L" or "R"
     \param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
-    \param[in,out] C The matrix to be multiplied by Q, modified in place
+    \param[in,out] C The matrix to be left-multiplied by Q, modified in place
   */
   template<class Mat1, class Mat2>
-  void computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C);
+  void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C);
 };
 
 template<class Mat1, class Mat2>
 void
-QRDecomposition::computeAndMultByQ(Mat1 &A, const char *side, const char *trans, Mat2 &C)
+QRDecomposition::computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C)
 {
   assert(A.getRows() == rows && A.getCols() == cols);
-  assert(C.getRows() == rows);
+  assert(C.getRows() == rows && C.getCols() == cols2);
 
   lapack_int m = rows, n = cols, lda = A.getLd();
   lapack_int info;
   dgeqrf(&m, &n, A.getData(), &lda, tau, work, &lwork, &info);
   assert(info == 0);
 
-  n = C.getCols();
+  n = cols2;
   lapack_int k = mind, ldc = C.getLd();
-  dormqr(side, trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
-         work, &lwork, &info);
+  dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
+         work2, &lwork2, &info);
   assert(info == 0);
 }

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

@@ -27,7 +27,7 @@ main(int argc, char **argv)
 {
   size_t m = 4, n = 3;
   Matrix S(m, n), Q(m), A(m, n), B(m);
-  QRDecomposition QRD(m, n);
+  QRDecomposition QRD(m, n, m);
 
   for (size_t i = 0; i < m; i++)
     for (size_t j = 0; j < n; j++)
@@ -37,7 +37,7 @@ main(int argc, char **argv)
 
   mat::set_identity(Q);
 
-  QRD.computeAndMultByQ(S, "L", "N", Q);
+  QRD.computeAndLeftMultByQ(S, "N", Q);
 
   std::cout << "Matrix Q:" << std::endl << Q << std::endl;