00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #ifndef EIGEN_QR_H
00028 #define EIGEN_QR_H
00029
00030 namespace Eigen {
00031
00057 template<typename _MatrixType> class HouseholderQR
00058 {
00059 public:
00060
00061 typedef _MatrixType MatrixType;
00062 enum {
00063 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065 Options = MatrixType::Options,
00066 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00067 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00068 };
00069 typedef typename MatrixType::Scalar Scalar;
00070 typedef typename MatrixType::RealScalar RealScalar;
00071 typedef typename MatrixType::Index Index;
00072 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00073 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00074 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00075 typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
00076
00083 HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
00084
00091 HouseholderQR(Index rows, Index cols)
00092 : m_qr(rows, cols),
00093 m_hCoeffs((std::min)(rows,cols)),
00094 m_temp(cols),
00095 m_isInitialized(false) {}
00096
00097 HouseholderQR(const MatrixType& matrix)
00098 : m_qr(matrix.rows(), matrix.cols()),
00099 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00100 m_temp(matrix.cols()),
00101 m_isInitialized(false)
00102 {
00103 compute(matrix);
00104 }
00105
00123 template<typename Rhs>
00124 inline const internal::solve_retval<HouseholderQR, Rhs>
00125 solve(const MatrixBase<Rhs>& b) const
00126 {
00127 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00128 return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
00129 }
00130
00131 HouseholderSequenceType householderQ() const
00132 {
00133 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00134 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00135 }
00136
00140 const MatrixType& matrixQR() const
00141 {
00142 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00143 return m_qr;
00144 }
00145
00146 HouseholderQR& compute(const MatrixType& matrix);
00147
00161 typename MatrixType::RealScalar absDeterminant() const;
00162
00175 typename MatrixType::RealScalar logAbsDeterminant() const;
00176
00177 inline Index rows() const { return m_qr.rows(); }
00178 inline Index cols() const { return m_qr.cols(); }
00179 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00180
00181 protected:
00182 MatrixType m_qr;
00183 HCoeffsType m_hCoeffs;
00184 RowVectorType m_temp;
00185 bool m_isInitialized;
00186 };
00187
00188 template<typename MatrixType>
00189 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00190 {
00191 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00192 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00193 return internal::abs(m_qr.diagonal().prod());
00194 }
00195
00196 template<typename MatrixType>
00197 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00198 {
00199 eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00200 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00201 return m_qr.diagonal().cwiseAbs().array().log().sum();
00202 }
00203
00204 namespace internal {
00205
00207 template<typename MatrixQR, typename HCoeffs>
00208 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00209 {
00210 typedef typename MatrixQR::Index Index;
00211 typedef typename MatrixQR::Scalar Scalar;
00212 typedef typename MatrixQR::RealScalar RealScalar;
00213 Index rows = mat.rows();
00214 Index cols = mat.cols();
00215 Index size = (std::min)(rows,cols);
00216
00217 eigen_assert(hCoeffs.size() == size);
00218
00219 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00220 TempType tempVector;
00221 if(tempData==0)
00222 {
00223 tempVector.resize(cols);
00224 tempData = tempVector.data();
00225 }
00226
00227 for(Index k = 0; k < size; ++k)
00228 {
00229 Index remainingRows = rows - k;
00230 Index remainingCols = cols - k - 1;
00231
00232 RealScalar beta;
00233 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00234 mat.coeffRef(k,k) = beta;
00235
00236
00237 mat.bottomRightCorner(remainingRows, remainingCols)
00238 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00239 }
00240 }
00241
00243 template<typename MatrixQR, typename HCoeffs>
00244 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
00245 typename MatrixQR::Index maxBlockSize=32,
00246 typename MatrixQR::Scalar* tempData = 0)
00247 {
00248 typedef typename MatrixQR::Index Index;
00249 typedef typename MatrixQR::Scalar Scalar;
00250 typedef typename MatrixQR::RealScalar RealScalar;
00251 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00252
00253 Index rows = mat.rows();
00254 Index cols = mat.cols();
00255 Index size = (std::min)(rows, cols);
00256
00257 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00258 TempType tempVector;
00259 if(tempData==0)
00260 {
00261 tempVector.resize(cols);
00262 tempData = tempVector.data();
00263 }
00264
00265 Index blockSize = (std::min)(maxBlockSize,size);
00266
00267 Index k = 0;
00268 for (k = 0; k < size; k += blockSize)
00269 {
00270 Index bs = (std::min)(size-k,blockSize);
00271 Index tcols = cols - k - bs;
00272 Index brows = rows-k;
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282 BlockType A11_21 = mat.block(k,k,brows,bs);
00283 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00284
00285 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00286
00287 if(tcols)
00288 {
00289 BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00290 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
00291 }
00292 }
00293 }
00294
00295 template<typename _MatrixType, typename Rhs>
00296 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
00297 : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
00298 {
00299 EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
00300
00301 template<typename Dest> void evalTo(Dest& dst) const
00302 {
00303 const Index rows = dec().rows(), cols = dec().cols();
00304 const Index rank = (std::min)(rows, cols);
00305 eigen_assert(rhs().rows() == rows);
00306
00307 typename Rhs::PlainObject c(rhs());
00308
00309
00310 c.applyOnTheLeft(householderSequence(
00311 dec().matrixQR().leftCols(rank),
00312 dec().hCoeffs().head(rank)).transpose()
00313 );
00314
00315 dec().matrixQR()
00316 .topLeftCorner(rank, rank)
00317 .template triangularView<Upper>()
00318 .solveInPlace(c.topRows(rank));
00319
00320 dst.topRows(rank) = c.topRows(rank);
00321 dst.bottomRows(cols-rank).setZero();
00322 }
00323 };
00324
00325 }
00326
00327 template<typename MatrixType>
00328 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00329 {
00330 Index rows = matrix.rows();
00331 Index cols = matrix.cols();
00332 Index size = (std::min)(rows,cols);
00333
00334 m_qr = matrix;
00335 m_hCoeffs.resize(size);
00336
00337 m_temp.resize(cols);
00338
00339 internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
00340
00341 m_isInitialized = true;
00342 return *this;
00343 }
00344
00349 template<typename Derived>
00350 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00351 MatrixBase<Derived>::householderQr() const
00352 {
00353 return HouseholderQR<PlainObject>(eval());
00354 }
00355
00356 }
00357
00358 #endif // EIGEN_QR_H