| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H |
| #define EIGEN_SIMPLICIAL_CHOLESKY_H |
| |
| namespace Eigen { |
| |
| enum SimplicialCholeskyMode { |
| SimplicialCholeskyLLT, |
| SimplicialCholeskyLDLT |
| }; |
| |
| namespace internal { |
| template<typename CholMatrixType, typename InputMatrixType> |
| struct simplicial_cholesky_grab_input { |
| typedef CholMatrixType const * ConstCholMatrixPtr; |
| static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp) |
| { |
| tmp = input; |
| pmat = &tmp; |
| } |
| }; |
| |
| template<typename MatrixType> |
| struct simplicial_cholesky_grab_input<MatrixType,MatrixType> { |
| typedef MatrixType const * ConstMatrixPtr; |
| static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/) |
| { |
| pmat = &input; |
| } |
| }; |
| } // end namespace internal |
| |
| /** \ingroup SparseCholesky_Module |
| * \brief A base class for direct sparse Cholesky factorizations |
| * |
| * This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are |
| * selfadjoint and positive definite. These factorizations allow for solving A.X = B where |
| * X and B can be either dense or sparse. |
| * |
| * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| * such that the factorized matrix is P A P^-1. |
| * |
| * \tparam Derived the type of the derived class, that is the actual factorization type. |
| * |
| */ |
| template<typename Derived> |
| class SimplicialCholeskyBase : public SparseSolverBase<Derived> |
| { |
| typedef SparseSolverBase<Derived> Base; |
| using Base::m_isInitialized; |
| |
| public: |
| typedef typename internal::traits<Derived>::MatrixType MatrixType; |
| typedef typename internal::traits<Derived>::OrderingType OrderingType; |
| enum { UpLo = internal::traits<Derived>::UpLo }; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; |
| typedef CholMatrixType const * ConstCholMatrixPtr; |
| typedef Matrix<Scalar,Dynamic,1> VectorType; |
| typedef Matrix<StorageIndex,Dynamic,1> VectorI; |
| |
| enum { |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
| }; |
| |
| public: |
| |
| using Base::derived; |
| |
| /** Default constructor */ |
| SimplicialCholeskyBase() |
| : m_info(Success), |
| m_factorizationIsOk(false), |
| m_analysisIsOk(false), |
| m_shiftOffset(0), |
| m_shiftScale(1) |
| {} |
| |
| explicit SimplicialCholeskyBase(const MatrixType& matrix) |
| : m_info(Success), |
| m_factorizationIsOk(false), |
| m_analysisIsOk(false), |
| m_shiftOffset(0), |
| m_shiftScale(1) |
| { |
| derived().compute(matrix); |
| } |
| |
| ~SimplicialCholeskyBase() |
| { |
| } |
| |
| Derived& derived() { return *static_cast<Derived*>(this); } |
| const Derived& derived() const { return *static_cast<const Derived*>(this); } |
| |
| inline Index cols() const { return m_matrix.cols(); } |
| inline Index rows() const { return m_matrix.rows(); } |
| |
| /** \brief Reports whether previous computation was successful. |
| * |
| * \returns \c Success if computation was successful, |
| * \c NumericalIssue if the matrix.appears to be negative. |
| */ |
| ComputationInfo info() const |
| { |
| eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| return m_info; |
| } |
| |
| /** \returns the permutation P |
| * \sa permutationPinv() */ |
| const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const |
| { return m_P; } |
| |
| /** \returns the inverse P^-1 of the permutation P |
| * \sa permutationP() */ |
| const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const |
| { return m_Pinv; } |
| |
| /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization. |
| * |
| * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n |
| * \c d_ii = \a offset + \a scale * \c d_ii |
| * |
| * The default is the identity transformation with \a offset=0, and \a scale=1. |
| * |
| * \returns a reference to \c *this. |
| */ |
| Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) |
| { |
| m_shiftOffset = offset; |
| m_shiftScale = scale; |
| return derived(); |
| } |
| |
| #ifndef EIGEN_PARSED_BY_DOXYGEN |
| /** \internal */ |
| template<typename Stream> |
| void dumpMemory(Stream& s) |
| { |
| int total = 0; |
| s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; |
| s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; |
| s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; |
| s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; |
| } |
| |
| /** \internal */ |
| template<typename Rhs,typename Dest> |
| void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const |
| { |
| eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); |
| eigen_assert(m_matrix.rows()==b.rows()); |
| |
| if(m_info!=Success) |
| return; |
| |
| if(m_P.size()>0) |
| dest = m_P * b; |
| else |
| dest = b; |
| |
| if(m_matrix.nonZeros()>0) // otherwise L==I |
| derived().matrixL().solveInPlace(dest); |
| |
| if(m_diag.size()>0) |
| dest = m_diag.asDiagonal().inverse() * dest; |
| |
| if (m_matrix.nonZeros()>0) // otherwise U==I |
| derived().matrixU().solveInPlace(dest); |
| |
| if(m_P.size()>0) |
| dest = m_Pinv * dest; |
| } |
| |
| template<typename Rhs,typename Dest> |
| void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const |
| { |
| internal::solve_sparse_through_dense_panels(derived(), b, dest); |
| } |
| |
| #endif // EIGEN_PARSED_BY_DOXYGEN |
| |
| protected: |
| |
| /** Computes the sparse Cholesky decomposition of \a matrix */ |
| template<bool DoLDLT> |
| void compute(const MatrixType& matrix) |
| { |
| eigen_assert(matrix.rows()==matrix.cols()); |
| Index size = matrix.cols(); |
| CholMatrixType tmp(size,size); |
| ConstCholMatrixPtr pmat; |
| ordering(matrix, pmat, tmp); |
| analyzePattern_preordered(*pmat, DoLDLT); |
| factorize_preordered<DoLDLT>(*pmat); |
| } |
| |
| template<bool DoLDLT> |
| void factorize(const MatrixType& a) |
| { |
| eigen_assert(a.rows()==a.cols()); |
| Index size = a.cols(); |
| CholMatrixType tmp(size,size); |
| ConstCholMatrixPtr pmat; |
| |
| if(m_P.size()==0 && (UpLo&Upper)==Upper) |
| { |
| // If there is no ordering, try to directly use the input matrix without any copy |
| internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp); |
| } |
| else |
| { |
| tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); |
| pmat = &tmp; |
| } |
| |
| factorize_preordered<DoLDLT>(*pmat); |
| } |
| |
| template<bool DoLDLT> |
| void factorize_preordered(const CholMatrixType& a); |
| |
| void analyzePattern(const MatrixType& a, bool doLDLT) |
| { |
| eigen_assert(a.rows()==a.cols()); |
| Index size = a.cols(); |
| CholMatrixType tmp(size,size); |
| ConstCholMatrixPtr pmat; |
| ordering(a, pmat, tmp); |
| analyzePattern_preordered(*pmat,doLDLT); |
| } |
| void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); |
| |
| void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap); |
| |
| /** keeps off-diagonal entries; drops diagonal entries */ |
| struct keep_diag { |
| inline bool operator() (const Index& row, const Index& col, const Scalar&) const |
| { |
| return row!=col; |
| } |
| }; |
| |
| mutable ComputationInfo m_info; |
| bool m_factorizationIsOk; |
| bool m_analysisIsOk; |
| |
| CholMatrixType m_matrix; |
| VectorType m_diag; // the diagonal coefficients (LDLT mode) |
| VectorI m_parent; // elimination tree |
| VectorI m_nonZerosPerCol; |
| PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation |
| PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation |
| |
| RealScalar m_shiftOffset; |
| RealScalar m_shiftScale; |
| }; |
| |
| template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT; |
| template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT; |
| template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky; |
| |
| namespace internal { |
| |
| template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > |
| { |
| typedef _MatrixType MatrixType; |
| typedef _Ordering OrderingType; |
| enum { UpLo = _UpLo }; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; |
| typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL; |
| typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; |
| static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } |
| static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } |
| }; |
| |
| template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > |
| { |
| typedef _MatrixType MatrixType; |
| typedef _Ordering OrderingType; |
| enum { UpLo = _UpLo }; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; |
| typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL; |
| typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; |
| static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } |
| static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } |
| }; |
| |
| template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > |
| { |
| typedef _MatrixType MatrixType; |
| typedef _Ordering OrderingType; |
| enum { UpLo = _UpLo }; |
| }; |
| |
| } |
| |
| /** \ingroup SparseCholesky_Module |
| * \class SimplicialLLT |
| * \brief A direct sparse LLT Cholesky factorizations |
| * |
| * This class provides a LL^T Cholesky factorizations of sparse matrices that are |
| * selfadjoint and positive definite. The factorization allows for solving A.X = B where |
| * X and B can be either dense or sparse. |
| * |
| * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| * such that the factorized matrix is P A P^-1. |
| * |
| * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower |
| * or Upper. Default is Lower. |
| * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> |
| * |
| * \implsparsesolverconcept |
| * |
| * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering |
| */ |
| template<typename _MatrixType, int _UpLo, typename _Ordering> |
| class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > |
| { |
| public: |
| typedef _MatrixType MatrixType; |
| enum { UpLo = _UpLo }; |
| typedef SimplicialCholeskyBase<SimplicialLLT> Base; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; |
| typedef Matrix<Scalar,Dynamic,1> VectorType; |
| typedef internal::traits<SimplicialLLT> Traits; |
| typedef typename Traits::MatrixL MatrixL; |
| typedef typename Traits::MatrixU MatrixU; |
| public: |
| /** Default constructor */ |
| SimplicialLLT() : Base() {} |
| /** Constructs and performs the LLT factorization of \a matrix */ |
| explicit SimplicialLLT(const MatrixType& matrix) |
| : Base(matrix) {} |
| |
| /** \returns an expression of the factor L */ |
| inline const MatrixL matrixL() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); |
| return Traits::getL(Base::m_matrix); |
| } |
| |
| /** \returns an expression of the factor U (= L^*) */ |
| inline const MatrixU matrixU() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); |
| return Traits::getU(Base::m_matrix); |
| } |
| |
| /** Computes the sparse Cholesky decomposition of \a matrix */ |
| SimplicialLLT& compute(const MatrixType& matrix) |
| { |
| Base::template compute<false>(matrix); |
| return *this; |
| } |
| |
| /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| * |
| * This function is particularly useful when solving for several problems having the same structure. |
| * |
| * \sa factorize() |
| */ |
| void analyzePattern(const MatrixType& a) |
| { |
| Base::analyzePattern(a, false); |
| } |
| |
| /** Performs a numeric decomposition of \a matrix |
| * |
| * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| * |
| * \sa analyzePattern() |
| */ |
| void factorize(const MatrixType& a) |
| { |
| Base::template factorize<false>(a); |
| } |
| |
| /** \returns the determinant of the underlying matrix from the current factorization */ |
| Scalar determinant() const |
| { |
| Scalar detL = Base::m_matrix.diagonal().prod(); |
| return numext::abs2(detL); |
| } |
| }; |
| |
| /** \ingroup SparseCholesky_Module |
| * \class SimplicialLDLT |
| * \brief A direct sparse LDLT Cholesky factorizations without square root. |
| * |
| * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are |
| * selfadjoint and positive definite. The factorization allows for solving A.X = B where |
| * X and B can be either dense or sparse. |
| * |
| * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization |
| * such that the factorized matrix is P A P^-1. |
| * |
| * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> |
| * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower |
| * or Upper. Default is Lower. |
| * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> |
| * |
| * \implsparsesolverconcept |
| * |
| * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering |
| */ |
| template<typename _MatrixType, int _UpLo, typename _Ordering> |
| class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > |
| { |
| public: |
| typedef _MatrixType MatrixType; |
| enum { UpLo = _UpLo }; |
| typedef SimplicialCholeskyBase<SimplicialLDLT> Base; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; |
| typedef Matrix<Scalar,Dynamic,1> VectorType; |
| typedef internal::traits<SimplicialLDLT> Traits; |
| typedef typename Traits::MatrixL MatrixL; |
| typedef typename Traits::MatrixU MatrixU; |
| public: |
| /** Default constructor */ |
| SimplicialLDLT() : Base() {} |
| |
| /** Constructs and performs the LLT factorization of \a matrix */ |
| explicit SimplicialLDLT(const MatrixType& matrix) |
| : Base(matrix) {} |
| |
| /** \returns a vector expression of the diagonal D */ |
| inline const VectorType vectorD() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| return Base::m_diag; |
| } |
| /** \returns an expression of the factor L */ |
| inline const MatrixL matrixL() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| return Traits::getL(Base::m_matrix); |
| } |
| |
| /** \returns an expression of the factor U (= L^*) */ |
| inline const MatrixU matrixU() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); |
| return Traits::getU(Base::m_matrix); |
| } |
| |
| /** Computes the sparse Cholesky decomposition of \a matrix */ |
| SimplicialLDLT& compute(const MatrixType& matrix) |
| { |
| Base::template compute<true>(matrix); |
| return *this; |
| } |
| |
| /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| * |
| * This function is particularly useful when solving for several problems having the same structure. |
| * |
| * \sa factorize() |
| */ |
| void analyzePattern(const MatrixType& a) |
| { |
| Base::analyzePattern(a, true); |
| } |
| |
| /** Performs a numeric decomposition of \a matrix |
| * |
| * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| * |
| * \sa analyzePattern() |
| */ |
| void factorize(const MatrixType& a) |
| { |
| Base::template factorize<true>(a); |
| } |
| |
| /** \returns the determinant of the underlying matrix from the current factorization */ |
| Scalar determinant() const |
| { |
| return Base::m_diag.prod(); |
| } |
| }; |
| |
| /** \deprecated use SimplicialLDLT or class SimplicialLLT |
| * \ingroup SparseCholesky_Module |
| * \class SimplicialCholesky |
| * |
| * \sa class SimplicialLDLT, class SimplicialLLT |
| */ |
| template<typename _MatrixType, int _UpLo, typename _Ordering> |
| class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > |
| { |
| public: |
| typedef _MatrixType MatrixType; |
| enum { UpLo = _UpLo }; |
| typedef SimplicialCholeskyBase<SimplicialCholesky> Base; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; |
| typedef Matrix<Scalar,Dynamic,1> VectorType; |
| typedef internal::traits<SimplicialCholesky> Traits; |
| typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; |
| typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; |
| public: |
| SimplicialCholesky() : Base(), m_LDLT(true) {} |
| |
| explicit SimplicialCholesky(const MatrixType& matrix) |
| : Base(), m_LDLT(true) |
| { |
| compute(matrix); |
| } |
| |
| SimplicialCholesky& setMode(SimplicialCholeskyMode mode) |
| { |
| switch(mode) |
| { |
| case SimplicialCholeskyLLT: |
| m_LDLT = false; |
| break; |
| case SimplicialCholeskyLDLT: |
| m_LDLT = true; |
| break; |
| default: |
| break; |
| } |
| |
| return *this; |
| } |
| |
| inline const VectorType vectorD() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); |
| return Base::m_diag; |
| } |
| inline const CholMatrixType rawMatrix() const { |
| eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); |
| return Base::m_matrix; |
| } |
| |
| /** Computes the sparse Cholesky decomposition of \a matrix */ |
| SimplicialCholesky& compute(const MatrixType& matrix) |
| { |
| if(m_LDLT) |
| Base::template compute<true>(matrix); |
| else |
| Base::template compute<false>(matrix); |
| return *this; |
| } |
| |
| /** Performs a symbolic decomposition on the sparcity of \a matrix. |
| * |
| * This function is particularly useful when solving for several problems having the same structure. |
| * |
| * \sa factorize() |
| */ |
| void analyzePattern(const MatrixType& a) |
| { |
| Base::analyzePattern(a, m_LDLT); |
| } |
| |
| /** Performs a numeric decomposition of \a matrix |
| * |
| * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. |
| * |
| * \sa analyzePattern() |
| */ |
| void factorize(const MatrixType& a) |
| { |
| if(m_LDLT) |
| Base::template factorize<true>(a); |
| else |
| Base::template factorize<false>(a); |
| } |
| |
| /** \internal */ |
| template<typename Rhs,typename Dest> |
| void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const |
| { |
| eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); |
| eigen_assert(Base::m_matrix.rows()==b.rows()); |
| |
| if(Base::m_info!=Success) |
| return; |
| |
| if(Base::m_P.size()>0) |
| dest = Base::m_P * b; |
| else |
| dest = b; |
| |
| if(Base::m_matrix.nonZeros()>0) // otherwise L==I |
| { |
| if(m_LDLT) |
| LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); |
| else |
| LLTTraits::getL(Base::m_matrix).solveInPlace(dest); |
| } |
| |
| if(Base::m_diag.size()>0) |
| dest = Base::m_diag.asDiagonal().inverse() * dest; |
| |
| if (Base::m_matrix.nonZeros()>0) // otherwise I==I |
| { |
| if(m_LDLT) |
| LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); |
| else |
| LLTTraits::getU(Base::m_matrix).solveInPlace(dest); |
| } |
| |
| if(Base::m_P.size()>0) |
| dest = Base::m_Pinv * dest; |
| } |
| |
| /** \internal */ |
| template<typename Rhs,typename Dest> |
| void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const |
| { |
| internal::solve_sparse_through_dense_panels(*this, b, dest); |
| } |
| |
| Scalar determinant() const |
| { |
| if(m_LDLT) |
| { |
| return Base::m_diag.prod(); |
| } |
| else |
| { |
| Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); |
| return numext::abs2(detL); |
| } |
| } |
| |
| protected: |
| bool m_LDLT; |
| }; |
| |
| template<typename Derived> |
| void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap) |
| { |
| eigen_assert(a.rows()==a.cols()); |
| const Index size = a.rows(); |
| pmat = ≈ |
| // Note that ordering methods compute the inverse permutation |
| if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) |
| { |
| { |
| CholMatrixType C; |
| C = a.template selfadjointView<UpLo>(); |
| |
| OrderingType ordering; |
| ordering(C,m_Pinv); |
| } |
| |
| if(m_Pinv.size()>0) m_P = m_Pinv.inverse(); |
| else m_P.resize(0); |
| |
| ap.resize(size,size); |
| ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); |
| } |
| else |
| { |
| m_Pinv.resize(0); |
| m_P.resize(0); |
| if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor) |
| { |
| // we have to transpose the lower part to to the upper one |
| ap.resize(size,size); |
| ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>(); |
| } |
| else |
| internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap); |
| } |
| } |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_SIMPLICIAL_CHOLESKY_H |
