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 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 template<typename _MatrixType> class FullPivHouseholderQR
00051 {
00052 public:
00053
00054 typedef _MatrixType MatrixType;
00055 enum {
00056 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00057 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00058 Options = MatrixType::Options,
00059 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00060 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00061 };
00062 typedef typename MatrixType::Scalar Scalar;
00063 typedef typename MatrixType::RealScalar RealScalar;
00064 typedef typename MatrixType::Index Index;
00065 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00066 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00067 typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00068 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00069 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00070 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00071 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00072
00073
00074
00075
00076
00077
00078 FullPivHouseholderQR()
00079 : m_qr(),
00080 m_hCoeffs(),
00081 m_rows_transpositions(),
00082 m_cols_transpositions(),
00083 m_cols_permutation(),
00084 m_temp(),
00085 m_isInitialized(false) {}
00086
00087
00088
00089
00090
00091
00092
00093 FullPivHouseholderQR(Index rows, Index cols)
00094 : m_qr(rows, cols),
00095 m_hCoeffs(std::min(rows,cols)),
00096 m_rows_transpositions(rows),
00097 m_cols_transpositions(cols),
00098 m_cols_permutation(cols),
00099 m_temp(std::min(rows,cols)),
00100 m_isInitialized(false) {}
00101
00102 FullPivHouseholderQR(const MatrixType& matrix)
00103 : m_qr(matrix.rows(), matrix.cols()),
00104 m_hCoeffs(std::min(matrix.rows(), matrix.cols())),
00105 m_rows_transpositions(matrix.rows()),
00106 m_cols_transpositions(matrix.cols()),
00107 m_cols_permutation(matrix.cols()),
00108 m_temp(std::min(matrix.rows(), matrix.cols())),
00109 m_isInitialized(false)
00110 {
00111 compute(matrix);
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 template<typename Rhs>
00132 inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00133 solve(const MatrixBase<Rhs>& b) const
00134 {
00135 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00136 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00137 }
00138
00139 MatrixQType matrixQ(void) const;
00140
00141
00142
00143 const MatrixType& matrixQR() const
00144 {
00145 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00146 return m_qr;
00147 }
00148
00149 FullPivHouseholderQR& compute(const MatrixType& matrix);
00150
00151 const PermutationType& colsPermutation() const
00152 {
00153 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00154 return m_cols_permutation;
00155 }
00156
00157 const IntColVectorType& rowsTranspositions() const
00158 {
00159 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00160 return m_rows_transpositions;
00161 }
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176 typename MatrixType::RealScalar absDeterminant() const;
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190 typename MatrixType::RealScalar logAbsDeterminant() const;
00191
00192
00193
00194
00195
00196
00197 inline Index rank() const
00198 {
00199 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00200 return m_rank;
00201 }
00202
00203
00204
00205
00206
00207
00208 inline Index dimensionOfKernel() const
00209 {
00210 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00211 return m_qr.cols() - m_rank;
00212 }
00213
00214
00215
00216
00217
00218
00219
00220 inline bool isInjective() const
00221 {
00222 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00223 return m_rank == m_qr.cols();
00224 }
00225
00226
00227
00228
00229
00230
00231
00232 inline bool isSurjective() const
00233 {
00234 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00235 return m_rank == m_qr.rows();
00236 }
00237
00238
00239
00240
00241
00242
00243 inline bool isInvertible() const
00244 {
00245 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00246 return isInjective() && isSurjective();
00247 }
00248
00249
00250
00251
00252
00253 inline const
00254 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
00255 inverse() const
00256 {
00257 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00258 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
00259 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00260 }
00261
00262 inline Index rows() const { return m_qr.rows(); }
00263 inline Index cols() const { return m_qr.cols(); }
00264 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00265
00266 protected:
00267 MatrixType m_qr;
00268 HCoeffsType m_hCoeffs;
00269 IntColVectorType m_rows_transpositions;
00270 IntRowVectorType m_cols_transpositions;
00271 PermutationType m_cols_permutation;
00272 RowVectorType m_temp;
00273 bool m_isInitialized;
00274 RealScalar m_precision;
00275 Index m_rank;
00276 Index m_det_pq;
00277 };
00278
00279 template<typename MatrixType>
00280 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00281 {
00282 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00283 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00284 return internal::abs(m_qr.diagonal().prod());
00285 }
00286
00287 template<typename MatrixType>
00288 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00289 {
00290 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00291 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00292 return m_qr.diagonal().cwiseAbs().array().log().sum();
00293 }
00294
00295 template<typename MatrixType>
00296 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00297 {
00298 Index rows = matrix.rows();
00299 Index cols = matrix.cols();
00300 Index size = std::min(rows,cols);
00301 m_rank = size;
00302
00303 m_qr = matrix;
00304 m_hCoeffs.resize(size);
00305
00306 m_temp.resize(cols);
00307
00308 m_precision = NumTraits<Scalar>::epsilon() * size;
00309
00310 m_rows_transpositions.resize(matrix.rows());
00311 m_cols_transpositions.resize(matrix.cols());
00312 Index number_of_transpositions = 0;
00313
00314 RealScalar biggest(0);
00315
00316 for (Index k = 0; k < size; ++k)
00317 {
00318 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00319 RealScalar biggest_in_corner;
00320
00321 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00322 .cwiseAbs()
00323 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00324 row_of_biggest_in_corner += k;
00325 col_of_biggest_in_corner += k;
00326 if(k==0) biggest = biggest_in_corner;
00327
00328
00329 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00330 {
00331 m_rank = k;
00332 for(Index i = k; i < size; i++)
00333 {
00334 m_rows_transpositions.coeffRef(i) = i;
00335 m_cols_transpositions.coeffRef(i) = i;
00336 m_hCoeffs.coeffRef(i) = Scalar(0);
00337 }
00338 break;
00339 }
00340
00341 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00342 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00343 if(k != row_of_biggest_in_corner) {
00344 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00345 ++number_of_transpositions;
00346 }
00347 if(k != col_of_biggest_in_corner) {
00348 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00349 ++number_of_transpositions;
00350 }
00351
00352 RealScalar beta;
00353 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00354 m_qr.coeffRef(k,k) = beta;
00355
00356 m_qr.bottomRightCorner(rows-k, cols-k-1)
00357 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00358 }
00359
00360 m_cols_permutation.setIdentity(cols);
00361 for(Index k = 0; k < size; ++k)
00362 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00363
00364 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00365 m_isInitialized = true;
00366
00367 return *this;
00368 }
00369
00370 namespace internal {
00371
00372 template<typename _MatrixType, typename Rhs>
00373 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00374 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00375 {
00376 EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00377
00378 template<typename Dest> void evalTo(Dest& dst) const
00379 {
00380 const Index rows = dec().rows(), cols = dec().cols();
00381 eigen_assert(rhs().rows() == rows);
00382
00383
00384
00385 if(dec().rank()==0)
00386 {
00387 dst.setZero();
00388 return;
00389 }
00390
00391 typename Rhs::PlainObject c(rhs());
00392
00393 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00394 for (Index k = 0; k < dec().rank(); ++k)
00395 {
00396 Index remainingSize = rows-k;
00397 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00398 c.bottomRightCorner(remainingSize, rhs().cols())
00399 .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00400 dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00401 }
00402
00403 if(!dec().isSurjective())
00404 {
00405
00406 RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
00407 RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00408
00409 const RealScalar m_precision = NumTraits<Scalar>::epsilon() * std::min(rows,cols);
00410
00411 if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00412 return;
00413 }
00414 dec().matrixQR()
00415 .topLeftCorner(dec().rank(), dec().rank())
00416 .template triangularView<Upper>()
00417 .solveInPlace(c.topRows(dec().rank()));
00418
00419 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00420 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00421 }
00422 };
00423
00424 }
00425
00426
00427 template<typename MatrixType>
00428 typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<MatrixType>::matrixQ() const
00429 {
00430 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00431
00432
00433
00434 Index rows = m_qr.rows();
00435 Index cols = m_qr.cols();
00436 Index size = std::min(rows,cols);
00437 MatrixQType res = MatrixQType::Identity(rows, rows);
00438 Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
00439 for (Index k = size-1; k >= 0; k--)
00440 {
00441 res.block(k, k, rows-k, rows-k)
00442 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
00443 res.row(k).swap(res.row(m_rows_transpositions.coeff(k)));
00444 }
00445 return res;
00446 }
00447
00448
00449
00450
00451
00452 template<typename Derived>
00453 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00454 MatrixBase<Derived>::fullPivHouseholderQr() const
00455 {
00456 return FullPivHouseholderQR<PlainObject>(eval());
00457 }
00458
00459 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H