QRDecomposition.hh 2.48 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
/*
 * Copyright (C) 2010 Dynare Team
 *
 * This file is part of Dynare.
 *
 * Dynare 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 3 of the License, or
 * (at your option) any later version.
 *
 * Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <algorithm> // For std::min()

#include <dynlapack.h>

#include "Vector.hh"
#include "Matrix.hh"
#include "BlasBindings.hh"

class QRDecomposition
{
private:
  const size_t rows, cols, mind;
32 33 34 35
  //! Number of columns of the matrix to be left-multiplied by Q
  const size_t cols2;
  lapack_int lwork, lwork2;
  double *work, *work2, *tau;
36 37 38
public:
  /*!
    \todo Replace heuristic choice for workspace size by a query to determine the optimal size
39 40 41
    \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
42
   */
43
  QRDecomposition(size_t rows_arg, size_t cols_arg, size_t cols2_arg);
44
  virtual ~QRDecomposition();
45
  //! Performs the QR decomposition of a matrix, and left-multiplies another matrix by Q
46 47 48
  /*!
    \param[in,out] A On input, the matrix to be decomposed. On output, equals to the output of dgeqrf
    \param[in] trans Specifies whether Q should be transposed before the multiplication, either "T" or "N"
49
    \param[in,out] C The matrix to be left-multiplied by Q, modified in place
50 51
  */
  template<class Mat1, class Mat2>
52
  void computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C);
53 54 55 56
};

template<class Mat1, class Mat2>
void
57
QRDecomposition::computeAndLeftMultByQ(Mat1 &A, const char *trans, Mat2 &C)
58 59
{
  assert(A.getRows() == rows && A.getCols() == cols);
60
  assert(C.getRows() == rows && C.getCols() == cols2);
61 62 63 64 65 66

  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);

67
  n = cols2;
68
  lapack_int k = mind, ldc = C.getLd();
69 70
  dormqr("L", trans, &m, &n, &k, A.getData(), &lda, tau, C.getData(), &ldc,
         work2, &lwork2, &info);
71 72
  assert(info == 0);
}