CholmodSupport.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_CHOLMODSUPPORT_H
00011 #define EIGEN_CHOLMODSUPPORT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Scalar, typename CholmodType>
00018 void cholmod_configure_matrix(CholmodType& mat)
00019 {
00020   if (internal::is_same<Scalar,float>::value)
00021   {
00022     mat.xtype = CHOLMOD_REAL;
00023     mat.dtype = CHOLMOD_SINGLE;
00024   }
00025   else if (internal::is_same<Scalar,double>::value)
00026   {
00027     mat.xtype = CHOLMOD_REAL;
00028     mat.dtype = CHOLMOD_DOUBLE;
00029   }
00030   else if (internal::is_same<Scalar,std::complex<float> >::value)
00031   {
00032     mat.xtype = CHOLMOD_COMPLEX;
00033     mat.dtype = CHOLMOD_SINGLE;
00034   }
00035   else if (internal::is_same<Scalar,std::complex<double> >::value)
00036   {
00037     mat.xtype = CHOLMOD_COMPLEX;
00038     mat.dtype = CHOLMOD_DOUBLE;
00039   }
00040   else
00041   {
00042     eigen_assert(false && "Scalar type not supported by CHOLMOD");
00043   }
00044 }
00045 
00046 } // namespace internal
00047 
00051 template<typename _Scalar, int _Options, typename _Index>
00052 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
00053 {
00054   cholmod_sparse res;
00055   res.nzmax   = mat.nonZeros();
00056   res.nrow    = mat.rows();;
00057   res.ncol    = mat.cols();
00058   res.p       = mat.outerIndexPtr();
00059   res.i       = mat.innerIndexPtr();
00060   res.x       = mat.valuePtr();
00061   res.sorted  = 1;
00062   if(mat.isCompressed())
00063   {
00064     res.packed  = 1;
00065   }
00066   else
00067   {
00068     res.packed  = 0;
00069     res.nz = mat.innerNonZeroPtr();
00070   }
00071 
00072   res.dtype   = 0;
00073   res.stype   = -1;
00074   
00075   if (internal::is_same<_Index,int>::value)
00076   {
00077     res.itype = CHOLMOD_INT;
00078   }
00079   else if (internal::is_same<_Index,UF_long>::value)
00080   {
00081     res.itype = CHOLMOD_LONG;
00082   }
00083   else
00084   {
00085     eigen_assert(false && "Index type not supported yet");
00086   }
00087 
00088   // setup res.xtype
00089   internal::cholmod_configure_matrix<_Scalar>(res);
00090   
00091   res.stype = 0;
00092   
00093   return res;
00094 }
00095 
00096 template<typename _Scalar, int _Options, typename _Index>
00097 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
00098 {
00099   cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
00100   return res;
00101 }
00102 
00105 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
00106 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
00107 {
00108   cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
00109   
00110   if(UpLo==Upper) res.stype =  1;
00111   if(UpLo==Lower) res.stype = -1;
00112 
00113   return res;
00114 }
00115 
00118 template<typename Derived>
00119 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
00120 {
00121   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00122   typedef typename Derived::Scalar Scalar;
00123 
00124   cholmod_dense res;
00125   res.nrow   = mat.rows();
00126   res.ncol   = mat.cols();
00127   res.nzmax  = res.nrow * res.ncol;
00128   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
00129   res.x      = (void*)(mat.derived().data());
00130   res.z      = 0;
00131 
00132   internal::cholmod_configure_matrix<Scalar>(res);
00133 
00134   return res;
00135 }
00136 
00139 template<typename Scalar, int Flags, typename Index>
00140 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
00141 {
00142   return MappedSparseMatrix<Scalar,Flags,Index>
00143          (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol],
00144           static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) );
00145 }
00146 
00147 enum CholmodMode {
00148   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
00149 };
00150 
00151 
00157 template<typename _MatrixType, int _UpLo, typename Derived>
00158 class CholmodBase : internal::noncopyable
00159 {
00160   public:
00161     typedef _MatrixType MatrixType;
00162     enum { UpLo = _UpLo };
00163     typedef typename MatrixType::Scalar Scalar;
00164     typedef typename MatrixType::RealScalar RealScalar;
00165     typedef MatrixType CholMatrixType;
00166     typedef typename MatrixType::Index Index;
00167 
00168   public:
00169 
00170     CholmodBase()
00171       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00172     {
00173       cholmod_start(&m_cholmod);
00174     }
00175 
00176     CholmodBase(const MatrixType& matrix)
00177       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00178     {
00179       m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
00180       cholmod_start(&m_cholmod);
00181       compute(matrix);
00182     }
00183 
00184     ~CholmodBase()
00185     {
00186       if(m_cholmodFactor)
00187         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00188       cholmod_finish(&m_cholmod);
00189     }
00190     
00191     inline Index cols() const { return m_cholmodFactor->n; }
00192     inline Index rows() const { return m_cholmodFactor->n; }
00193     
00194     Derived& derived() { return *static_cast<Derived*>(this); }
00195     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00196     
00202     ComputationInfo info() const
00203     {
00204       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00205       return m_info;
00206     }
00207 
00209     Derived& compute(const MatrixType& matrix)
00210     {
00211       analyzePattern(matrix);
00212       factorize(matrix);
00213       return derived();
00214     }
00215     
00220     template<typename Rhs>
00221     inline const internal::solve_retval<CholmodBase, Rhs>
00222     solve(const MatrixBase<Rhs>& b) const
00223     {
00224       eigen_assert(m_isInitialized && "LLT is not initialized.");
00225       eigen_assert(rows()==b.rows()
00226                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00227       return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
00228     }
00229     
00234     template<typename Rhs>
00235     inline const internal::sparse_solve_retval<CholmodBase, Rhs>
00236     solve(const SparseMatrixBase<Rhs>& b) const
00237     {
00238       eigen_assert(m_isInitialized && "LLT is not initialized.");
00239       eigen_assert(rows()==b.rows()
00240                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00241       return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
00242     }
00243     
00250     void analyzePattern(const MatrixType& matrix)
00251     {
00252       if(m_cholmodFactor)
00253       {
00254         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00255         m_cholmodFactor = 0;
00256       }
00257       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00258       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00259       
00260       this->m_isInitialized = true;
00261       this->m_info = Success;
00262       m_analysisIsOk = true;
00263       m_factorizationIsOk = false;
00264     }
00265     
00272     void factorize(const MatrixType& matrix)
00273     {
00274       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00275       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00276       cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
00277       
00278       // If the factorization failed, minor is the column at which it did. On success minor == n.
00279       this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
00280       m_factorizationIsOk = true;
00281     }
00282     
00285     cholmod_common& cholmod() { return m_cholmod; }
00286     
00287     #ifndef EIGEN_PARSED_BY_DOXYGEN
00288 
00289     template<typename Rhs,typename Dest>
00290     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00291     {
00292       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00293       const Index size = m_cholmodFactor->n;
00294       EIGEN_UNUSED_VARIABLE(size);
00295       eigen_assert(size==b.rows());
00296 
00297       // note: cd stands for Cholmod Dense
00298       Rhs& b_ref(b.const_cast_derived());
00299       cholmod_dense b_cd = viewAsCholmod(b_ref);
00300       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
00301       if(!x_cd)
00302       {
00303         this->m_info = NumericalIssue;
00304       }
00305       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00306       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
00307       cholmod_free_dense(&x_cd, &m_cholmod);
00308     }
00309     
00311     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
00312     void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00313     {
00314       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00315       const Index size = m_cholmodFactor->n;
00316       EIGEN_UNUSED_VARIABLE(size);
00317       eigen_assert(size==b.rows());
00318 
00319       // note: cs stands for Cholmod Sparse
00320       cholmod_sparse b_cs = viewAsCholmod(b);
00321       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
00322       if(!x_cs)
00323       {
00324         this->m_info = NumericalIssue;
00325       }
00326       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00327       dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
00328       cholmod_free_sparse(&x_cs, &m_cholmod);
00329     }
00330     #endif // EIGEN_PARSED_BY_DOXYGEN
00331     
00332     
00342     Derived& setShift(const RealScalar& offset)
00343     {
00344       m_shiftOffset[0] = offset;
00345       return derived();
00346     }
00347     
00348     template<typename Stream>
00349     void dumpMemory(Stream& /*s*/)
00350     {}
00351     
00352   protected:
00353     mutable cholmod_common m_cholmod;
00354     cholmod_factor* m_cholmodFactor;
00355     RealScalar m_shiftOffset[2];
00356     mutable ComputationInfo m_info;
00357     bool m_isInitialized;
00358     int m_factorizationIsOk;
00359     int m_analysisIsOk;
00360 };
00361 
00380 template<typename _MatrixType, int _UpLo = Lower>
00381 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
00382 {
00383     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
00384     using Base::m_cholmod;
00385     
00386   public:
00387     
00388     typedef _MatrixType MatrixType;
00389     
00390     CholmodSimplicialLLT() : Base() { init(); }
00391 
00392     CholmodSimplicialLLT(const MatrixType& matrix) : Base()
00393     {
00394       init();
00395       compute(matrix);
00396     }
00397 
00398     ~CholmodSimplicialLLT() {}
00399   protected:
00400     void init()
00401     {
00402       m_cholmod.final_asis = 0;
00403       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00404       m_cholmod.final_ll = 1;
00405     }
00406 };
00407 
00408 
00427 template<typename _MatrixType, int _UpLo = Lower>
00428 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
00429 {
00430     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
00431     using Base::m_cholmod;
00432     
00433   public:
00434     
00435     typedef _MatrixType MatrixType;
00436     
00437     CholmodSimplicialLDLT() : Base() { init(); }
00438 
00439     CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
00440     {
00441       init();
00442       compute(matrix);
00443     }
00444 
00445     ~CholmodSimplicialLDLT() {}
00446   protected:
00447     void init()
00448     {
00449       m_cholmod.final_asis = 1;
00450       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00451     }
00452 };
00453 
00472 template<typename _MatrixType, int _UpLo = Lower>
00473 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
00474 {
00475     typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
00476     using Base::m_cholmod;
00477     
00478   public:
00479     
00480     typedef _MatrixType MatrixType;
00481     
00482     CholmodSupernodalLLT() : Base() { init(); }
00483 
00484     CholmodSupernodalLLT(const MatrixType& matrix) : Base()
00485     {
00486       init();
00487       compute(matrix);
00488     }
00489 
00490     ~CholmodSupernodalLLT() {}
00491   protected:
00492     void init()
00493     {
00494       m_cholmod.final_asis = 1;
00495       m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00496     }
00497 };
00498 
00519 template<typename _MatrixType, int _UpLo = Lower>
00520 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
00521 {
00522     typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
00523     using Base::m_cholmod;
00524     
00525   public:
00526     
00527     typedef _MatrixType MatrixType;
00528     
00529     CholmodDecomposition() : Base() { init(); }
00530 
00531     CholmodDecomposition(const MatrixType& matrix) : Base()
00532     {
00533       init();
00534       compute(matrix);
00535     }
00536 
00537     ~CholmodDecomposition() {}
00538     
00539     void setMode(CholmodMode mode)
00540     {
00541       switch(mode)
00542       {
00543         case CholmodAuto:
00544           m_cholmod.final_asis = 1;
00545           m_cholmod.supernodal = CHOLMOD_AUTO;
00546           break;
00547         case CholmodSimplicialLLt:
00548           m_cholmod.final_asis = 0;
00549           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00550           m_cholmod.final_ll = 1;
00551           break;
00552         case CholmodSupernodalLLt:
00553           m_cholmod.final_asis = 1;
00554           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00555           break;
00556         case CholmodLDLt:
00557           m_cholmod.final_asis = 1;
00558           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00559           break;
00560         default:
00561           break;
00562       }
00563     }
00564   protected:
00565     void init()
00566     {
00567       m_cholmod.final_asis = 1;
00568       m_cholmod.supernodal = CHOLMOD_AUTO;
00569     }
00570 };
00571 
00572 namespace internal {
00573   
00574 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00575 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00576   : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00577 {
00578   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00579   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00580 
00581   template<typename Dest> void evalTo(Dest& dst) const
00582   {
00583     dec()._solve(rhs(),dst);
00584   }
00585 };
00586 
00587 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00588 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00589   : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00590 {
00591   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00592   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00593 
00594   template<typename Dest> void evalTo(Dest& dst) const
00595   {
00596     dec()._solve(rhs(),dst);
00597   }
00598 };
00599 
00600 } // end namespace internal
00601 
00602 } // end namespace Eigen
00603 
00604 #endif // EIGEN_CHOLMODSUPPORT_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:49