Bugzilla – Attachment 714 Details for
Bug 707
In-Place decompositions
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
[patch]
Patch for inplace decomposition via SolverName<Ref<Matrix>>
bug_707_inplace_decomposition.diff (text/plain), 16.52 KB, created by
Gael Guennebaud
on 2016-06-23 16:32:39 UTC
(
hide
)
Description:
Patch for inplace decomposition via SolverName<Ref<Matrix>>
Filename:
MIME Type:
Creator:
Gael Guennebaud
Created:
2016-06-23 16:32:39 UTC
Size:
16.52 KB
patch
obsolete
>diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h >--- a/Eigen/src/Cholesky/LDLT.h >+++ b/Eigen/src/Cholesky/LDLT.h >@@ -47,26 +47,25 @@ namespace internal { > */ > template<typename _MatrixType, int _UpLo> class LDLT > { > public: > typedef _MatrixType MatrixType; > enum { > RowsAtCompileTime = MatrixType::RowsAtCompileTime, > ColsAtCompileTime = MatrixType::ColsAtCompileTime, >- Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here! > MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, > MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, > UpLo = _UpLo > }; > typedef typename MatrixType::Scalar Scalar; > typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; > typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 > typedef typename MatrixType::StorageIndex StorageIndex; >- typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType; >+ typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType; > > typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; > typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; > > typedef internal::LDLT_Traits<MatrixType,UpLo> Traits; > > /** \brief Default Constructor. > * >@@ -92,29 +91,47 @@ template<typename _MatrixType, int _UpLo > m_temporary(size), > m_sign(internal::ZeroSign), > m_isInitialized(false) > {} > > /** \brief Constructor with decomposition > * > * This calculates the decomposition for the input \a matrix. >+ * > * \sa LDLT(Index size) > */ > template<typename InputType> > explicit LDLT(const EigenBase<InputType>& matrix) > : m_matrix(matrix.rows(), matrix.cols()), > m_transpositions(matrix.rows()), > m_temporary(matrix.rows()), > m_sign(internal::ZeroSign), > m_isInitialized(false) > { > compute(matrix.derived()); > } > >+ /** \brief Constructs a LDLT factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. >+ * >+ * \sa LDLT(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit LDLT(EigenBase<InputType>& matrix) >+ : m_matrix(matrix.derived()), >+ m_transpositions(matrix.rows()), >+ m_temporary(matrix.rows()), >+ m_sign(internal::ZeroSign), >+ m_isInitialized(false) >+ { >+ compute(matrix.derived()); >+ } >+ > /** Clear any existing decomposition > * \sa rankUpdate(w,sigma) > */ > void setZero() > { > m_isInitialized = false; > } > >diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h >--- a/Eigen/src/Cholesky/LLT.h >+++ b/Eigen/src/Cholesky/LLT.h >@@ -49,17 +49,16 @@ template<typename MatrixType, int UpLo> > */ > template<typename _MatrixType, int _UpLo> class LLT > { > public: > typedef _MatrixType MatrixType; > enum { > RowsAtCompileTime = MatrixType::RowsAtCompileTime, > ColsAtCompileTime = MatrixType::ColsAtCompileTime, >- Options = MatrixType::Options, > MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime > }; > typedef typename MatrixType::Scalar Scalar; > typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; > typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 > typedef typename MatrixType::StorageIndex StorageIndex; > > enum { >@@ -90,16 +89,31 @@ template<typename _MatrixType, int _UpLo > template<typename InputType> > explicit LLT(const EigenBase<InputType>& matrix) > : m_matrix(matrix.rows(), matrix.cols()), > m_isInitialized(false) > { > compute(matrix.derived()); > } > >+ /** \brief Constructs a LDLT factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when >+ * \c MatrixType is a Eigen::Ref. >+ * >+ * \sa LLT(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit LLT(EigenBase<InputType>& matrix) >+ : m_matrix(matrix.derived()), >+ m_isInitialized(false) >+ { >+ compute(matrix.derived()); >+ } >+ > /** \returns a view of the upper triangular matrix U */ > inline typename Traits::MatrixU matrixU() const > { > eigen_assert(m_isInitialized && "LLT is not initialized."); > return Traits::getU(m_matrix); > } > > /** \returns a view of the lower triangular matrix L */ >diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h >--- a/Eigen/src/LU/FullPivLU.h >+++ b/Eigen/src/LU/FullPivLU.h >@@ -92,16 +92,25 @@ template<typename _MatrixType> class Ful > /** Constructor. > * > * \param matrix the matrix of which to compute the LU decomposition. > * It is required to be nonzero. > */ > template<typename InputType> > explicit FullPivLU(const EigenBase<InputType>& matrix); > >+ /** \brief Constructs a LU factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. >+ * >+ * \sa FullPivLU(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit FullPivLU(EigenBase<InputType>& matrix); >+ > /** Computes the LU decomposition of the given matrix. > * > * \param matrix the matrix of which to compute the LU decomposition. > * It is required to be nonzero. > * > * \returns a reference to *this > */ > template<typename InputType> >@@ -454,16 +463,30 @@ FullPivLU<MatrixType>::FullPivLU(const E > m_isInitialized(false), > m_usePrescribedThreshold(false) > { > compute(matrix.derived()); > } > > template<typename MatrixType> > template<typename InputType> >+FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix) >+ : m_lu(matrix.derived()), >+ m_p(matrix.rows()), >+ m_q(matrix.cols()), >+ m_rowsTranspositions(matrix.rows()), >+ m_colsTranspositions(matrix.cols()), >+ m_isInitialized(false), >+ m_usePrescribedThreshold(false) >+{ >+ compute(matrix.derived()); >+} >+ >+template<typename MatrixType> >+template<typename InputType> > FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const EigenBase<InputType>& matrix) > { > check_template_parameters(); > > // the permutations are stored as int indices, so just to be sure: > eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest()); > > m_lu = matrix.derived(); >diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h >--- a/Eigen/src/LU/PartialPivLU.h >+++ b/Eigen/src/LU/PartialPivLU.h >@@ -97,16 +97,33 @@ template<typename _MatrixType> class Par > * \param matrix the matrix of which to compute the LU decomposition. > * > * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). > * If you need to deal with non-full rank, use class FullPivLU instead. > */ > template<typename InputType> > explicit PartialPivLU(const EigenBase<InputType>& matrix); > >+ /** Constructor for inplace decomposition >+ * >+ * \param matrix the matrix of which to compute the LU decomposition. >+ * >+ * If \c MatrixType is an Eigen::Ref, then the storage of \a matrix will be shared >+ * between \a matrix and \c *this and the decomposition will take place in-place. >+ * The memory of \a matrix will be used througrough the lifetime of \c *this. In >+ * particular, further calls to \c this->compute(A) will still operate on the memory >+ * of \a matrix meaning. This also implies that the sizes of \c A must match the >+ * ones of \a matrix. >+ * >+ * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). >+ * If you need to deal with non-full rank, use class FullPivLU instead. >+ */ >+ template<typename InputType> >+ explicit PartialPivLU(EigenBase<InputType>& matrix); >+ > template<typename InputType> > PartialPivLU& compute(const EigenBase<InputType>& matrix); > > /** \returns the LU decomposition matrix: the upper-triangular part is U, the > * unit-lower-triangular part is L (at least for square matrices; in the non-square > * case, special care is needed, see the documentation of class FullPivLU). > * > * \sa matrixL(), matrixU() >@@ -279,17 +296,30 @@ PartialPivLU<MatrixType>::PartialPivLU(I > m_det_p(0), > m_isInitialized(false) > { > } > > template<typename MatrixType> > template<typename InputType> > PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix) >- : m_lu(matrix.rows(), matrix.rows()), >+ : m_lu(matrix.rows(),matrix.cols()), >+ m_p(matrix.rows()), >+ m_rowsTranspositions(matrix.rows()), >+ m_l1_norm(0), >+ m_det_p(0), >+ m_isInitialized(false) >+{ >+ compute(matrix.derived()); >+} >+ >+template<typename MatrixType> >+template<typename InputType> >+PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix) >+ : m_lu(matrix.derived()), > m_p(matrix.rows()), > m_rowsTranspositions(matrix.rows()), > m_l1_norm(0), > m_det_p(0), > m_isInitialized(false) > { > compute(matrix.derived()); > } >diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h >--- a/Eigen/src/QR/ColPivHouseholderQR.h >+++ b/Eigen/src/QR/ColPivHouseholderQR.h >@@ -46,25 +46,23 @@ template<typename _MatrixType> struct tr > template<typename _MatrixType> class ColPivHouseholderQR > { > public: > > typedef _MatrixType MatrixType; > enum { > RowsAtCompileTime = MatrixType::RowsAtCompileTime, > ColsAtCompileTime = MatrixType::ColsAtCompileTime, >- Options = MatrixType::Options, > MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, > MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime > }; > typedef typename MatrixType::Scalar Scalar; > typedef typename MatrixType::RealScalar RealScalar; > // FIXME should be int > typedef typename MatrixType::StorageIndex StorageIndex; >- typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; > typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; > typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; > typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; > typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; > typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; > typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; > typedef typename MatrixType::PlainObject PlainObject; > >@@ -130,16 +128,37 @@ template<typename _MatrixType> class Col > m_colNormsUpdated(matrix.cols()), > m_colNormsDirect(matrix.cols()), > m_isInitialized(false), > m_usePrescribedThreshold(false) > { > compute(matrix.derived()); > } > >+ /** \brief Constructs a QR factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. >+ * >+ * \sa ColPivHouseholderQR(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) >+ : m_qr(matrix.derived()), >+ m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), >+ m_colsPermutation(PermIndexType(matrix.cols())), >+ m_colsTranspositions(matrix.cols()), >+ m_temp(matrix.cols()), >+ m_colNormsUpdated(matrix.cols()), >+ m_colNormsDirect(matrix.cols()), >+ m_isInitialized(false), >+ m_usePrescribedThreshold(false) >+ { >+ compute(matrix.derived()); >+ } >+ > /** This method finds a solution x to the equation Ax=b, where A is the matrix of which > * *this is the QR decomposition, if any exists. > * > * \param b the right-hand-side of the equation to solve. > * > * \returns a solution. > * > * \note The case where b is a matrix is not yet implemented. Also, this >diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h >--- a/Eigen/src/QR/FullPivHouseholderQR.h >+++ b/Eigen/src/QR/FullPivHouseholderQR.h >@@ -55,17 +55,16 @@ struct traits<FullPivHouseholderQRMatrix > template<typename _MatrixType> class FullPivHouseholderQR > { > public: > > typedef _MatrixType MatrixType; > enum { > RowsAtCompileTime = MatrixType::RowsAtCompileTime, > ColsAtCompileTime = MatrixType::ColsAtCompileTime, >- Options = MatrixType::Options, > MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, > MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime > }; > typedef typename MatrixType::Scalar Scalar; > typedef typename MatrixType::RealScalar RealScalar; > // FIXME should be int > typedef typename MatrixType::StorageIndex StorageIndex; > typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType; >@@ -130,16 +129,36 @@ template<typename _MatrixType> class Ful > m_cols_permutation(matrix.cols()), > m_temp(matrix.cols()), > m_isInitialized(false), > m_usePrescribedThreshold(false) > { > compute(matrix.derived()); > } > >+ /** \brief Constructs a QR factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. >+ * >+ * \sa FullPivHouseholderQR(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit FullPivHouseholderQR(EigenBase<InputType>& matrix) >+ : m_qr(matrix.derived()), >+ m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), >+ m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), >+ m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), >+ m_cols_permutation(matrix.cols()), >+ m_temp(matrix.cols()), >+ m_isInitialized(false), >+ m_usePrescribedThreshold(false) >+ { >+ compute(matrix.derived()); >+ } >+ > /** This method finds a solution x to the equation Ax=b, where A is the matrix of which > * \c *this is the QR decomposition. > * > * \param b the right-hand-side of the equation to solve. > * > * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A, > * and an arbitrary solution otherwise. > * >diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h >--- a/Eigen/src/QR/HouseholderQR.h >+++ b/Eigen/src/QR/HouseholderQR.h >@@ -42,17 +42,16 @@ namespace Eigen { > template<typename _MatrixType> class HouseholderQR > { > public: > > typedef _MatrixType MatrixType; > enum { > RowsAtCompileTime = MatrixType::RowsAtCompileTime, > ColsAtCompileTime = MatrixType::ColsAtCompileTime, >- Options = MatrixType::Options, > MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, > MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime > }; > typedef typename MatrixType::Scalar Scalar; > typedef typename MatrixType::RealScalar RealScalar; > // FIXME should be int > typedef typename MatrixType::StorageIndex StorageIndex; > typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType; >@@ -97,16 +96,34 @@ template<typename _MatrixType> class Hou > : m_qr(matrix.rows(), matrix.cols()), > m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), > m_temp(matrix.cols()), > m_isInitialized(false) > { > compute(matrix.derived()); > } > >+ >+ /** \brief Constructs a QR factorization from a given matrix >+ * >+ * This overloaded constructor is provided for inplace solving when >+ * \c MatrixType is a Eigen::Ref. >+ * >+ * \sa HouseholderQR(const EigenBase&) >+ */ >+ template<typename InputType> >+ explicit HouseholderQR(EigenBase<InputType>& matrix) >+ : m_qr(matrix.derived()), >+ m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), >+ m_temp(matrix.cols()), >+ m_isInitialized(false) >+ { >+ compute(matrix.derived()); >+ } >+ > /** This method finds a solution x to the equation Ax=b, where A is the matrix of which > * *this is the QR decomposition, if any exists. > * > * \param b the right-hand-side of the equation to solve. > * > * \returns a solution. > * > * \note The case where b is a matrix is not yet implemented. Also, this
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Diff
View Attachment As Raw
Actions:
View
|
Diff
Attachments on
bug 707
: 714 |
715