Eigen/src/SparseCholesky/SimplicialCholesky.h - fp2-dev/platform/external/eigen - Gitiles

 // 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
 };

 /** \ingroup SparseCholesky_Module
   * \brief A direct sparse Cholesky factorizations
   *
   * These classes provide LL^T and LDL^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.
   *
   */
 template<typename Derived>
 class SimplicialCholeskyBase : internal::noncopyable
 {
   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::Index Index;
     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
     typedef Matrix<Scalar,Dynamic,1> VectorType;

   public:

     /** Default constructor */
     SimplicialCholeskyBase()
       : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
     {}

     SimplicialCholeskyBase(const MatrixType& matrix)
       : m_info(Success), m_isInitialized(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 succesful,
       *          \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 solution x of \f$ A x = b \f$ using the current decomposition of A.
       *
       * \sa compute()
       */
     template<typename Rhs>
     inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
     solve(const MatrixBase<Rhs>& b) const
     {
       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
       eigen_assert(rows()==b.rows()
                 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
       return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
     }

     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
       *
       * \sa compute()
       */
     template<typename Rhs>
     inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
     solve(const SparseMatrixBase<Rhs>& b) const
     {
       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
       eigen_assert(rows()==b.rows()
                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
       return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
     }

     /** \returns the permutation P
       * \sa permutationPinv() */
     const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
     { return m_P; }

     /** \returns the inverse P^-1 of the permutation P
       * \sa permutationP() */
     const PermutationMatrix<Dynamic,Dynamic,Index>& 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(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;
     }

 #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 ap(size,size);
       ordering(matrix, ap);
       analyzePattern_preordered(ap, DoLDLT);
       factorize_preordered<DoLDLT>(ap);
     }

     template<bool DoLDLT>
     void factorize(const MatrixType& a)
     {
       eigen_assert(a.rows()==a.cols());
       int size = a.cols();
       CholMatrixType ap(size,size);
       ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
       factorize_preordered<DoLDLT>(ap);
     }

     template<bool DoLDLT>
     void factorize_preordered(const CholMatrixType& a);

     void analyzePattern(const MatrixType& a, bool doLDLT)
     {
       eigen_assert(a.rows()==a.cols());
       int size = a.cols();
       CholMatrixType ap(size,size);
       ordering(a, ap);
       analyzePattern_preordered(ap,doLDLT);
     }
     void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);

     void ordering(const MatrixType& a, 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_isInitialized;
     bool m_factorizationIsOk;
     bool m_analysisIsOk;

     CholMatrixType m_matrix;
     VectorType m_diag;                                // the diagonal coefficients (LDLT mode)
     VectorXi m_parent;                                // elimination tree
     VectorXi m_nonZerosPerCol;
     PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
     PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation

     RealScalar m_shiftOffset;
     RealScalar m_shiftScale;
 };

 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > 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::Index                          Index;
   typedef SparseMatrix<Scalar, ColMajor, Index>               CholMatrixType;
   typedef SparseTriangularView<CholMatrixType, Eigen::Lower>  MatrixL;
   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper>   MatrixU;
   static inline MatrixL getL(const MatrixType& m) { return m; }
   static inline MatrixU getU(const MatrixType& m) { return 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::Index                              Index;
   typedef SparseMatrix<Scalar, ColMajor, Index>                   CholMatrixType;
   typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower>  MatrixL;
   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
   static inline MatrixL getL(const MatrixType& m) { return m; }
   static inline MatrixU getU(const MatrixType& m) { return 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<>
   *
   * \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::Index Index;
     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 */
     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<>
   *
   * \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::Index Index;
     typedef SparseMatrix<Scalar,ColMajor,Index> 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 */
     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::Index Index;
     typedef SparseMatrix<Scalar,ColMajor,Index> 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) {}

     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(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;
     }

     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, CholMatrixType& ap)
 {
   eigen_assert(a.rows()==a.cols());
   const Index size = a.rows();
   // Note that amd compute the inverse permutation
   {
     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);
 }

 namespace internal {

 template<typename Derived, typename Rhs>
 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
   : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
 {
   typedef SimplicialCholeskyBase<Derived> Dec;
   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)

   template<typename Dest> void evalTo(Dest& dst) const
   {
     dec().derived()._solve(rhs(),dst);
   }
 };

 template<typename Derived, typename Rhs>
 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
   : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
 {
   typedef SimplicialCholeskyBase<Derived> Dec;
   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)

   template<typename Dest> void evalTo(Dest& dst) const
   {
     this->defaultEvalTo(dst);
   }
 };

 } // end namespace internal

 } // end namespace Eigen

 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
	// 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
	};

	/** \ingroup SparseCholesky_Module
	* \brief A direct sparse Cholesky factorizations
	*
	* These classes provide LL^T and LDL^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.
	*
	*/
	template<typename Derived>
	class SimplicialCholeskyBase : internal::noncopyable
	{
	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::Index Index;
	typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
	typedef Matrix<Scalar,Dynamic,1> VectorType;

	public:

	/** Default constructor */
	SimplicialCholeskyBase()
	: m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
	{}

	SimplicialCholeskyBase(const MatrixType& matrix)
	: m_info(Success), m_isInitialized(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 succesful,
	* \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 solution x of \f$ A x = b \f$ using the current decomposition of A.
	*
	* \sa compute()
	*/
	template<typename Rhs>
	inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
	solve(const MatrixBase<Rhs>& b) const
	{
	eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
	eigen_assert(rows()==b.rows()
	&& "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
	return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
	}

	/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
	*
	* \sa compute()
	*/
	template<typename Rhs>
	inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
	solve(const SparseMatrixBase<Rhs>& b) const
	{
	eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
	eigen_assert(rows()==b.rows()
	&& "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
	return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
	}

	/** \returns the permutation P
	* \sa permutationPinv() */
	const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
	{ return m_P; }

	/** \returns the inverse P^-1 of the permutation P
	* \sa permutationP() */
	const PermutationMatrix<Dynamic,Dynamic,Index>& 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(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;
	}

	#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 ap(size,size);
	ordering(matrix, ap);
	analyzePattern_preordered(ap, DoLDLT);
	factorize_preordered<DoLDLT>(ap);
	}

	template<bool DoLDLT>
	void factorize(const MatrixType& a)
	{
	eigen_assert(a.rows()==a.cols());
	int size = a.cols();
	CholMatrixType ap(size,size);
	ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
	factorize_preordered<DoLDLT>(ap);
	}

	template<bool DoLDLT>
	void factorize_preordered(const CholMatrixType& a);

	void analyzePattern(const MatrixType& a, bool doLDLT)
	{
	eigen_assert(a.rows()==a.cols());
	int size = a.cols();
	CholMatrixType ap(size,size);
	ordering(a, ap);
	analyzePattern_preordered(ap,doLDLT);
	}
	void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);

	void ordering(const MatrixType& a, 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_isInitialized;
	bool m_factorizationIsOk;
	bool m_analysisIsOk;

	CholMatrixType m_matrix;
	VectorType m_diag; // the diagonal coefficients (LDLT mode)
	VectorXi m_parent; // elimination tree
	VectorXi m_nonZerosPerCol;
	PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
	PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation

	RealScalar m_shiftOffset;
	RealScalar m_shiftScale;
	};

	template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
	template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
	template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > 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::Index Index;
	typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
	typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
	typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
	static inline MatrixL getL(const MatrixType& m) { return m; }
	static inline MatrixU getU(const MatrixType& m) { return 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::Index Index;
	typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
	typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
	typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
	static inline MatrixL getL(const MatrixType& m) { return m; }
	static inline MatrixU getU(const MatrixType& m) { return 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<>
	*
	* \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::Index Index;
	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 */
	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<>
	*
	* \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::Index Index;
	typedef SparseMatrix<Scalar,ColMajor,Index> 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 */
	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::Index Index;
	typedef SparseMatrix<Scalar,ColMajor,Index> 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) {}

	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(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;
	}

	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, CholMatrixType& ap)
	{
	eigen_assert(a.rows()==a.cols());
	const Index size = a.rows();
	// Note that amd compute the inverse permutation
	{
	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);
	}

	namespace internal {

	template<typename Derived, typename Rhs>
	struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
	: solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
	{
	typedef SimplicialCholeskyBase<Derived> Dec;
	EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)

	template<typename Dest> void evalTo(Dest& dst) const
	{
	dec().derived()._solve(rhs(),dst);
	}
	};

	template<typename Derived, typename Rhs>
	struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
	: sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
	{
	typedef SimplicialCholeskyBase<Derived> Dec;
	EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)

	template<typename Dest> void evalTo(Dest& dst) const
	{
	this->defaultEvalTo(dst);
	}
	};

	} // end namespace internal

	} // end namespace Eigen

	#endif // EIGEN_SIMPLICIAL_CHOLESKY_H