10 #ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
11 #define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
13 #include "./InternalHeaderCheck.h"
18 template <
typename MatrixType_>
19 struct traits<CompleteOrthogonalDecomposition<MatrixType_> >
20 : traits<MatrixType_> {
21 typedef MatrixXpr XprKind;
22 typedef SolverStorage StorageKind;
23 typedef int StorageIndex;
53 :
public SolverBase<CompleteOrthogonalDecomposition<MatrixType_> >
56 typedef MatrixType_ MatrixType;
59 template<
typename Derived>
60 friend struct internal::solve_assertion;
64 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
70 typedef typename internal::plain_row_type<MatrixType, Index>::type
72 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
73 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
76 MatrixType, internal::remove_all_t<
77 typename HCoeffsType::ConjugateReturnType>>
79 typedef typename MatrixType::PlainObject PlainObject;
101 : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
119 template <
typename InputType>
121 : m_cpqr(matrix.rows(), matrix.cols()),
122 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
123 m_temp(matrix.cols())
134 template<
typename InputType>
137 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
138 m_temp(matrix.cols())
143 #ifdef EIGEN_PARSED_BY_DOXYGEN
153 template <
typename Rhs>
164 MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
165 applyZOnTheLeftInPlace<false>(Z);
187 template <
typename InputType>
190 m_cpqr.compute(matrix);
282 eigen_assert(m_cpqr.m_isInitialized &&
"CompleteOrthogonalDecomposition is not initialized.");
286 inline Index rows()
const {
return m_cpqr.rows(); }
287 inline Index cols()
const {
return m_cpqr.cols(); }
301 const HCoeffsType&
zCoeffs()
const {
return m_zCoeffs; }
369 eigen_assert(m_cpqr.m_isInitialized &&
"Decomposition is not initialized.");
373 #ifndef EIGEN_PARSED_BY_DOXYGEN
374 template <
typename RhsType,
typename DstType>
375 void _solve_impl(
const RhsType& rhs, DstType& dst)
const;
377 template<
bool Conjugate,
typename RhsType,
typename DstType>
378 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
382 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
384 template<
bool Transpose_,
typename Rhs>
385 void _check_solve_assertion(
const Rhs& b)
const {
386 EIGEN_ONLY_USED_FOR_DEBUG(b);
387 eigen_assert(m_cpqr.m_isInitialized &&
"CompleteOrthogonalDecomposition is not initialized.");
388 eigen_assert((Transpose_?
derived().cols():
derived().rows())==b.rows() &&
"CompleteOrthogonalDecomposition::solve(): invalid number of rows of the right hand side matrix b");
397 template <
bool Conjugate,
typename Rhs>
402 template <
typename Rhs>
405 ColPivHouseholderQR<MatrixType> m_cpqr;
406 HCoeffsType m_zCoeffs;
407 RowVectorType m_temp;
410 template <
typename MatrixType>
411 typename MatrixType::RealScalar
413 return m_cpqr.absDeterminant();
416 template <
typename MatrixType>
417 typename MatrixType::RealScalar
419 return m_cpqr.logAbsDeterminant();
429 template <
typename MatrixType>
435 const Index rank = m_cpqr.rank();
436 const Index cols = m_cpqr.cols();
437 const Index rows = m_cpqr.rows();
438 m_zCoeffs.resize((std::min)(rows, cols));
453 for (
Index k = rank - 1; k >= 0; --k) {
458 m_cpqr.m_qr.col(k).head(k + 1).swap(
459 m_cpqr.m_qr.col(rank - 1).head(k + 1));
466 .tail(cols - rank + 1)
467 .makeHouseholderInPlace(m_zCoeffs(k), beta);
468 m_cpqr.m_qr(k, rank - 1) = beta;
471 m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
472 .applyHouseholderOnTheRight(
473 m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
478 m_cpqr.m_qr.col(k).head(k + 1).swap(
479 m_cpqr.m_qr.col(rank - 1).head(k + 1));
485 template <
typename MatrixType>
486 template <
bool Conjugate,
typename Rhs>
489 const Index cols = this->cols();
490 const Index nrhs = rhs.cols();
491 const Index rank = this->rank();
493 for (
Index k = rank-1; k >= 0; --k) {
495 rhs.row(k).swap(rhs.row(rank - 1));
497 rhs.middleRows(rank - 1, cols - rank + 1)
498 .applyHouseholderOnTheLeft(
499 matrixQTZ().row(k).tail(cols - rank).transpose().
template conjugateIf<!Conjugate>(), zCoeffs().
template conjugateIf<Conjugate>()(k),
502 rhs.row(k).swap(rhs.row(rank - 1));
507 template <
typename MatrixType>
508 template <
typename Rhs>
511 const Index cols = this->cols();
512 const Index nrhs = rhs.cols();
513 const Index rank = this->rank();
515 for (
Index k = 0; k < rank; ++k) {
517 rhs.row(k).swap(rhs.row(rank - 1));
519 rhs.middleRows(rank - 1, cols - rank + 1)
520 .applyHouseholderOnTheLeft(
521 matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
524 rhs.row(k).swap(rhs.row(rank - 1));
529 #ifndef EIGEN_PARSED_BY_DOXYGEN
530 template <
typename MatrixType_>
531 template <
typename RhsType,
typename DstType>
533 const RhsType& rhs, DstType& dst)
const {
534 const Index rank = this->rank();
541 typename RhsType::PlainObject c(rhs);
542 c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
545 dst.topRows(rank) = matrixT()
546 .topLeftCorner(rank, rank)
547 .template triangularView<Upper>()
548 .solve(c.topRows(rank));
550 const Index cols = this->cols();
554 dst.bottomRows(cols - rank).setZero();
555 applyZAdjointOnTheLeftInPlace(dst);
559 dst = colsPermutation() * dst;
562 template<
typename MatrixType_>
563 template<
bool Conjugate,
typename RhsType,
typename DstType>
564 void CompleteOrthogonalDecomposition<MatrixType_>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
566 const Index rank = this->rank();
573 typename RhsType::PlainObject c(colsPermutation().transpose()*rhs);
576 applyZOnTheLeftInPlace<!Conjugate>(c);
579 matrixT().topLeftCorner(rank, rank)
580 .template triangularView<Upper>()
581 .transpose().template conjugateIf<Conjugate>()
582 .solveInPlace(c.topRows(rank));
584 dst.topRows(rank) = c.topRows(rank);
585 dst.bottomRows(rows()-rank).setZero();
587 dst.applyOnTheLeft(householderQ().setLength(rank).
template conjugateIf<!Conjugate>() );
593 template<
typename MatrixType>
594 struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
595 : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
600 template<
typename DstXprType,
typename MatrixType>
601 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
603 typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
604 typedef Inverse<CodType> SrcXprType;
605 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
607 typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType;
608 dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols()));
615 template <
typename MatrixType>
616 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
618 return m_cpqr.householderQ();
625 template <
typename Derived>
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:336
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:396
bool isInvertible() const
Definition: ColPivHouseholderQR.h:312
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:652
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:191
Index rank() const
Definition: ColPivHouseholderQR.h:257
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:380
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:216
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:405
bool isSurjective() const
Definition: ColPivHouseholderQR.h:300
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:274
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:355
bool isInjective() const
Definition: ColPivHouseholderQR.h:287
Complete orthogonal decomposition (COD) of a matrix.
Definition: CompleteOrthogonalDecomposition.h:54
void applyZAdjointOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:509
Index rank() const
Definition: CompleteOrthogonalDecomposition.h:237
RealScalar threshold() const
Definition: CompleteOrthogonalDecomposition.h:344
const HCoeffsType & zCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:301
Index nonzeroPivots() const
Definition: CompleteOrthogonalDecomposition.h:353
const PermutationType & colsPermutation() const
Definition: CompleteOrthogonalDecomposition.h:196
bool isInvertible() const
Definition: CompleteOrthogonalDecomposition.h:273
ComputationInfo info() const
Reports whether the complete orthogonal decomposition was successful.
Definition: CompleteOrthogonalDecomposition.h:368
bool isInjective() const
Definition: CompleteOrthogonalDecomposition.h:255
CompleteOrthogonalDecomposition & setThreshold(Default_t)
Definition: CompleteOrthogonalDecomposition.h:335
MatrixType matrixZ() const
Definition: CompleteOrthogonalDecomposition.h:163
bool isSurjective() const
Definition: CompleteOrthogonalDecomposition.h:264
const Solve< CompleteOrthogonalDecomposition, Rhs > solve(const MatrixBase< Rhs > &b) const
CompleteOrthogonalDecomposition()
Default Constructor.
Definition: CompleteOrthogonalDecomposition.h:92
const MatrixType & matrixT() const
Definition: CompleteOrthogonalDecomposition.h:185
CompleteOrthogonalDecomposition & setThreshold(const RealScalar &threshold)
Definition: CompleteOrthogonalDecomposition.h:322
CompleteOrthogonalDecomposition(const EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:120
Index dimensionOfKernel() const
Definition: CompleteOrthogonalDecomposition.h:246
const Inverse< CompleteOrthogonalDecomposition > pseudoInverse() const
Definition: CompleteOrthogonalDecomposition.h:280
MatrixType::RealScalar absDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:412
RealScalar maxPivot() const
Definition: CompleteOrthogonalDecomposition.h:358
HouseholderSequenceType householderQ(void) const
Definition: CompleteOrthogonalDecomposition.h:617
MatrixType::RealScalar logAbsDeterminant() const
Definition: CompleteOrthogonalDecomposition.h:418
CompleteOrthogonalDecomposition(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: CompleteOrthogonalDecomposition.h:100
void computeInPlace()
Definition: CompleteOrthogonalDecomposition.h:430
const MatrixType & matrixQTZ() const
Definition: CompleteOrthogonalDecomposition.h:172
CompleteOrthogonalDecomposition(EigenBase< InputType > &matrix)
Constructs a complete orthogonal decomposition from a given matrix.
Definition: CompleteOrthogonalDecomposition.h:135
const HCoeffsType & hCoeffs() const
Definition: CompleteOrthogonalDecomposition.h:294
void applyZOnTheLeftInPlace(Rhs &rhs) const
Definition: CompleteOrthogonalDecomposition.h:487
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const CompleteOrthogonalDecomposition< PlainObject > completeOrthogonalDecomposition() const
Definition: CompleteOrthogonalDecomposition.h:627
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
CompleteOrthogonalDecomposition< MatrixType_ > & derived()
Definition: EigenBase.h:48
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231