1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
11#define EIGEN_SIMPLICIAL_CHOLESKY_H
12
13namespace Eigen {
14
15enum SimplicialCholeskyMode {
16  SimplicialCholeskyLLT,
17  SimplicialCholeskyLDLT
18};
19
20/** \ingroup SparseCholesky_Module
21  * \brief A direct sparse Cholesky factorizations
22  *
23  * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are
24  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
25  * X and B can be either dense or sparse.
26  *
27  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
28  * such that the factorized matrix is P A P^-1.
29  *
30  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
31  * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
32  *               or Upper. Default is Lower.
33  *
34  */
35template<typename Derived>
36class SimplicialCholeskyBase : internal::noncopyable
37{
38  public:
39    typedef typename internal::traits<Derived>::MatrixType MatrixType;
40    typedef typename internal::traits<Derived>::OrderingType OrderingType;
41    enum { UpLo = internal::traits<Derived>::UpLo };
42    typedef typename MatrixType::Scalar Scalar;
43    typedef typename MatrixType::RealScalar RealScalar;
44    typedef typename MatrixType::Index Index;
45    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
46    typedef Matrix<Scalar,Dynamic,1> VectorType;
47
48  public:
49
50    /** Default constructor */
51    SimplicialCholeskyBase()
52      : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
53    {}
54
55    SimplicialCholeskyBase(const MatrixType& matrix)
56      : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
57    {
58      derived().compute(matrix);
59    }
60
61    ~SimplicialCholeskyBase()
62    {
63    }
64
65    Derived& derived() { return *static_cast<Derived*>(this); }
66    const Derived& derived() const { return *static_cast<const Derived*>(this); }
67
68    inline Index cols() const { return m_matrix.cols(); }
69    inline Index rows() const { return m_matrix.rows(); }
70
71    /** \brief Reports whether previous computation was successful.
72      *
73      * \returns \c Success if computation was succesful,
74      *          \c NumericalIssue if the matrix.appears to be negative.
75      */
76    ComputationInfo info() const
77    {
78      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
79      return m_info;
80    }
81
82    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
83      *
84      * \sa compute()
85      */
86    template<typename Rhs>
87    inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
88    solve(const MatrixBase<Rhs>& b) const
89    {
90      eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
91      eigen_assert(rows()==b.rows()
92                && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
93      return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
94    }
95
96    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
97      *
98      * \sa compute()
99      */
100    template<typename Rhs>
101    inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
102    solve(const SparseMatrixBase<Rhs>& b) const
103    {
104      eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
105      eigen_assert(rows()==b.rows()
106                && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
107      return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
108    }
109
110    /** \returns the permutation P
111      * \sa permutationPinv() */
112    const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
113    { return m_P; }
114
115    /** \returns the inverse P^-1 of the permutation P
116      * \sa permutationP() */
117    const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
118    { return m_Pinv; }
119
120    /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
121      *
122      * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n
123      * \c d_ii = \a offset + \a scale * \c d_ii
124      *
125      * The default is the identity transformation with \a offset=0, and \a scale=1.
126      *
127      * \returns a reference to \c *this.
128      */
129    Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
130    {
131      m_shiftOffset = offset;
132      m_shiftScale = scale;
133      return derived();
134    }
135
136#ifndef EIGEN_PARSED_BY_DOXYGEN
137    /** \internal */
138    template<typename Stream>
139    void dumpMemory(Stream& s)
140    {
141      int total = 0;
142      s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
143      s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
144      s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
145      s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
146      s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
147      s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
148      s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
149    }
150
151    /** \internal */
152    template<typename Rhs,typename Dest>
153    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
154    {
155      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
156      eigen_assert(m_matrix.rows()==b.rows());
157
158      if(m_info!=Success)
159        return;
160
161      if(m_P.size()>0)
162        dest = m_P * b;
163      else
164        dest = b;
165
166      if(m_matrix.nonZeros()>0) // otherwise L==I
167        derived().matrixL().solveInPlace(dest);
168
169      if(m_diag.size()>0)
170        dest = m_diag.asDiagonal().inverse() * dest;
171
172      if (m_matrix.nonZeros()>0) // otherwise U==I
173        derived().matrixU().solveInPlace(dest);
174
175      if(m_P.size()>0)
176        dest = m_Pinv * dest;
177    }
178
179#endif // EIGEN_PARSED_BY_DOXYGEN
180
181  protected:
182
183    /** Computes the sparse Cholesky decomposition of \a matrix */
184    template<bool DoLDLT>
185    void compute(const MatrixType& matrix)
186    {
187      eigen_assert(matrix.rows()==matrix.cols());
188      Index size = matrix.cols();
189      CholMatrixType ap(size,size);
190      ordering(matrix, ap);
191      analyzePattern_preordered(ap, DoLDLT);
192      factorize_preordered<DoLDLT>(ap);
193    }
194
195    template<bool DoLDLT>
196    void factorize(const MatrixType& a)
197    {
198      eigen_assert(a.rows()==a.cols());
199      int size = a.cols();
200      CholMatrixType ap(size,size);
201      ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
202      factorize_preordered<DoLDLT>(ap);
203    }
204
205    template<bool DoLDLT>
206    void factorize_preordered(const CholMatrixType& a);
207
208    void analyzePattern(const MatrixType& a, bool doLDLT)
209    {
210      eigen_assert(a.rows()==a.cols());
211      int size = a.cols();
212      CholMatrixType ap(size,size);
213      ordering(a, ap);
214      analyzePattern_preordered(ap,doLDLT);
215    }
216    void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
217
218    void ordering(const MatrixType& a, CholMatrixType& ap);
219
220    /** keeps off-diagonal entries; drops diagonal entries */
221    struct keep_diag {
222      inline bool operator() (const Index& row, const Index& col, const Scalar&) const
223      {
224        return row!=col;
225      }
226    };
227
228    mutable ComputationInfo m_info;
229    bool m_isInitialized;
230    bool m_factorizationIsOk;
231    bool m_analysisIsOk;
232
233    CholMatrixType m_matrix;
234    VectorType m_diag;                                // the diagonal coefficients (LDLT mode)
235    VectorXi m_parent;                                // elimination tree
236    VectorXi m_nonZerosPerCol;
237    PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
238    PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation
239
240    RealScalar m_shiftOffset;
241    RealScalar m_shiftScale;
242};
243
244template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
245template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
246template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
247
248namespace internal {
249
250template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
251{
252  typedef _MatrixType MatrixType;
253  typedef _Ordering OrderingType;
254  enum { UpLo = _UpLo };
255  typedef typename MatrixType::Scalar                         Scalar;
256  typedef typename MatrixType::Index                          Index;
257  typedef SparseMatrix<Scalar, ColMajor, Index>               CholMatrixType;
258  typedef SparseTriangularView<CholMatrixType, Eigen::Lower>  MatrixL;
259  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper>   MatrixU;
260  static inline MatrixL getL(const MatrixType& m) { return m; }
261  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
262};
263
264template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
265{
266  typedef _MatrixType MatrixType;
267  typedef _Ordering OrderingType;
268  enum { UpLo = _UpLo };
269  typedef typename MatrixType::Scalar                             Scalar;
270  typedef typename MatrixType::Index                              Index;
271  typedef SparseMatrix<Scalar, ColMajor, Index>                   CholMatrixType;
272  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower>  MatrixL;
273  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
274  static inline MatrixL getL(const MatrixType& m) { return m; }
275  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
276};
277
278template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
279{
280  typedef _MatrixType MatrixType;
281  typedef _Ordering OrderingType;
282  enum { UpLo = _UpLo };
283};
284
285}
286
287/** \ingroup SparseCholesky_Module
288  * \class SimplicialLLT
289  * \brief A direct sparse LLT Cholesky factorizations
290  *
291  * This class provides a LL^T Cholesky factorizations of sparse matrices that are
292  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
293  * X and B can be either dense or sparse.
294  *
295  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
296  * such that the factorized matrix is P A P^-1.
297  *
298  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
299  * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
300  *               or Upper. Default is Lower.
301  * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
302  *
303  * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
304  */
305template<typename _MatrixType, int _UpLo, typename _Ordering>
306    class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
307{
308public:
309    typedef _MatrixType MatrixType;
310    enum { UpLo = _UpLo };
311    typedef SimplicialCholeskyBase<SimplicialLLT> Base;
312    typedef typename MatrixType::Scalar Scalar;
313    typedef typename MatrixType::RealScalar RealScalar;
314    typedef typename MatrixType::Index Index;
315    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
316    typedef Matrix<Scalar,Dynamic,1> VectorType;
317    typedef internal::traits<SimplicialLLT> Traits;
318    typedef typename Traits::MatrixL  MatrixL;
319    typedef typename Traits::MatrixU  MatrixU;
320public:
321    /** Default constructor */
322    SimplicialLLT() : Base() {}
323    /** Constructs and performs the LLT factorization of \a matrix */
324    SimplicialLLT(const MatrixType& matrix)
325        : Base(matrix) {}
326
327    /** \returns an expression of the factor L */
328    inline const MatrixL matrixL() const {
329        eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
330        return Traits::getL(Base::m_matrix);
331    }
332
333    /** \returns an expression of the factor U (= L^*) */
334    inline const MatrixU matrixU() const {
335        eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
336        return Traits::getU(Base::m_matrix);
337    }
338
339    /** Computes the sparse Cholesky decomposition of \a matrix */
340    SimplicialLLT& compute(const MatrixType& matrix)
341    {
342      Base::template compute<false>(matrix);
343      return *this;
344    }
345
346    /** Performs a symbolic decomposition on the sparcity of \a matrix.
347      *
348      * This function is particularly useful when solving for several problems having the same structure.
349      *
350      * \sa factorize()
351      */
352    void analyzePattern(const MatrixType& a)
353    {
354      Base::analyzePattern(a, false);
355    }
356
357    /** Performs a numeric decomposition of \a matrix
358      *
359      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
360      *
361      * \sa analyzePattern()
362      */
363    void factorize(const MatrixType& a)
364    {
365      Base::template factorize<false>(a);
366    }
367
368    /** \returns the determinant of the underlying matrix from the current factorization */
369    Scalar determinant() const
370    {
371      Scalar detL = Base::m_matrix.diagonal().prod();
372      return numext::abs2(detL);
373    }
374};
375
376/** \ingroup SparseCholesky_Module
377  * \class SimplicialLDLT
378  * \brief A direct sparse LDLT Cholesky factorizations without square root.
379  *
380  * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
381  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
382  * X and B can be either dense or sparse.
383  *
384  * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
385  * such that the factorized matrix is P A P^-1.
386  *
387  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
388  * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
389  *               or Upper. Default is Lower.
390  * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
391  *
392  * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
393  */
394template<typename _MatrixType, int _UpLo, typename _Ordering>
395    class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
396{
397public:
398    typedef _MatrixType MatrixType;
399    enum { UpLo = _UpLo };
400    typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
401    typedef typename MatrixType::Scalar Scalar;
402    typedef typename MatrixType::RealScalar RealScalar;
403    typedef typename MatrixType::Index Index;
404    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
405    typedef Matrix<Scalar,Dynamic,1> VectorType;
406    typedef internal::traits<SimplicialLDLT> Traits;
407    typedef typename Traits::MatrixL  MatrixL;
408    typedef typename Traits::MatrixU  MatrixU;
409public:
410    /** Default constructor */
411    SimplicialLDLT() : Base() {}
412
413    /** Constructs and performs the LLT factorization of \a matrix */
414    SimplicialLDLT(const MatrixType& matrix)
415        : Base(matrix) {}
416
417    /** \returns a vector expression of the diagonal D */
418    inline const VectorType vectorD() const {
419        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
420        return Base::m_diag;
421    }
422    /** \returns an expression of the factor L */
423    inline const MatrixL matrixL() const {
424        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
425        return Traits::getL(Base::m_matrix);
426    }
427
428    /** \returns an expression of the factor U (= L^*) */
429    inline const MatrixU matrixU() const {
430        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
431        return Traits::getU(Base::m_matrix);
432    }
433
434    /** Computes the sparse Cholesky decomposition of \a matrix */
435    SimplicialLDLT& compute(const MatrixType& matrix)
436    {
437      Base::template compute<true>(matrix);
438      return *this;
439    }
440
441    /** Performs a symbolic decomposition on the sparcity of \a matrix.
442      *
443      * This function is particularly useful when solving for several problems having the same structure.
444      *
445      * \sa factorize()
446      */
447    void analyzePattern(const MatrixType& a)
448    {
449      Base::analyzePattern(a, true);
450    }
451
452    /** Performs a numeric decomposition of \a matrix
453      *
454      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
455      *
456      * \sa analyzePattern()
457      */
458    void factorize(const MatrixType& a)
459    {
460      Base::template factorize<true>(a);
461    }
462
463    /** \returns the determinant of the underlying matrix from the current factorization */
464    Scalar determinant() const
465    {
466      return Base::m_diag.prod();
467    }
468};
469
470/** \deprecated use SimplicialLDLT or class SimplicialLLT
471  * \ingroup SparseCholesky_Module
472  * \class SimplicialCholesky
473  *
474  * \sa class SimplicialLDLT, class SimplicialLLT
475  */
476template<typename _MatrixType, int _UpLo, typename _Ordering>
477    class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
478{
479public:
480    typedef _MatrixType MatrixType;
481    enum { UpLo = _UpLo };
482    typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
483    typedef typename MatrixType::Scalar Scalar;
484    typedef typename MatrixType::RealScalar RealScalar;
485    typedef typename MatrixType::Index Index;
486    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
487    typedef Matrix<Scalar,Dynamic,1> VectorType;
488    typedef internal::traits<SimplicialCholesky> Traits;
489    typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
490    typedef internal::traits<SimplicialLLT<MatrixType,UpLo>  > LLTTraits;
491  public:
492    SimplicialCholesky() : Base(), m_LDLT(true) {}
493
494    SimplicialCholesky(const MatrixType& matrix)
495      : Base(), m_LDLT(true)
496    {
497      compute(matrix);
498    }
499
500    SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
501    {
502      switch(mode)
503      {
504      case SimplicialCholeskyLLT:
505        m_LDLT = false;
506        break;
507      case SimplicialCholeskyLDLT:
508        m_LDLT = true;
509        break;
510      default:
511        break;
512      }
513
514      return *this;
515    }
516
517    inline const VectorType vectorD() const {
518        eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
519        return Base::m_diag;
520    }
521    inline const CholMatrixType rawMatrix() const {
522        eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
523        return Base::m_matrix;
524    }
525
526    /** Computes the sparse Cholesky decomposition of \a matrix */
527    SimplicialCholesky& compute(const MatrixType& matrix)
528    {
529      if(m_LDLT)
530        Base::template compute<true>(matrix);
531      else
532        Base::template compute<false>(matrix);
533      return *this;
534    }
535
536    /** Performs a symbolic decomposition on the sparcity of \a matrix.
537      *
538      * This function is particularly useful when solving for several problems having the same structure.
539      *
540      * \sa factorize()
541      */
542    void analyzePattern(const MatrixType& a)
543    {
544      Base::analyzePattern(a, m_LDLT);
545    }
546
547    /** Performs a numeric decomposition of \a matrix
548      *
549      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
550      *
551      * \sa analyzePattern()
552      */
553    void factorize(const MatrixType& a)
554    {
555      if(m_LDLT)
556        Base::template factorize<true>(a);
557      else
558        Base::template factorize<false>(a);
559    }
560
561    /** \internal */
562    template<typename Rhs,typename Dest>
563    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
564    {
565      eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
566      eigen_assert(Base::m_matrix.rows()==b.rows());
567
568      if(Base::m_info!=Success)
569        return;
570
571      if(Base::m_P.size()>0)
572        dest = Base::m_P * b;
573      else
574        dest = b;
575
576      if(Base::m_matrix.nonZeros()>0) // otherwise L==I
577      {
578        if(m_LDLT)
579          LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
580        else
581          LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
582      }
583
584      if(Base::m_diag.size()>0)
585        dest = Base::m_diag.asDiagonal().inverse() * dest;
586
587      if (Base::m_matrix.nonZeros()>0) // otherwise I==I
588      {
589        if(m_LDLT)
590          LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
591        else
592          LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
593      }
594
595      if(Base::m_P.size()>0)
596        dest = Base::m_Pinv * dest;
597    }
598
599    Scalar determinant() const
600    {
601      if(m_LDLT)
602      {
603        return Base::m_diag.prod();
604      }
605      else
606      {
607        Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
608        return numext::abs2(detL);
609      }
610    }
611
612  protected:
613    bool m_LDLT;
614};
615
616template<typename Derived>
617void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
618{
619  eigen_assert(a.rows()==a.cols());
620  const Index size = a.rows();
621  // Note that amd compute the inverse permutation
622  {
623    CholMatrixType C;
624    C = a.template selfadjointView<UpLo>();
625
626    OrderingType ordering;
627    ordering(C,m_Pinv);
628  }
629
630  if(m_Pinv.size()>0)
631    m_P = m_Pinv.inverse();
632  else
633    m_P.resize(0);
634
635  ap.resize(size,size);
636  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
637}
638
639namespace internal {
640
641template<typename Derived, typename Rhs>
642struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
643  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
644{
645  typedef SimplicialCholeskyBase<Derived> Dec;
646  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
647
648  template<typename Dest> void evalTo(Dest& dst) const
649  {
650    dec().derived()._solve(rhs(),dst);
651  }
652};
653
654template<typename Derived, typename Rhs>
655struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
656  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
657{
658  typedef SimplicialCholeskyBase<Derived> Dec;
659  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
660
661  template<typename Dest> void evalTo(Dest& dst) const
662  {
663    this->defaultEvalTo(dst);
664  }
665};
666
667} // end namespace internal
668
669} // end namespace Eigen
670
671#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
672