vcglib: CholmodSupport.h Source File

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_CHOLMODSUPPORT_H
00026 #define EIGEN_CHOLMODSUPPORT_H
00027 
00028 template<typename Scalar, typename CholmodType>
00029 void ei_cholmod_configure_matrix(CholmodType& mat)
00030 {
00031   if (ei_is_same_type<Scalar,float>::ret)
00032   {
00033     mat.xtype = CHOLMOD_REAL;
00034     mat.dtype = 1;
00035   }
00036   else if (ei_is_same_type<Scalar,double>::ret)
00037   {
00038     mat.xtype = CHOLMOD_REAL;
00039     mat.dtype = 0;
00040   }
00041   else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
00042   {
00043     mat.xtype = CHOLMOD_COMPLEX;
00044     mat.dtype = 1;
00045   }
00046   else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
00047   {
00048     mat.xtype = CHOLMOD_COMPLEX;
00049     mat.dtype = 0;
00050   }
00051   else
00052   {
00053     ei_assert(false && "Scalar type not supported by CHOLMOD");
00054   }
00055 }
00056 
00057 template<typename Derived>
00058 cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
00059 {
00060   typedef typename Derived::Scalar Scalar;
00061   cholmod_sparse res;
00062   res.nzmax   = nonZeros();
00063   res.nrow    = rows();;
00064   res.ncol    = cols();
00065   res.p       = derived()._outerIndexPtr();
00066   res.i       = derived()._innerIndexPtr();
00067   res.x       = derived()._valuePtr();
00068   res.xtype = CHOLMOD_REAL;
00069   res.itype = CHOLMOD_INT;
00070   res.sorted = 1;
00071   res.packed = 1;
00072   res.dtype = 0;
00073   res.stype = -1;
00074 
00075   ei_cholmod_configure_matrix<Scalar>(res);
00076 
00077   if (Derived::Flags & SelfAdjoint)
00078   {
00079     if (Derived::Flags & UpperTriangular)
00080       res.stype = 1;
00081     else if (Derived::Flags & LowerTriangular)
00082       res.stype = -1;
00083     else
00084       res.stype = 0;
00085   }
00086   else
00087     res.stype = 0;
00088 
00089   return res;
00090 }
00091 
00092 template<typename Derived>
00093 cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
00094 {
00095   EIGEN_STATIC_ASSERT((ei_traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00096   typedef typename Derived::Scalar Scalar;
00097 
00098   cholmod_dense res;
00099   res.nrow   = mat.rows();
00100   res.ncol   = mat.cols();
00101   res.nzmax  = res.nrow * res.ncol;
00102   res.d      = mat.derived().stride();
00103   res.x      = mat.derived().data();
00104   res.z      = 0;
00105 
00106   ei_cholmod_configure_matrix<Scalar>(res);
00107 
00108   return res;
00109 }
00110 
00111 template<typename Scalar, int Flags>
00112 MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(cholmod_sparse& cm)
00113 {
00114   m_innerSize = cm.nrow;
00115   m_outerSize = cm.ncol;
00116   m_outerIndex = reinterpret_cast<int*>(cm.p);
00117   m_innerIndices = reinterpret_cast<int*>(cm.i);
00118   m_values = reinterpret_cast<Scalar*>(cm.x);
00119   m_nnz = m_outerIndex[cm.ncol];
00120 }
00121 
00122 template<typename MatrixType>
00123 class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
00124 {
00125   protected:
00126     typedef SparseLLT<MatrixType> Base;
00127     typedef typename Base::Scalar Scalar;
00128     typedef typename Base::RealScalar RealScalar;
00129     using Base::MatrixLIsDirty;
00130     using Base::SupernodalFactorIsDirty;
00131     using Base::m_flags;
00132     using Base::m_matrix;
00133     using Base::m_status;
00134 
00135   public:
00136 
00137     SparseLLT(int flags = 0)
00138       : Base(flags), m_cholmodFactor(0)
00139     {
00140       cholmod_start(&m_cholmod);
00141     }
00142 
00143     SparseLLT(const MatrixType& matrix, int flags = 0)
00144       : Base(flags), m_cholmodFactor(0)
00145     {
00146       cholmod_start(&m_cholmod);
00147       compute(matrix);
00148     }
00149 
00150     ~SparseLLT()
00151     {
00152       if (m_cholmodFactor)
00153         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00154       cholmod_finish(&m_cholmod);
00155     }
00156 
00157     inline const typename Base::CholMatrixType& matrixL(void) const;
00158 
00159     template<typename Derived>
00160     void solveInPlace(MatrixBase<Derived> &b) const;
00161 
00162     void compute(const MatrixType& matrix);
00163 
00164   protected:
00165     mutable cholmod_common m_cholmod;
00166     cholmod_factor* m_cholmodFactor;
00167 };
00168 
00169 template<typename MatrixType>
00170 void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
00171 {
00172   if (m_cholmodFactor)
00173   {
00174     cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00175     m_cholmodFactor = 0;
00176   }
00177 
00178   cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
00179   m_cholmod.supernodal = CHOLMOD_AUTO;
00180   // TODO
00181   if (m_flags&IncompleteFactorization)
00182   {
00183     m_cholmod.nmethods = 1;
00184     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
00185     m_cholmod.postorder = 0;
00186   }
00187   else
00188   {
00189     m_cholmod.nmethods = 1;
00190     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
00191     m_cholmod.postorder = 0;
00192   }
00193   m_cholmod.final_ll = 1;
00194   m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00195   cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
00196 
00197   m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
00198 }
00199 
00200 template<typename MatrixType>
00201 inline const typename SparseLLT<MatrixType>::CholMatrixType&
00202 SparseLLT<MatrixType,Cholmod>::matrixL() const
00203 {
00204   if (m_status & MatrixLIsDirty)
00205   {
00206     ei_assert(!(m_status & SupernodalFactorIsDirty));
00207 
00208     cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
00209     const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
00210     free(cmRes);
00211 
00212     m_status = (m_status & ~MatrixLIsDirty);
00213   }
00214   return m_matrix;
00215 }
00216 
00217 template<typename MatrixType>
00218 template<typename Derived>
00219 void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
00220 {
00221   const int size = m_cholmodFactor->n;
00222   ei_assert(size==b.rows());
00223 
00224   // this uses Eigen's triangular sparse solver
00225 //   if (m_status & MatrixLIsDirty)
00226 //     matrixL();
00227 //   Base::solveInPlace(b);
00228   // as long as our own triangular sparse solver is not fully optimal,
00229   // let's use CHOLMOD's one:
00230   cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
00231   cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
00232   b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
00233   cholmod_free_dense(&x, &m_cholmod);
00234 }
00235 
00236 #endif // EIGEN_CHOLMODSUPPORT_H