Jaesung Oh / Eigen

Eigne Matrix Class Library

Dependents:   MPC_current_control HydraulicControlBoard_SW AHRS Test_ekf ... more

src/QR/HouseholderQR.h@0:13a5d365ba16, 2016-10-13 (annotated)

Committer:
ykuroda
Date:
Thu Oct 13 04:07:23 2016 +0000
Revision:
0:13a5d365ba16
First commint, Eigne Matrix Class Library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
ykuroda 0:13a5d365ba16 6 // Copyright (C) 2010 Vincent Lejeune
ykuroda 0:13a5d365ba16 7 //
ykuroda 0:13a5d365ba16 8 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 9 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 11
ykuroda 0:13a5d365ba16 12 #ifndef EIGEN_QR_H
ykuroda 0:13a5d365ba16 13 #define EIGEN_QR_H
ykuroda 0:13a5d365ba16 14
ykuroda 0:13a5d365ba16 15 namespace Eigen {
ykuroda 0:13a5d365ba16 16
ykuroda 0:13a5d365ba16 17 /** \ingroup QR_Module
ykuroda 0:13a5d365ba16 18 *
ykuroda 0:13a5d365ba16 19 *
ykuroda 0:13a5d365ba16 20 * \class HouseholderQR
ykuroda 0:13a5d365ba16 21 *
ykuroda 0:13a5d365ba16 22 * \brief Householder QR decomposition of a matrix
ykuroda 0:13a5d365ba16 23 *
ykuroda 0:13a5d365ba16 24 * \param MatrixType the type of the matrix of which we are computing the QR decomposition
ykuroda 0:13a5d365ba16 25 *
ykuroda 0:13a5d365ba16 26 * This class performs a QR decomposition of a matrix \b A into matrices \b Q and \b R
ykuroda 0:13a5d365ba16 27 * such that
ykuroda 0:13a5d365ba16 28 * \f[
ykuroda 0:13a5d365ba16 29 * \mathbf{A} = \mathbf{Q} \, \mathbf{R}
ykuroda 0:13a5d365ba16 30 * \f]
ykuroda 0:13a5d365ba16 31 * by using Householder transformations. Here, \b Q a unitary matrix and \b R an upper triangular matrix.
ykuroda 0:13a5d365ba16 32 * The result is stored in a compact way compatible with LAPACK.
ykuroda 0:13a5d365ba16 33 *
ykuroda 0:13a5d365ba16 34 * Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
ykuroda 0:13a5d365ba16 35 * If you want that feature, use FullPivHouseholderQR or ColPivHouseholderQR instead.
ykuroda 0:13a5d365ba16 36 *
ykuroda 0:13a5d365ba16 37 * This Householder QR decomposition is faster, but less numerically stable and less feature-full than
ykuroda 0:13a5d365ba16 38 * FullPivHouseholderQR or ColPivHouseholderQR.
ykuroda 0:13a5d365ba16 39 *
ykuroda 0:13a5d365ba16 40 * \sa MatrixBase::householderQr()
ykuroda 0:13a5d365ba16 41 */
ykuroda 0:13a5d365ba16 42 template<typename _MatrixType> class HouseholderQR
ykuroda 0:13a5d365ba16 43 {
ykuroda 0:13a5d365ba16 44 public:
ykuroda 0:13a5d365ba16 45
ykuroda 0:13a5d365ba16 46 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 47 enum {
ykuroda 0:13a5d365ba16 48 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 49 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 50 Options = MatrixType::Options,
ykuroda 0:13a5d365ba16 51 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 52 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 53 };
ykuroda 0:13a5d365ba16 54 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 55 typedef typename MatrixType::RealScalar RealScalar;
ykuroda 0:13a5d365ba16 56 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 57 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
ykuroda 0:13a5d365ba16 58 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
ykuroda 0:13a5d365ba16 59 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
ykuroda 0:13a5d365ba16 60 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
ykuroda 0:13a5d365ba16 61
ykuroda 0:13a5d365ba16 62 /**
ykuroda 0:13a5d365ba16 63 * \brief Default Constructor.
ykuroda 0:13a5d365ba16 64 *
ykuroda 0:13a5d365ba16 65 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 66 * perform decompositions via HouseholderQR::compute(const MatrixType&).
ykuroda 0:13a5d365ba16 67 */
ykuroda 0:13a5d365ba16 68 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
ykuroda 0:13a5d365ba16 69
ykuroda 0:13a5d365ba16 70 /** \brief Default Constructor with memory preallocation
ykuroda 0:13a5d365ba16 71 *
ykuroda 0:13a5d365ba16 72 * Like the default constructor but with preallocation of the internal data
ykuroda 0:13a5d365ba16 73 * according to the specified problem \a size.
ykuroda 0:13a5d365ba16 74 * \sa HouseholderQR()
ykuroda 0:13a5d365ba16 75 */
ykuroda 0:13a5d365ba16 76 HouseholderQR(Index rows, Index cols)
ykuroda 0:13a5d365ba16 77 : m_qr(rows, cols),
ykuroda 0:13a5d365ba16 78 m_hCoeffs((std::min)(rows,cols)),
ykuroda 0:13a5d365ba16 79 m_temp(cols),
ykuroda 0:13a5d365ba16 80 m_isInitialized(false) {}
ykuroda 0:13a5d365ba16 81
ykuroda 0:13a5d365ba16 82 /** \brief Constructs a QR factorization from a given matrix
ykuroda 0:13a5d365ba16 83 *
ykuroda 0:13a5d365ba16 84 * This constructor computes the QR factorization of the matrix \a matrix by calling
ykuroda 0:13a5d365ba16 85 * the method compute(). It is a short cut for:
ykuroda 0:13a5d365ba16 86 *
ykuroda 0:13a5d365ba16 87 * \code
ykuroda 0:13a5d365ba16 88 * HouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
ykuroda 0:13a5d365ba16 89 * qr.compute(matrix);
ykuroda 0:13a5d365ba16 90 * \endcode
ykuroda 0:13a5d365ba16 91 *
ykuroda 0:13a5d365ba16 92 * \sa compute()
ykuroda 0:13a5d365ba16 93 */
ykuroda 0:13a5d365ba16 94 HouseholderQR(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 95 : m_qr(matrix.rows(), matrix.cols()),
ykuroda 0:13a5d365ba16 96 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
ykuroda 0:13a5d365ba16 97 m_temp(matrix.cols()),
ykuroda 0:13a5d365ba16 98 m_isInitialized(false)
ykuroda 0:13a5d365ba16 99 {
ykuroda 0:13a5d365ba16 100 compute(matrix);
ykuroda 0:13a5d365ba16 101 }
ykuroda 0:13a5d365ba16 102
ykuroda 0:13a5d365ba16 103 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
ykuroda 0:13a5d365ba16 104 * *this is the QR decomposition, if any exists.
ykuroda 0:13a5d365ba16 105 *
ykuroda 0:13a5d365ba16 106 * \param b the right-hand-side of the equation to solve.
ykuroda 0:13a5d365ba16 107 *
ykuroda 0:13a5d365ba16 108 * \returns a solution.
ykuroda 0:13a5d365ba16 109 *
ykuroda 0:13a5d365ba16 110 * \note The case where b is a matrix is not yet implemented. Also, this
ykuroda 0:13a5d365ba16 111 * code is space inefficient.
ykuroda 0:13a5d365ba16 112 *
ykuroda 0:13a5d365ba16 113 * \note_about_checking_solutions
ykuroda 0:13a5d365ba16 114 *
ykuroda 0:13a5d365ba16 115 * \note_about_arbitrary_choice_of_solution
ykuroda 0:13a5d365ba16 116 *
ykuroda 0:13a5d365ba16 117 * Example: \include HouseholderQR_solve.cpp
ykuroda 0:13a5d365ba16 118 * Output: \verbinclude HouseholderQR_solve.out
ykuroda 0:13a5d365ba16 119 */
ykuroda 0:13a5d365ba16 120 template<typename Rhs>
ykuroda 0:13a5d365ba16 121 inline const internal::solve_retval<HouseholderQR, Rhs>
ykuroda 0:13a5d365ba16 122 solve(const MatrixBase<Rhs>& b) const
ykuroda 0:13a5d365ba16 123 {
ykuroda 0:13a5d365ba16 124 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
ykuroda 0:13a5d365ba16 125 return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
ykuroda 0:13a5d365ba16 126 }
ykuroda 0:13a5d365ba16 127
ykuroda 0:13a5d365ba16 128 /** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations.
ykuroda 0:13a5d365ba16 129 *
ykuroda 0:13a5d365ba16 130 * The returned expression can directly be used to perform matrix products. It can also be assigned to a dense Matrix object.
ykuroda 0:13a5d365ba16 131 * Here is an example showing how to recover the full or thin matrix Q, as well as how to perform matrix products using operator*:
ykuroda 0:13a5d365ba16 132 *
ykuroda 0:13a5d365ba16 133 * Example: \include HouseholderQR_householderQ.cpp
ykuroda 0:13a5d365ba16 134 * Output: \verbinclude HouseholderQR_householderQ.out
ykuroda 0:13a5d365ba16 135 */
ykuroda 0:13a5d365ba16 136 HouseholderSequenceType householderQ() const
ykuroda 0:13a5d365ba16 137 {
ykuroda 0:13a5d365ba16 138 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
ykuroda 0:13a5d365ba16 139 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
ykuroda 0:13a5d365ba16 140 }
ykuroda 0:13a5d365ba16 141
ykuroda 0:13a5d365ba16 142 /** \returns a reference to the matrix where the Householder QR decomposition is stored
ykuroda 0:13a5d365ba16 143 * in a LAPACK-compatible way.
ykuroda 0:13a5d365ba16 144 */
ykuroda 0:13a5d365ba16 145 const MatrixType& matrixQR() const
ykuroda 0:13a5d365ba16 146 {
ykuroda 0:13a5d365ba16 147 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
ykuroda 0:13a5d365ba16 148 return m_qr;
ykuroda 0:13a5d365ba16 149 }
ykuroda 0:13a5d365ba16 150
ykuroda 0:13a5d365ba16 151 HouseholderQR& compute(const MatrixType& matrix);
ykuroda 0:13a5d365ba16 152
ykuroda 0:13a5d365ba16 153 /** \returns the absolute value of the determinant of the matrix of which
ykuroda 0:13a5d365ba16 154 * *this is the QR decomposition. It has only linear complexity
ykuroda 0:13a5d365ba16 155 * (that is, O(n) where n is the dimension of the square matrix)
ykuroda 0:13a5d365ba16 156 * as the QR decomposition has already been computed.
ykuroda 0:13a5d365ba16 157 *
ykuroda 0:13a5d365ba16 158 * \note This is only for square matrices.
ykuroda 0:13a5d365ba16 159 *
ykuroda 0:13a5d365ba16 160 * \warning a determinant can be very big or small, so for matrices
ykuroda 0:13a5d365ba16 161 * of large enough dimension, there is a risk of overflow/underflow.
ykuroda 0:13a5d365ba16 162 * One way to work around that is to use logAbsDeterminant() instead.
ykuroda 0:13a5d365ba16 163 *
ykuroda 0:13a5d365ba16 164 * \sa logAbsDeterminant(), MatrixBase::determinant()
ykuroda 0:13a5d365ba16 165 */
ykuroda 0:13a5d365ba16 166 typename MatrixType::RealScalar absDeterminant() const;
ykuroda 0:13a5d365ba16 167
ykuroda 0:13a5d365ba16 168 /** \returns the natural log of the absolute value of the determinant of the matrix of which
ykuroda 0:13a5d365ba16 169 * *this is the QR decomposition. It has only linear complexity
ykuroda 0:13a5d365ba16 170 * (that is, O(n) where n is the dimension of the square matrix)
ykuroda 0:13a5d365ba16 171 * as the QR decomposition has already been computed.
ykuroda 0:13a5d365ba16 172 *
ykuroda 0:13a5d365ba16 173 * \note This is only for square matrices.
ykuroda 0:13a5d365ba16 174 *
ykuroda 0:13a5d365ba16 175 * \note This method is useful to work around the risk of overflow/underflow that's inherent
ykuroda 0:13a5d365ba16 176 * to determinant computation.
ykuroda 0:13a5d365ba16 177 *
ykuroda 0:13a5d365ba16 178 * \sa absDeterminant(), MatrixBase::determinant()
ykuroda 0:13a5d365ba16 179 */
ykuroda 0:13a5d365ba16 180 typename MatrixType::RealScalar logAbsDeterminant() const;
ykuroda 0:13a5d365ba16 181
ykuroda 0:13a5d365ba16 182 inline Index rows() const { return m_qr.rows(); }
ykuroda 0:13a5d365ba16 183 inline Index cols() const { return m_qr.cols(); }
ykuroda 0:13a5d365ba16 184
ykuroda 0:13a5d365ba16 185 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
ykuroda 0:13a5d365ba16 186 *
ykuroda 0:13a5d365ba16 187 * For advanced uses only.
ykuroda 0:13a5d365ba16 188 */
ykuroda 0:13a5d365ba16 189 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
ykuroda 0:13a5d365ba16 190
ykuroda 0:13a5d365ba16 191 protected:
ykuroda 0:13a5d365ba16 192
ykuroda 0:13a5d365ba16 193 static void check_template_parameters()
ykuroda 0:13a5d365ba16 194 {
ykuroda 0:13a5d365ba16 195 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
ykuroda 0:13a5d365ba16 196 }
ykuroda 0:13a5d365ba16 197
ykuroda 0:13a5d365ba16 198 MatrixType m_qr;
ykuroda 0:13a5d365ba16 199 HCoeffsType m_hCoeffs;
ykuroda 0:13a5d365ba16 200 RowVectorType m_temp;
ykuroda 0:13a5d365ba16 201 bool m_isInitialized;
ykuroda 0:13a5d365ba16 202 };
ykuroda 0:13a5d365ba16 203
ykuroda 0:13a5d365ba16 204 template<typename MatrixType>
ykuroda 0:13a5d365ba16 205 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
ykuroda 0:13a5d365ba16 206 {
ykuroda 0:13a5d365ba16 207 using std::abs;
ykuroda 0:13a5d365ba16 208 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
ykuroda 0:13a5d365ba16 209 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
ykuroda 0:13a5d365ba16 210 return abs(m_qr.diagonal().prod());
ykuroda 0:13a5d365ba16 211 }
ykuroda 0:13a5d365ba16 212
ykuroda 0:13a5d365ba16 213 template<typename MatrixType>
ykuroda 0:13a5d365ba16 214 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
ykuroda 0:13a5d365ba16 215 {
ykuroda 0:13a5d365ba16 216 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
ykuroda 0:13a5d365ba16 217 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
ykuroda 0:13a5d365ba16 218 return m_qr.diagonal().cwiseAbs().array().log().sum();
ykuroda 0:13a5d365ba16 219 }
ykuroda 0:13a5d365ba16 220
ykuroda 0:13a5d365ba16 221 namespace internal {
ykuroda 0:13a5d365ba16 222
ykuroda 0:13a5d365ba16 223 /** \internal */
ykuroda 0:13a5d365ba16 224 template<typename MatrixQR, typename HCoeffs>
ykuroda 0:13a5d365ba16 225 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
ykuroda 0:13a5d365ba16 226 {
ykuroda 0:13a5d365ba16 227 typedef typename MatrixQR::Index Index;
ykuroda 0:13a5d365ba16 228 typedef typename MatrixQR::Scalar Scalar;
ykuroda 0:13a5d365ba16 229 typedef typename MatrixQR::RealScalar RealScalar;
ykuroda 0:13a5d365ba16 230 Index rows = mat.rows();
ykuroda 0:13a5d365ba16 231 Index cols = mat.cols();
ykuroda 0:13a5d365ba16 232 Index size = (std::min)(rows,cols);
ykuroda 0:13a5d365ba16 233
ykuroda 0:13a5d365ba16 234 eigen_assert(hCoeffs.size() == size);
ykuroda 0:13a5d365ba16 235
ykuroda 0:13a5d365ba16 236 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
ykuroda 0:13a5d365ba16 237 TempType tempVector;
ykuroda 0:13a5d365ba16 238 if(tempData==0)
ykuroda 0:13a5d365ba16 239 {
ykuroda 0:13a5d365ba16 240 tempVector.resize(cols);
ykuroda 0:13a5d365ba16 241 tempData = tempVector.data();
ykuroda 0:13a5d365ba16 242 }
ykuroda 0:13a5d365ba16 243
ykuroda 0:13a5d365ba16 244 for(Index k = 0; k < size; ++k)
ykuroda 0:13a5d365ba16 245 {
ykuroda 0:13a5d365ba16 246 Index remainingRows = rows - k;
ykuroda 0:13a5d365ba16 247 Index remainingCols = cols - k - 1;
ykuroda 0:13a5d365ba16 248
ykuroda 0:13a5d365ba16 249 RealScalar beta;
ykuroda 0:13a5d365ba16 250 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
ykuroda 0:13a5d365ba16 251 mat.coeffRef(k,k) = beta;
ykuroda 0:13a5d365ba16 252
ykuroda 0:13a5d365ba16 253 // apply H to remaining part of m_qr from the left
ykuroda 0:13a5d365ba16 254 mat.bottomRightCorner(remainingRows, remainingCols)
ykuroda 0:13a5d365ba16 255 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
ykuroda 0:13a5d365ba16 256 }
ykuroda 0:13a5d365ba16 257 }
ykuroda 0:13a5d365ba16 258
ykuroda 0:13a5d365ba16 259 /** \internal */
ykuroda 0:13a5d365ba16 260 template<typename MatrixQR, typename HCoeffs,
ykuroda 0:13a5d365ba16 261 typename MatrixQRScalar = typename MatrixQR::Scalar,
ykuroda 0:13a5d365ba16 262 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
ykuroda 0:13a5d365ba16 263 struct householder_qr_inplace_blocked
ykuroda 0:13a5d365ba16 264 {
ykuroda 0:13a5d365ba16 265 // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
ykuroda 0:13a5d365ba16 266 static void run(MatrixQR& mat, HCoeffs& hCoeffs,
ykuroda 0:13a5d365ba16 267 typename MatrixQR::Index maxBlockSize=32,
ykuroda 0:13a5d365ba16 268 typename MatrixQR::Scalar* tempData = 0)
ykuroda 0:13a5d365ba16 269 {
ykuroda 0:13a5d365ba16 270 typedef typename MatrixQR::Index Index;
ykuroda 0:13a5d365ba16 271 typedef typename MatrixQR::Scalar Scalar;
ykuroda 0:13a5d365ba16 272 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
ykuroda 0:13a5d365ba16 273
ykuroda 0:13a5d365ba16 274 Index rows = mat.rows();
ykuroda 0:13a5d365ba16 275 Index cols = mat.cols();
ykuroda 0:13a5d365ba16 276 Index size = (std::min)(rows, cols);
ykuroda 0:13a5d365ba16 277
ykuroda 0:13a5d365ba16 278 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
ykuroda 0:13a5d365ba16 279 TempType tempVector;
ykuroda 0:13a5d365ba16 280 if(tempData==0)
ykuroda 0:13a5d365ba16 281 {
ykuroda 0:13a5d365ba16 282 tempVector.resize(cols);
ykuroda 0:13a5d365ba16 283 tempData = tempVector.data();
ykuroda 0:13a5d365ba16 284 }
ykuroda 0:13a5d365ba16 285
ykuroda 0:13a5d365ba16 286 Index blockSize = (std::min)(maxBlockSize,size);
ykuroda 0:13a5d365ba16 287
ykuroda 0:13a5d365ba16 288 Index k = 0;
ykuroda 0:13a5d365ba16 289 for (k = 0; k < size; k += blockSize)
ykuroda 0:13a5d365ba16 290 {
ykuroda 0:13a5d365ba16 291 Index bs = (std::min)(size-k,blockSize); // actual size of the block
ykuroda 0:13a5d365ba16 292 Index tcols = cols - k - bs; // trailing columns
ykuroda 0:13a5d365ba16 293 Index brows = rows-k; // rows of the block
ykuroda 0:13a5d365ba16 294
ykuroda 0:13a5d365ba16 295 // partition the matrix:
ykuroda 0:13a5d365ba16 296 // A00 | A01 | A02
ykuroda 0:13a5d365ba16 297 // mat = A10 | A11 | A12
ykuroda 0:13a5d365ba16 298 // A20 | A21 | A22
ykuroda 0:13a5d365ba16 299 // and performs the qr dec of [A11^T A12^T]^T
ykuroda 0:13a5d365ba16 300 // and update [A21^T A22^T]^T using level 3 operations.
ykuroda 0:13a5d365ba16 301 // Finally, the algorithm continue on A22
ykuroda 0:13a5d365ba16 302
ykuroda 0:13a5d365ba16 303 BlockType A11_21 = mat.block(k,k,brows,bs);
ykuroda 0:13a5d365ba16 304 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
ykuroda 0:13a5d365ba16 305
ykuroda 0:13a5d365ba16 306 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
ykuroda 0:13a5d365ba16 307
ykuroda 0:13a5d365ba16 308 if(tcols)
ykuroda 0:13a5d365ba16 309 {
ykuroda 0:13a5d365ba16 310 BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
ykuroda 0:13a5d365ba16 311 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
ykuroda 0:13a5d365ba16 312 }
ykuroda 0:13a5d365ba16 313 }
ykuroda 0:13a5d365ba16 314 }
ykuroda 0:13a5d365ba16 315 };
ykuroda 0:13a5d365ba16 316
ykuroda 0:13a5d365ba16 317 template<typename _MatrixType, typename Rhs>
ykuroda 0:13a5d365ba16 318 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
ykuroda 0:13a5d365ba16 319 : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
ykuroda 0:13a5d365ba16 320 {
ykuroda 0:13a5d365ba16 321 EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
ykuroda 0:13a5d365ba16 322
ykuroda 0:13a5d365ba16 323 template<typename Dest> void evalTo(Dest& dst) const
ykuroda 0:13a5d365ba16 324 {
ykuroda 0:13a5d365ba16 325 const Index rows = dec().rows(), cols = dec().cols();
ykuroda 0:13a5d365ba16 326 const Index rank = (std::min)(rows, cols);
ykuroda 0:13a5d365ba16 327 eigen_assert(rhs().rows() == rows);
ykuroda 0:13a5d365ba16 328
ykuroda 0:13a5d365ba16 329 typename Rhs::PlainObject c(rhs());
ykuroda 0:13a5d365ba16 330
ykuroda 0:13a5d365ba16 331 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
ykuroda 0:13a5d365ba16 332 c.applyOnTheLeft(householderSequence(
ykuroda 0:13a5d365ba16 333 dec().matrixQR().leftCols(rank),
ykuroda 0:13a5d365ba16 334 dec().hCoeffs().head(rank)).transpose()
ykuroda 0:13a5d365ba16 335 );
ykuroda 0:13a5d365ba16 336
ykuroda 0:13a5d365ba16 337 dec().matrixQR()
ykuroda 0:13a5d365ba16 338 .topLeftCorner(rank, rank)
ykuroda 0:13a5d365ba16 339 .template triangularView<Upper>()
ykuroda 0:13a5d365ba16 340 .solveInPlace(c.topRows(rank));
ykuroda 0:13a5d365ba16 341
ykuroda 0:13a5d365ba16 342 dst.topRows(rank) = c.topRows(rank);
ykuroda 0:13a5d365ba16 343 dst.bottomRows(cols-rank).setZero();
ykuroda 0:13a5d365ba16 344 }
ykuroda 0:13a5d365ba16 345 };
ykuroda 0:13a5d365ba16 346
ykuroda 0:13a5d365ba16 347 } // end namespace internal
ykuroda 0:13a5d365ba16 348
ykuroda 0:13a5d365ba16 349 /** Performs the QR factorization of the given matrix \a matrix. The result of
ykuroda 0:13a5d365ba16 350 * the factorization is stored into \c *this, and a reference to \c *this
ykuroda 0:13a5d365ba16 351 * is returned.
ykuroda 0:13a5d365ba16 352 *
ykuroda 0:13a5d365ba16 353 * \sa class HouseholderQR, HouseholderQR(const MatrixType&)
ykuroda 0:13a5d365ba16 354 */
ykuroda 0:13a5d365ba16 355 template<typename MatrixType>
ykuroda 0:13a5d365ba16 356 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 357 {
ykuroda 0:13a5d365ba16 358 check_template_parameters();
ykuroda 0:13a5d365ba16 359
ykuroda 0:13a5d365ba16 360 Index rows = matrix.rows();
ykuroda 0:13a5d365ba16 361 Index cols = matrix.cols();
ykuroda 0:13a5d365ba16 362 Index size = (std::min)(rows,cols);
ykuroda 0:13a5d365ba16 363
ykuroda 0:13a5d365ba16 364 m_qr = matrix;
ykuroda 0:13a5d365ba16 365 m_hCoeffs.resize(size);
ykuroda 0:13a5d365ba16 366
ykuroda 0:13a5d365ba16 367 m_temp.resize(cols);
ykuroda 0:13a5d365ba16 368
ykuroda 0:13a5d365ba16 369 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
ykuroda 0:13a5d365ba16 370
ykuroda 0:13a5d365ba16 371 m_isInitialized = true;
ykuroda 0:13a5d365ba16 372 return *this;
ykuroda 0:13a5d365ba16 373 }
ykuroda 0:13a5d365ba16 374
ykuroda 0:13a5d365ba16 375 /** \return the Householder QR decomposition of \c *this.
ykuroda 0:13a5d365ba16 376 *
ykuroda 0:13a5d365ba16 377 * \sa class HouseholderQR
ykuroda 0:13a5d365ba16 378 */
ykuroda 0:13a5d365ba16 379 template<typename Derived>
ykuroda 0:13a5d365ba16 380 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
ykuroda 0:13a5d365ba16 381 MatrixBase<Derived>::householderQr() const
ykuroda 0:13a5d365ba16 382 {
ykuroda 0:13a5d365ba16 383 return HouseholderQR<PlainObject>(eval());
ykuroda 0:13a5d365ba16 384 }
ykuroda 0:13a5d365ba16 385
ykuroda 0:13a5d365ba16 386 } // end namespace Eigen
ykuroda 0:13a5d365ba16 387
ykuroda 0:13a5d365ba16 388 #endif // EIGEN_QR_H