Eigen/src/CholmodSupport/CholmodSupport.h - mirror - Git at Google

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2008-2010 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_CHOLMODSUPPORT_H
 #define EIGEN_CHOLMODSUPPORT_H

 namespace Eigen {

 namespace internal {

 template<typename Scalar, typename CholmodType>
 void cholmod_configure_matrix(CholmodType& mat)
 {
   if (internal::is_same<Scalar,float>::value)
   {
     mat.xtype = CHOLMOD_REAL;
     mat.dtype = CHOLMOD_SINGLE;
   }
   else if (internal::is_same<Scalar,double>::value)
   {
     mat.xtype = CHOLMOD_REAL;
     mat.dtype = CHOLMOD_DOUBLE;
   }
   else if (internal::is_same<Scalar,std::complex<float> >::value)
   {
     mat.xtype = CHOLMOD_COMPLEX;
     mat.dtype = CHOLMOD_SINGLE;
   }
   else if (internal::is_same<Scalar,std::complex<double> >::value)
   {
     mat.xtype = CHOLMOD_COMPLEX;
     mat.dtype = CHOLMOD_DOUBLE;
   }
   else
   {
     eigen_assert(false && "Scalar type not supported by CHOLMOD");
   }
 }

 } // namespace internal

 /** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
   * Note that the data are shared.
   */
 template<typename _Scalar, int _Options, typename _StorageIndex>
 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_StorageIndex>& mat)
 {
   cholmod_sparse res;
   res.nzmax   = mat.nonZeros();
   res.nrow    = mat.rows();;
   res.ncol    = mat.cols();
   res.p       = mat.outerIndexPtr();
   res.i       = mat.innerIndexPtr();
   res.x       = mat.valuePtr();
   res.z       = 0;
   res.sorted  = 1;
   if(mat.isCompressed())
   {
     res.packed  = 1;
     res.nz = 0;
   }
   else
   {
     res.packed  = 0;
     res.nz = mat.innerNonZeroPtr();
   }

   res.dtype   = 0;
   res.stype   = -1;

   if (internal::is_same<_StorageIndex,int>::value)
   {
     res.itype = CHOLMOD_INT;
   }
   else if (internal::is_same<_StorageIndex,UF_long>::value)
   {
     res.itype = CHOLMOD_LONG;
   }
   else
   {
     eigen_assert(false && "Index type not supported yet");
   }

   // setup res.xtype
   internal::cholmod_configure_matrix<_Scalar>(res);

   res.stype = 0;

   return res;
 }

 template<typename _Scalar, int _Options, typename _Index>
 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
 {
   cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
   return res;
 }

 /** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
   * The data are not copied but shared. */
 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
 {
   cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());

   if(UpLo==Upper) res.stype =  1;
   if(UpLo==Lower) res.stype = -1;

   return res;
 }

 /** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
   * The data are not copied but shared. */
 template<typename Derived>
 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
 {
   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
   typedef typename Derived::Scalar Scalar;

   cholmod_dense res;
   res.nrow   = mat.rows();
   res.ncol   = mat.cols();
   res.nzmax  = res.nrow * res.ncol;
   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
   res.x      = (void*)(mat.derived().data());
   res.z      = 0;

   internal::cholmod_configure_matrix<Scalar>(res);

   return res;
 }

 /** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
   * The data are not copied but shared. */
 template<typename Scalar, int Flags, typename StorageIndex>
 MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm)
 {
   return MappedSparseMatrix<Scalar,Flags,StorageIndex>
          (cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol],
           static_cast<StorageIndex*>(cm.p), static_cast<StorageIndex*>(cm.i),static_cast<Scalar*>(cm.x) );
 }

 enum CholmodMode {
   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
 };


 /** \ingroup CholmodSupport_Module
   * \class CholmodBase
   * \brief The base class for the direct Cholesky factorization of Cholmod
   * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT
   */
 template<typename _MatrixType, int _UpLo, typename Derived>
 class CholmodBase : public SparseSolverBase<Derived>
 {
   protected:
     typedef SparseSolverBase<Derived> Base;
     using Base::derived;
     using Base::m_isInitialized;
   public:
     typedef _MatrixType MatrixType;
     enum { UpLo = _UpLo };
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef MatrixType CholMatrixType;
     typedef typename MatrixType::StorageIndex StorageIndex;

   public:

     CholmodBase()
       : m_cholmodFactor(0), m_info(Success)
     {
       m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
       cholmod_start(&m_cholmod);
     }

     explicit CholmodBase(const MatrixType& matrix)
       : m_cholmodFactor(0), m_info(Success)
     {
       m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
       cholmod_start(&m_cholmod);
       compute(matrix);
     }

     ~CholmodBase()
     {
       if(m_cholmodFactor)
         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
       cholmod_finish(&m_cholmod);
     }

     inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
     inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }

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

     /** Computes the sparse Cholesky decomposition of \a matrix */
     Derived& compute(const MatrixType& matrix)
     {
       analyzePattern(matrix);
       factorize(matrix);
       return derived();
     }

     /** Performs a symbolic decomposition on the sparsity pattern of \a matrix.
       *
       * This function is particularly useful when solving for several problems having the same structure.
       *
       * \sa factorize()
       */
     void analyzePattern(const MatrixType& matrix)
     {
       if(m_cholmodFactor)
       {
         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
         m_cholmodFactor = 0;
       }
       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);

       this->m_isInitialized = true;
       this->m_info = Success;
       m_analysisIsOk = true;
       m_factorizationIsOk = false;
     }

     /** Performs a numeric decomposition of \a matrix
       *
       * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been performed.
       *
       * \sa analyzePattern()
       */
     void factorize(const MatrixType& matrix)
     {
       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
       cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);

       // If the factorization failed, minor is the column at which it did. On success minor == n.
       this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
       m_factorizationIsOk = true;
     }

     /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
      *  See the Cholmod user guide for details. */
     cholmod_common& cholmod() { return m_cholmod; }

     #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal */
     template<typename Rhs,typename Dest>
     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
     {
       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
       const Index size = m_cholmodFactor->n;
       EIGEN_UNUSED_VARIABLE(size);
       eigen_assert(size==b.rows());

       // note: cd stands for Cholmod Dense
       Rhs& b_ref(b.const_cast_derived());
       cholmod_dense b_cd = viewAsCholmod(b_ref);
       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
       if(!x_cd)
       {
         this->m_info = NumericalIssue;
         return;
       }
       // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
       cholmod_free_dense(&x_cd, &m_cholmod);
     }

     /** \internal */
     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
     void _solve_impl(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &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()");
       const Index size = m_cholmodFactor->n;
       EIGEN_UNUSED_VARIABLE(size);
       eigen_assert(size==b.rows());

       // note: cs stands for Cholmod Sparse
       cholmod_sparse b_cs = viewAsCholmod(b);
       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
       if(!x_cs)
       {
         this->m_info = NumericalIssue;
         return;
       }
       // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
       dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
       cholmod_free_sparse(&x_cs, &m_cholmod);
     }
     #endif // EIGEN_PARSED_BY_DOXYGEN


     /** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization.
       *
       * During the numerical factorization, an offset term is added to the diagonal coefficients:\n
       * \c d_ii = \a offset + \c d_ii
       *
       * The default is \a offset=0.
       *
       * \returns a reference to \c *this.
       */
     Derived& setShift(const RealScalar& offset)
     {
       m_shiftOffset[0] = offset;
       return derived();
     }

     template<typename Stream>
     void dumpMemory(Stream& /*s*/)
     {}

   protected:
     mutable cholmod_common m_cholmod;
     cholmod_factor* m_cholmodFactor;
     RealScalar m_shiftOffset[2];
     mutable ComputationInfo m_info;
     int m_factorizationIsOk;
     int m_analysisIsOk;
 };

 /** \ingroup CholmodSupport_Module
   * \class CholmodSimplicialLLT
   * \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod
   *
   * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization
   * using the Cholmod library.
   * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical interest.
   * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
   * X and B can be either dense or sparse.
   *
   * \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.
   *
   * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
   *
   * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT
   */
 template<typename _MatrixType, int _UpLo = Lower>
 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
 {
     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
     using Base::m_cholmod;

   public:

     typedef _MatrixType MatrixType;

     CholmodSimplicialLLT() : Base() { init(); }

     CholmodSimplicialLLT(const MatrixType& matrix) : Base()
     {
       init();
       this->compute(matrix);
     }

     ~CholmodSimplicialLLT() {}
   protected:
     void init()
     {
       m_cholmod.final_asis = 0;
       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
       m_cholmod.final_ll = 1;
     }
 };


 /** \ingroup CholmodSupport_Module
   * \class CholmodSimplicialLDLT
   * \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod
   *
   * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization
   * using the Cholmod library.
   * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical interest.
   * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
   * X and B can be either dense or sparse.
   *
   * \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.
   *
   * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
   *
   * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT
   */
 template<typename _MatrixType, int _UpLo = Lower>
 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
 {
     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
     using Base::m_cholmod;

   public:

     typedef _MatrixType MatrixType;

     CholmodSimplicialLDLT() : Base() { init(); }

     CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
     {
       init();
       this->compute(matrix);
     }

     ~CholmodSimplicialLDLT() {}
   protected:
     void init()
     {
       m_cholmod.final_asis = 1;
       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
     }
 };

 /** \ingroup CholmodSupport_Module
   * \class CholmodSupernodalLLT
   * \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod
   *
   * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization
   * using the Cholmod library.
   * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM.
   * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
   * X and B can be either dense or sparse.
   *
   * \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.
   *
   * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
   *
   * \sa \ref TutorialSparseDirectSolvers
   */
 template<typename _MatrixType, int _UpLo = Lower>
 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
 {
     typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
     using Base::m_cholmod;

   public:

     typedef _MatrixType MatrixType;

     CholmodSupernodalLLT() : Base() { init(); }

     CholmodSupernodalLLT(const MatrixType& matrix) : Base()
     {
       init();
       this->compute(matrix);
     }

     ~CholmodSupernodalLLT() {}
   protected:
     void init()
     {
       m_cholmod.final_asis = 1;
       m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
     }
 };

 /** \ingroup CholmodSupport_Module
   * \class CholmodDecomposition
   * \brief A general Cholesky factorization and solver based on Cholmod
   *
   * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
   * using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
   * X and B can be either dense or sparse.
   *
   * This variant permits to change the underlying Cholesky method at runtime.
   * On the other hand, it does not provide access to the result of the factorization.
   * The default is to let Cholmod automatically choose between a simplicial and supernodal factorization.
   *
   * \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.
   *
   * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
   *
   * \sa \ref TutorialSparseDirectSolvers
   */
 template<typename _MatrixType, int _UpLo = Lower>
 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
 {
     typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
     using Base::m_cholmod;

   public:

     typedef _MatrixType MatrixType;

     CholmodDecomposition() : Base() { init(); }

     CholmodDecomposition(const MatrixType& matrix) : Base()
     {
       init();
       this->compute(matrix);
     }

     ~CholmodDecomposition() {}

     void setMode(CholmodMode mode)
     {
       switch(mode)
       {
         case CholmodAuto:
           m_cholmod.final_asis = 1;
           m_cholmod.supernodal = CHOLMOD_AUTO;
           break;
         case CholmodSimplicialLLt:
           m_cholmod.final_asis = 0;
           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
           m_cholmod.final_ll = 1;
           break;
         case CholmodSupernodalLLt:
           m_cholmod.final_asis = 1;
           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
           break;
         case CholmodLDLt:
           m_cholmod.final_asis = 1;
           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
           break;
         default:
           break;
       }
     }
   protected:
     void init()
     {
       m_cholmod.final_asis = 1;
       m_cholmod.supernodal = CHOLMOD_AUTO;
     }
 };

 } // end namespace Eigen

 #endif // EIGEN_CHOLMODSUPPORT_H
	// This file is part of Eigen, a lightweight C++ template library
	// for linear algebra.
	//
	// Copyright (C) 2008-2010 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_CHOLMODSUPPORT_H
	#define EIGEN_CHOLMODSUPPORT_H

	namespace Eigen {

	namespace internal {

	template<typename Scalar, typename CholmodType>
	void cholmod_configure_matrix(CholmodType& mat)
	{
	if (internal::is_same<Scalar,float>::value)
	{
	mat.xtype = CHOLMOD_REAL;
	mat.dtype = CHOLMOD_SINGLE;
	}
	else if (internal::is_same<Scalar,double>::value)
	{
	mat.xtype = CHOLMOD_REAL;
	mat.dtype = CHOLMOD_DOUBLE;
	}
	else if (internal::is_same<Scalar,std::complex<float> >::value)
	{
	mat.xtype = CHOLMOD_COMPLEX;
	mat.dtype = CHOLMOD_SINGLE;
	}
	else if (internal::is_same<Scalar,std::complex<double> >::value)
	{
	mat.xtype = CHOLMOD_COMPLEX;
	mat.dtype = CHOLMOD_DOUBLE;
	}
	else
	{
	eigen_assert(false && "Scalar type not supported by CHOLMOD");
	}
	}

	} // namespace internal

	/** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
	* Note that the data are shared.
	*/
	template<typename _Scalar, int _Options, typename _StorageIndex>
	cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_StorageIndex>& mat)
	{
	cholmod_sparse res;
	res.nzmax = mat.nonZeros();
	res.nrow = mat.rows();;
	res.ncol = mat.cols();
	res.p = mat.outerIndexPtr();
	res.i = mat.innerIndexPtr();
	res.x = mat.valuePtr();
	res.z = 0;
	res.sorted = 1;
	if(mat.isCompressed())
	{
	res.packed = 1;
	res.nz = 0;
	}
	else
	{
	res.packed = 0;
	res.nz = mat.innerNonZeroPtr();
	}

	res.dtype = 0;
	res.stype = -1;

	if (internal::is_same<_StorageIndex,int>::value)
	{
	res.itype = CHOLMOD_INT;
	}
	else if (internal::is_same<_StorageIndex,UF_long>::value)
	{
	res.itype = CHOLMOD_LONG;
	}
	else
	{
	eigen_assert(false && "Index type not supported yet");
	}

	// setup res.xtype
	internal::cholmod_configure_matrix<_Scalar>(res);

	res.stype = 0;

	return res;
	}

	template<typename _Scalar, int _Options, typename _Index>
	const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
	{
	cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
	return res;
	}

	/** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
	* The data are not copied but shared. */
	template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
	cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
	{
	cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());

	if(UpLo==Upper) res.stype = 1;
	if(UpLo==Lower) res.stype = -1;

	return res;
	}

	/** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
	* The data are not copied but shared. */
	template<typename Derived>
	cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
	{
	EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
	typedef typename Derived::Scalar Scalar;

	cholmod_dense res;
	res.nrow = mat.rows();
	res.ncol = mat.cols();
	res.nzmax = res.nrow * res.ncol;
	res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
	res.x = (void*)(mat.derived().data());
	res.z = 0;

	internal::cholmod_configure_matrix<Scalar>(res);

	return res;
	}

	/** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
	* The data are not copied but shared. */
	template<typename Scalar, int Flags, typename StorageIndex>
	MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm)
	{
	return MappedSparseMatrix<Scalar,Flags,StorageIndex>
	(cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol],
	static_cast<StorageIndex>(cm.p), static_cast<StorageIndex>(cm.i),static_cast<Scalar*>(cm.x) );
	}

	enum CholmodMode {
	CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
	};


	/** \ingroup CholmodSupport_Module
	* \class CholmodBase
	* \brief The base class for the direct Cholesky factorization of Cholmod
	* \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT
	*/
	template<typename _MatrixType, int _UpLo, typename Derived>
	class CholmodBase : public SparseSolverBase<Derived>
	{
	protected:
	typedef SparseSolverBase<Derived> Base;
	using Base::derived;
	using Base::m_isInitialized;
	public:
	typedef _MatrixType MatrixType;
	enum { UpLo = _UpLo };
	typedef typename MatrixType::Scalar Scalar;
	typedef typename MatrixType::RealScalar RealScalar;
	typedef MatrixType CholMatrixType;
	typedef typename MatrixType::StorageIndex StorageIndex;

	public:

	CholmodBase()
	: m_cholmodFactor(0), m_info(Success)
	{
	m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
	cholmod_start(&m_cholmod);
	}

	explicit CholmodBase(const MatrixType& matrix)
	: m_cholmodFactor(0), m_info(Success)
	{
	m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
	cholmod_start(&m_cholmod);
	compute(matrix);
	}

	~CholmodBase()
	{
	if(m_cholmodFactor)
	cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
	cholmod_finish(&m_cholmod);
	}

	inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
	inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }

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

	/** Computes the sparse Cholesky decomposition of \a matrix */
	Derived& compute(const MatrixType& matrix)
	{
	analyzePattern(matrix);
	factorize(matrix);
	return derived();
	}

	/** Performs a symbolic decomposition on the sparsity pattern of \a matrix.
	*
	* This function is particularly useful when solving for several problems having the same structure.
	*
	* \sa factorize()
	*/
	void analyzePattern(const MatrixType& matrix)
	{
	if(m_cholmodFactor)
	{
	cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
	m_cholmodFactor = 0;
	}
	cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
	m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);

	this->m_isInitialized = true;
	this->m_info = Success;
	m_analysisIsOk = true;
	m_factorizationIsOk = false;
	}

	/** Performs a numeric decomposition of \a matrix
	*
	* The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been performed.
	*
	* \sa analyzePattern()
	*/
	void factorize(const MatrixType& matrix)
	{
	eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
	cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
	cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);

	// If the factorization failed, minor is the column at which it did. On success minor == n.
	this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
	m_factorizationIsOk = true;
	}

	/** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
	* See the Cholmod user guide for details. */
	cholmod_common& cholmod() { return m_cholmod; }

	#ifndef EIGEN_PARSED_BY_DOXYGEN
	/** \internal */
	template<typename Rhs,typename Dest>
	void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
	{
	eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
	const Index size = m_cholmodFactor->n;
	EIGEN_UNUSED_VARIABLE(size);
	eigen_assert(size==b.rows());

	// note: cd stands for Cholmod Dense
	Rhs& b_ref(b.const_cast_derived());
	cholmod_dense b_cd = viewAsCholmod(b_ref);
	cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
	if(!x_cd)
	{
	this->m_info = NumericalIssue;
	return;
	}
	// TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
	dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
	cholmod_free_dense(&x_cd, &m_cholmod);
	}

	/** \internal */
	template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
	void _solve_impl(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &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()");
	const Index size = m_cholmodFactor->n;
	EIGEN_UNUSED_VARIABLE(size);
	eigen_assert(size==b.rows());

	// note: cs stands for Cholmod Sparse
	cholmod_sparse b_cs = viewAsCholmod(b);
	cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
	if(!x_cs)
	{
	this->m_info = NumericalIssue;
	return;
	}
	// TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
	dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
	cholmod_free_sparse(&x_cs, &m_cholmod);
	}
	#endif // EIGEN_PARSED_BY_DOXYGEN


	/** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization.
	*
	* During the numerical factorization, an offset term is added to the diagonal coefficients:\n
	* \c d_ii = \a offset + \c d_ii
	*
	* The default is \a offset=0.
	*
	* \returns a reference to \c *this.
	*/
	Derived& setShift(const RealScalar& offset)
	{
	m_shiftOffset[0] = offset;
	return derived();
	}

	template<typename Stream>
	void dumpMemory(Stream& /s/)
	{}

	protected:
	mutable cholmod_common m_cholmod;
	cholmod_factor* m_cholmodFactor;
	RealScalar m_shiftOffset[2];
	mutable ComputationInfo m_info;
	int m_factorizationIsOk;
	int m_analysisIsOk;
	};

	/** \ingroup CholmodSupport_Module
	* \class CholmodSimplicialLLT
	* \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod
	*
	* This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization
	* using the Cholmod library.
	* This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical interest.
	* The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
	* X and B can be either dense or sparse.
	*
	* \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.
	*
	* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
	*
	* \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT
	*/
	template<typename _MatrixType, int _UpLo = Lower>
	class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
	{
	typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
	using Base::m_cholmod;

	public:

	typedef _MatrixType MatrixType;

	CholmodSimplicialLLT() : Base() { init(); }

	CholmodSimplicialLLT(const MatrixType& matrix) : Base()
	{
	init();
	this->compute(matrix);
	}

	~CholmodSimplicialLLT() {}
	protected:
	void init()
	{
	m_cholmod.final_asis = 0;
	m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
	m_cholmod.final_ll = 1;
	}
	};


	/** \ingroup CholmodSupport_Module
	* \class CholmodSimplicialLDLT
	* \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod
	*
	* This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization
	* using the Cholmod library.
	* This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical interest.
	* The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
	* X and B can be either dense or sparse.
	*
	* \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.
	*
	* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
	*
	* \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT
	*/
	template<typename _MatrixType, int _UpLo = Lower>
	class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
	{
	typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
	using Base::m_cholmod;

	public:

	typedef _MatrixType MatrixType;

	CholmodSimplicialLDLT() : Base() { init(); }

	CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
	{
	init();
	this->compute(matrix);
	}

	~CholmodSimplicialLDLT() {}
	protected:
	void init()
	{
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
	}
	};

	/** \ingroup CholmodSupport_Module
	* \class CholmodSupernodalLLT
	* \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod
	*
	* This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization
	* using the Cholmod library.
	* This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM.
	* The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
	* X and B can be either dense or sparse.
	*
	* \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.
	*
	* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
	*
	* \sa \ref TutorialSparseDirectSolvers
	*/
	template<typename _MatrixType, int _UpLo = Lower>
	class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
	{
	typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
	using Base::m_cholmod;

	public:

	typedef _MatrixType MatrixType;

	CholmodSupernodalLLT() : Base() { init(); }

	CholmodSupernodalLLT(const MatrixType& matrix) : Base()
	{
	init();
	this->compute(matrix);
	}

	~CholmodSupernodalLLT() {}
	protected:
	void init()
	{
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
	}
	};

	/** \ingroup CholmodSupport_Module
	* \class CholmodDecomposition
	* \brief A general Cholesky factorization and solver based on Cholmod
	*
	* This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
	* using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
	* X and B can be either dense or sparse.
	*
	* This variant permits to change the underlying Cholesky method at runtime.
	* On the other hand, it does not provide access to the result of the factorization.
	* The default is to let Cholmod automatically choose between a simplicial and supernodal factorization.
	*
	* \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.
	*
	* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
	*
	* \sa \ref TutorialSparseDirectSolvers
	*/
	template<typename _MatrixType, int _UpLo = Lower>
	class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
	{
	typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
	using Base::m_cholmod;

	public:

	typedef _MatrixType MatrixType;

	CholmodDecomposition() : Base() { init(); }

	CholmodDecomposition(const MatrixType& matrix) : Base()
	{
	init();
	this->compute(matrix);
	}

	~CholmodDecomposition() {}

	void setMode(CholmodMode mode)
	{
	switch(mode)
	{
	case CholmodAuto:
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_AUTO;
	break;
	case CholmodSimplicialLLt:
	m_cholmod.final_asis = 0;
	m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
	m_cholmod.final_ll = 1;
	break;
	case CholmodSupernodalLLt:
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
	break;
	case CholmodLDLt:
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
	break;
	default:
	break;
	}
	}
	protected:
	void init()
	{
	m_cholmod.final_asis = 1;
	m_cholmod.supernodal = CHOLMOD_AUTO;
	}
	};

	} // end namespace Eigen

	#endif // EIGEN_CHOLMODSUPPORT_H