src/SparseCholesky/SimplicialCholesky.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
0011 #define EIGEN_SIMPLICIAL_CHOLESKY_H
0012
0013 namespace RivetEigen {
0014
0015 enum SimplicialCholeskyMode {
0016   SimplicialCholeskyLLT,
0017   SimplicialCholeskyLDLT
0018 };
0019
0020 namespace internal {
0021   template<typename CholMatrixType, typename InputMatrixType>
0022   struct simplicial_cholesky_grab_input {
0023     typedef CholMatrixType const * ConstCholMatrixPtr;
0024     static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
0025     {
0026       tmp = input;
0027       pmat = &tmp;
0028     }
0029   };
0030
0031   template<typename MatrixType>
0032   struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
0033     typedef MatrixType const * ConstMatrixPtr;
0034     static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
0035     {
0036       pmat = &input;
0037     }
0038   };
0039 } // end namespace internal
0040
0041 /** \ingroup SparseCholesky_Module
0042   * \brief A base class for direct sparse Cholesky factorizations
0043   *
0044   * This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are
0045   * selfadjoint and positive definite. These factorizations allow for solving A.X = B where
0046   * X and B can be either dense or sparse.
0047   *
0048   * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
0049   * such that the factorized matrix is P A P^-1.
0050   *
0051   * \tparam Derived the type of the derived class, that is the actual factorization type.
0052   *
0053   */
0054 template<typename Derived>
0055 class SimplicialCholeskyBase : public SparseSolverBase<Derived>
0056 {
0057     typedef SparseSolverBase<Derived> Base;
0058     using Base::m_isInitialized;
0059
0060   public:
0061     typedef typename internal::traits<Derived>::MatrixType MatrixType;
0062     typedef typename internal::traits<Derived>::OrderingType OrderingType;
0063     enum { UpLo = internal::traits<Derived>::UpLo };
0064     typedef typename MatrixType::Scalar Scalar;
0065     typedef typename MatrixType::RealScalar RealScalar;
0066     typedef typename MatrixType::StorageIndex StorageIndex;
0067     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
0068     typedef CholMatrixType const * ConstCholMatrixPtr;
0069     typedef Matrix<Scalar,Dynamic,1> VectorType;
0070     typedef Matrix<StorageIndex,Dynamic,1> VectorI;
0071
0072     enum {
0073       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0074       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0075     };
0076
0077   public:
0078
0079     using Base::derived;
0080
0081     /** Default constructor */
0082     SimplicialCholeskyBase()
0083       : m_info(Success),
0084         m_factorizationIsOk(false),
0085         m_analysisIsOk(false),
0086         m_shiftOffset(0),
0087         m_shiftScale(1)
0088     {}
0089
0090     explicit SimplicialCholeskyBase(const MatrixType& matrix)
0091       : m_info(Success),
0092         m_factorizationIsOk(false),
0093         m_analysisIsOk(false),
0094         m_shiftOffset(0),
0095         m_shiftScale(1)
0096     {
0097       derived().compute(matrix);
0098     }
0099
0100     ~SimplicialCholeskyBase()
0101     {
0102     }
0103
0104     Derived& derived() { return *static_cast<Derived*>(this); }
0105     const Derived& derived() const { return *static_cast<const Derived*>(this); }
0106
0107     inline Index cols() const { return m_matrix.cols(); }
0108     inline Index rows() const { return m_matrix.rows(); }
0109
0110     /** \brief Reports whether previous computation was successful.
0111       *
0112       * \returns \c Success if computation was successful,
0113       *          \c NumericalIssue if the matrix.appears to be negative.
0114       */
0115     ComputationInfo info() const
0116     {
0117       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0118       return m_info;
0119     }
0120
0121     /** \returns the permutation P
0122       * \sa permutationPinv() */
0123     const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const
0124     { return m_P; }
0125
0126     /** \returns the inverse P^-1 of the permutation P
0127       * \sa permutationP() */
0128     const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const
0129     { return m_Pinv; }
0130
0131     /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
0132       *
0133       * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n
0134       * \c d_ii = \a offset + \a scale * \c d_ii
0135       *
0136       * The default is the identity transformation with \a offset=0, and \a scale=1.
0137       *
0138       * \returns a reference to \c *this.
0139       */
0140     Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
0141     {
0142       m_shiftOffset = offset;
0143       m_shiftScale = scale;
0144       return derived();
0145     }
0146
0147 #ifndef EIGEN_PARSED_BY_DOXYGEN
0148     /** \internal */
0149     template<typename Stream>
0150     void dumpMemory(Stream& s)
0151     {
0152       int total = 0;
0153       s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
0154       s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
0155       s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
0156       s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
0157       s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
0158       s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
0159       s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
0160     }
0161
0162     /** \internal */
0163     template<typename Rhs,typename Dest>
0164     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
0165     {
0166       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
0167       eigen_assert(m_matrix.rows()==b.rows());
0168
0169       if(m_info!=Success)
0170         return;
0171
0172       if(m_P.size()>0)
0173         dest = m_P * b;
0174       else
0175         dest = b;
0176
0177       if(m_matrix.nonZeros()>0) // otherwise L==I
0178         derived().matrixL().solveInPlace(dest);
0179
0180       if(m_diag.size()>0)
0181         dest = m_diag.asDiagonal().inverse() * dest;
0182
0183       if (m_matrix.nonZeros()>0) // otherwise U==I
0184         derived().matrixU().solveInPlace(dest);
0185
0186       if(m_P.size()>0)
0187         dest = m_Pinv * dest;
0188     }
0189
0190     template<typename Rhs,typename Dest>
0191     void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
0192     {
0193       internal::solve_sparse_through_dense_panels(derived(), b, dest);
0194     }
0195
0196 #endif // EIGEN_PARSED_BY_DOXYGEN
0197
0198   protected:
0199
0200     /** Computes the sparse Cholesky decomposition of \a matrix */
0201     template<bool DoLDLT>
0202     void compute(const MatrixType& matrix)
0203     {
0204       eigen_assert(matrix.rows()==matrix.cols());
0205       Index size = matrix.cols();
0206       CholMatrixType tmp(size,size);
0207       ConstCholMatrixPtr pmat;
0208       ordering(matrix, pmat, tmp);
0209       analyzePattern_preordered(*pmat, DoLDLT);
0210       factorize_preordered<DoLDLT>(*pmat);
0211     }
0212
0213     template<bool DoLDLT>
0214     void factorize(const MatrixType& a)
0215     {
0216       eigen_assert(a.rows()==a.cols());
0217       Index size = a.cols();
0218       CholMatrixType tmp(size,size);
0219       ConstCholMatrixPtr pmat;
0220
0221       if(m_P.size() == 0 && (int(UpLo) & int(Upper)) == Upper)
0222       {
0223         // If there is no ordering, try to directly use the input matrix without any copy
0224         internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
0225       }
0226       else
0227       {
0228         tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
0229         pmat = &tmp;
0230       }
0231
0232       factorize_preordered<DoLDLT>(*pmat);
0233     }
0234
0235     template<bool DoLDLT>
0236     void factorize_preordered(const CholMatrixType& a);
0237
0238     void analyzePattern(const MatrixType& a, bool doLDLT)
0239     {
0240       eigen_assert(a.rows()==a.cols());
0241       Index size = a.cols();
0242       CholMatrixType tmp(size,size);
0243       ConstCholMatrixPtr pmat;
0244       ordering(a, pmat, tmp);
0245       analyzePattern_preordered(*pmat,doLDLT);
0246     }
0247     void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
0248
0249     void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
0250
0251     /** keeps off-diagonal entries; drops diagonal entries */
0252     struct keep_diag {
0253       inline bool operator() (const Index& row, const Index& col, const Scalar&) const
0254       {
0255         return row!=col;
0256       }
0257     };
0258
0259     mutable ComputationInfo m_info;
0260     bool m_factorizationIsOk;
0261     bool m_analysisIsOk;
0262
0263     CholMatrixType m_matrix;
0264     VectorType m_diag;                                // the diagonal coefficients (LDLT mode)
0265     VectorI m_parent;                                 // elimination tree
0266     VectorI m_nonZerosPerCol;
0267     PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P;     // the permutation
0268     PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv;  // the inverse permutation
0269
0270     RealScalar m_shiftOffset;
0271     RealScalar m_shiftScale;
0272 };
0273
0274 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
0275 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
0276 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
0277
0278 namespace internal {
0279
0280 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
0281 {
0282   typedef _MatrixType MatrixType;
0283   typedef _Ordering OrderingType;
0284   enum { UpLo = _UpLo };
0285   typedef typename MatrixType::Scalar                         Scalar;
0286   typedef typename MatrixType::StorageIndex                   StorageIndex;
0287   typedef SparseMatrix<Scalar, ColMajor, StorageIndex>        CholMatrixType;
0288   typedef TriangularView<const CholMatrixType, RivetEigen::Lower>  MatrixL;
0289   typedef TriangularView<const typename CholMatrixType::AdjointReturnType, RivetEigen::Upper>   MatrixU;
0290   static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
0291   static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
0292 };
0293
0294 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
0295 {
0296   typedef _MatrixType MatrixType;
0297   typedef _Ordering OrderingType;
0298   enum { UpLo = _UpLo };
0299   typedef typename MatrixType::Scalar                             Scalar;
0300   typedef typename MatrixType::StorageIndex                       StorageIndex;
0301   typedef SparseMatrix<Scalar, ColMajor, StorageIndex>            CholMatrixType;
0302   typedef TriangularView<const CholMatrixType, RivetEigen::UnitLower>  MatrixL;
0303   typedef TriangularView<const typename CholMatrixType::AdjointReturnType, RivetEigen::UnitUpper> MatrixU;
0304   static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
0305   static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
0306 };
0307
0308 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
0309 {
0310   typedef _MatrixType MatrixType;
0311   typedef _Ordering OrderingType;
0312   enum { UpLo = _UpLo };
0313 };
0314
0315 }
0316
0317 /** \ingroup SparseCholesky_Module
0318   * \class SimplicialLLT
0319   * \brief A direct sparse LLT Cholesky factorizations
0320   *
0321   * This class provides a LL^T Cholesky factorizations of sparse matrices that are
0322   * selfadjoint and positive definite. The factorization allows for solving A.X = B where
0323   * X and B can be either dense or sparse.
0324   *
0325   * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
0326   * such that the factorized matrix is P A P^-1.
0327   *
0328   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0329   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
0330   *               or Upper. Default is Lower.
0331   * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
0332   *
0333   * \implsparsesolverconcept
0334   *
0335   * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
0336   */
0337 template<typename _MatrixType, int _UpLo, typename _Ordering>
0338     class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
0339 {
0340 public:
0341     typedef _MatrixType MatrixType;
0342     enum { UpLo = _UpLo };
0343     typedef SimplicialCholeskyBase<SimplicialLLT> Base;
0344     typedef typename MatrixType::Scalar Scalar;
0345     typedef typename MatrixType::RealScalar RealScalar;
0346     typedef typename MatrixType::StorageIndex StorageIndex;
0347     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
0348     typedef Matrix<Scalar,Dynamic,1> VectorType;
0349     typedef internal::traits<SimplicialLLT> Traits;
0350     typedef typename Traits::MatrixL  MatrixL;
0351     typedef typename Traits::MatrixU  MatrixU;
0352 public:
0353     /** Default constructor */
0354     SimplicialLLT() : Base() {}
0355     /** Constructs and performs the LLT factorization of \a matrix */
0356     explicit SimplicialLLT(const MatrixType& matrix)
0357         : Base(matrix) {}
0358
0359     /** \returns an expression of the factor L */
0360     inline const MatrixL matrixL() const {
0361         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
0362         return Traits::getL(Base::m_matrix);
0363     }
0364
0365     /** \returns an expression of the factor U (= L^*) */
0366     inline const MatrixU matrixU() const {
0367         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
0368         return Traits::getU(Base::m_matrix);
0369     }
0370
0371     /** Computes the sparse Cholesky decomposition of \a matrix */
0372     SimplicialLLT& compute(const MatrixType& matrix)
0373     {
0374       Base::template compute<false>(matrix);
0375       return *this;
0376     }
0377
0378     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0379       *
0380       * This function is particularly useful when solving for several problems having the same structure.
0381       *
0382       * \sa factorize()
0383       */
0384     void analyzePattern(const MatrixType& a)
0385     {
0386       Base::analyzePattern(a, false);
0387     }
0388
0389     /** Performs a numeric decomposition of \a matrix
0390       *
0391       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
0392       *
0393       * \sa analyzePattern()
0394       */
0395     void factorize(const MatrixType& a)
0396     {
0397       Base::template factorize<false>(a);
0398     }
0399
0400     /** \returns the determinant of the underlying matrix from the current factorization */
0401     Scalar determinant() const
0402     {
0403       Scalar detL = Base::m_matrix.diagonal().prod();
0404       return numext::abs2(detL);
0405     }
0406 };
0407
0408 /** \ingroup SparseCholesky_Module
0409   * \class SimplicialLDLT
0410   * \brief A direct sparse LDLT Cholesky factorizations without square root.
0411   *
0412   * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
0413   * selfadjoint and positive definite. The factorization allows for solving A.X = B where
0414   * X and B can be either dense or sparse.
0415   *
0416   * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
0417   * such that the factorized matrix is P A P^-1.
0418   *
0419   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0420   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
0421   *               or Upper. Default is Lower.
0422   * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
0423   *
0424   * \implsparsesolverconcept
0425   *
0426   * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
0427   */
0428 template<typename _MatrixType, int _UpLo, typename _Ordering>
0429     class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
0430 {
0431 public:
0432     typedef _MatrixType MatrixType;
0433     enum { UpLo = _UpLo };
0434     typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
0435     typedef typename MatrixType::Scalar Scalar;
0436     typedef typename MatrixType::RealScalar RealScalar;
0437     typedef typename MatrixType::StorageIndex StorageIndex;
0438     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
0439     typedef Matrix<Scalar,Dynamic,1> VectorType;
0440     typedef internal::traits<SimplicialLDLT> Traits;
0441     typedef typename Traits::MatrixL  MatrixL;
0442     typedef typename Traits::MatrixU  MatrixU;
0443 public:
0444     /** Default constructor */
0445     SimplicialLDLT() : Base() {}
0446
0447     /** Constructs and performs the LLT factorization of \a matrix */
0448     explicit SimplicialLDLT(const MatrixType& matrix)
0449         : Base(matrix) {}
0450
0451     /** \returns a vector expression of the diagonal D */
0452     inline const VectorType vectorD() const {
0453         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
0454         return Base::m_diag;
0455     }
0456     /** \returns an expression of the factor L */
0457     inline const MatrixL matrixL() const {
0458         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
0459         return Traits::getL(Base::m_matrix);
0460     }
0461
0462     /** \returns an expression of the factor U (= L^*) */
0463     inline const MatrixU matrixU() const {
0464         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
0465         return Traits::getU(Base::m_matrix);
0466     }
0467
0468     /** Computes the sparse Cholesky decomposition of \a matrix */
0469     SimplicialLDLT& compute(const MatrixType& matrix)
0470     {
0471       Base::template compute<true>(matrix);
0472       return *this;
0473     }
0474
0475     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0476       *
0477       * This function is particularly useful when solving for several problems having the same structure.
0478       *
0479       * \sa factorize()
0480       */
0481     void analyzePattern(const MatrixType& a)
0482     {
0483       Base::analyzePattern(a, true);
0484     }
0485
0486     /** Performs a numeric decomposition of \a matrix
0487       *
0488       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
0489       *
0490       * \sa analyzePattern()
0491       */
0492     void factorize(const MatrixType& a)
0493     {
0494       Base::template factorize<true>(a);
0495     }
0496
0497     /** \returns the determinant of the underlying matrix from the current factorization */
0498     Scalar determinant() const
0499     {
0500       return Base::m_diag.prod();
0501     }
0502 };
0503
0504 /** \deprecated use SimplicialLDLT or class SimplicialLLT
0505   * \ingroup SparseCholesky_Module
0506   * \class SimplicialCholesky
0507   *
0508   * \sa class SimplicialLDLT, class SimplicialLLT
0509   */
0510 template<typename _MatrixType, int _UpLo, typename _Ordering>
0511     class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
0512 {
0513 public:
0514     typedef _MatrixType MatrixType;
0515     enum { UpLo = _UpLo };
0516     typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
0517     typedef typename MatrixType::Scalar Scalar;
0518     typedef typename MatrixType::RealScalar RealScalar;
0519     typedef typename MatrixType::StorageIndex StorageIndex;
0520     typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
0521     typedef Matrix<Scalar,Dynamic,1> VectorType;
0522     typedef internal::traits<SimplicialCholesky> Traits;
0523     typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
0524     typedef internal::traits<SimplicialLLT<MatrixType,UpLo>  > LLTTraits;
0525   public:
0526     SimplicialCholesky() : Base(), m_LDLT(true) {}
0527
0528     explicit SimplicialCholesky(const MatrixType& matrix)
0529       : Base(), m_LDLT(true)
0530     {
0531       compute(matrix);
0532     }
0533
0534     SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
0535     {
0536       switch(mode)
0537       {
0538       case SimplicialCholeskyLLT:
0539         m_LDLT = false;
0540         break;
0541       case SimplicialCholeskyLDLT:
0542         m_LDLT = true;
0543         break;
0544       default:
0545         break;
0546       }
0547
0548       return *this;
0549     }
0550
0551     inline const VectorType vectorD() const {
0552         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
0553         return Base::m_diag;
0554     }
0555     inline const CholMatrixType rawMatrix() const {
0556         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
0557         return Base::m_matrix;
0558     }
0559
0560     /** Computes the sparse Cholesky decomposition of \a matrix */
0561     SimplicialCholesky& compute(const MatrixType& matrix)
0562     {
0563       if(m_LDLT)
0564         Base::template compute<true>(matrix);
0565       else
0566         Base::template compute<false>(matrix);
0567       return *this;
0568     }
0569
0570     /** Performs a symbolic decomposition on the sparcity of \a matrix.
0571       *
0572       * This function is particularly useful when solving for several problems having the same structure.
0573       *
0574       * \sa factorize()
0575       */
0576     void analyzePattern(const MatrixType& a)
0577     {
0578       Base::analyzePattern(a, m_LDLT);
0579     }
0580
0581     /** Performs a numeric decomposition of \a matrix
0582       *
0583       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
0584       *
0585       * \sa analyzePattern()
0586       */
0587     void factorize(const MatrixType& a)
0588     {
0589       if(m_LDLT)
0590         Base::template factorize<true>(a);
0591       else
0592         Base::template factorize<false>(a);
0593     }
0594
0595     /** \internal */
0596     template<typename Rhs,typename Dest>
0597     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
0598     {
0599       eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
0600       eigen_assert(Base::m_matrix.rows()==b.rows());
0601
0602       if(Base::m_info!=Success)
0603         return;
0604
0605       if(Base::m_P.size()>0)
0606         dest = Base::m_P * b;
0607       else
0608         dest = b;
0609
0610       if(Base::m_matrix.nonZeros()>0) // otherwise L==I
0611       {
0612         if(m_LDLT)
0613           LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
0614         else
0615           LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
0616       }
0617
0618       if(Base::m_diag.size()>0)
0619         dest = Base::m_diag.real().asDiagonal().inverse() * dest;
0620
0621       if (Base::m_matrix.nonZeros()>0) // otherwise I==I
0622       {
0623         if(m_LDLT)
0624           LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
0625         else
0626           LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
0627       }
0628
0629       if(Base::m_P.size()>0)
0630         dest = Base::m_Pinv * dest;
0631     }
0632
0633     /** \internal */
0634     template<typename Rhs,typename Dest>
0635     void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
0636     {
0637       internal::solve_sparse_through_dense_panels(*this, b, dest);
0638     }
0639
0640     Scalar determinant() const
0641     {
0642       if(m_LDLT)
0643       {
0644         return Base::m_diag.prod();
0645       }
0646       else
0647       {
0648         Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
0649         return numext::abs2(detL);
0650       }
0651     }
0652
0653   protected:
0654     bool m_LDLT;
0655 };
0656
0657 template<typename Derived>
0658 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
0659 {
0660   eigen_assert(a.rows()==a.cols());
0661   const Index size = a.rows();
0662   pmat = &ap;
0663   // Note that ordering methods compute the inverse permutation
0664   if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
0665   {
0666     {
0667       CholMatrixType C;
0668       C = a.template selfadjointView<UpLo>();
0669
0670       OrderingType ordering;
0671       ordering(C,m_Pinv);
0672     }
0673
0674     if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
0675     else                m_P.resize(0);
0676
0677     ap.resize(size,size);
0678     ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
0679   }
0680   else
0681   {
0682     m_Pinv.resize(0);
0683     m_P.resize(0);
0684     if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
0685     {
0686       // we have to transpose the lower part to to the upper one
0687       ap.resize(size,size);
0688       ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
0689     }
0690     else
0691       internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
0692   }
0693 }
0694
0695 } // end namespace RivetEigen
0696
0697 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H