diff options
Diffstat (limited to 'Eigen/src/SparseCholesky/SimplicialCholesky.h')
| -rw-r--r-- | Eigen/src/SparseCholesky/SimplicialCholesky.h | 697 |
1 files changed, 697 insertions, 0 deletions
diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h new file mode 100644 index 0000000..9f93e32 --- /dev/null +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -0,0 +1,697 @@ +// 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 && (int(UpLo) & int(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 CholMatrixType& m) { return MatrixL(m); } + static inline MatrixU getU(const CholMatrixType& 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 CholMatrixType& m) { return MatrixL(m); } + static inline MatrixU getU(const CholMatrixType& 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.real().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 |