Eigen  3.3.3
PaStiXSupport.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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_PASTIXSUPPORT_H
00011 #define EIGEN_PASTIXSUPPORT_H
00012 
00013 namespace Eigen { 
00014 
00015 #if defined(DCOMPLEX)
00016   #define PASTIX_COMPLEX  COMPLEX
00017   #define PASTIX_DCOMPLEX DCOMPLEX
00018 #else
00019   #define PASTIX_COMPLEX  std::complex<float>
00020   #define PASTIX_DCOMPLEX std::complex<double>
00021 #endif
00022 
00031 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
00032 template<typename _MatrixType, int Options> class PastixLLT;
00033 template<typename _MatrixType, int Options> class PastixLDLT;
00034 
00035 namespace internal
00036 {
00037     
00038   template<class Pastix> struct pastix_traits;
00039 
00040   template<typename _MatrixType>
00041   struct pastix_traits< PastixLU<_MatrixType> >
00042   {
00043     typedef _MatrixType MatrixType;
00044     typedef typename _MatrixType::Scalar Scalar;
00045     typedef typename _MatrixType::RealScalar RealScalar;
00046     typedef typename _MatrixType::StorageIndex StorageIndex;
00047   };
00048 
00049   template<typename _MatrixType, int Options>
00050   struct pastix_traits< PastixLLT<_MatrixType,Options> >
00051   {
00052     typedef _MatrixType MatrixType;
00053     typedef typename _MatrixType::Scalar Scalar;
00054     typedef typename _MatrixType::RealScalar RealScalar;
00055     typedef typename _MatrixType::StorageIndex StorageIndex;
00056   };
00057 
00058   template<typename _MatrixType, int Options>
00059   struct pastix_traits< PastixLDLT<_MatrixType,Options> >
00060   {
00061     typedef _MatrixType MatrixType;
00062     typedef typename _MatrixType::Scalar Scalar;
00063     typedef typename _MatrixType::RealScalar RealScalar;
00064     typedef typename _MatrixType::StorageIndex StorageIndex;
00065   };
00066   
00067   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm)
00068   {
00069     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00070     if (nbrhs == 0) {x = NULL; nbrhs=1;}
00071     s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
00072   }
00073   
00074   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm)
00075   {
00076     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00077     if (nbrhs == 0) {x = NULL; nbrhs=1;}
00078     d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
00079   }
00080   
00081   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm)
00082   {
00083     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00084     if (nbrhs == 0) {x = NULL; nbrhs=1;}
00085     c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); 
00086   }
00087   
00088   void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm)
00089   {
00090     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00091     if (nbrhs == 0) {x = NULL; nbrhs=1;}
00092     z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); 
00093   }
00094 
00095   // Convert the matrix  to Fortran-style Numbering
00096   template <typename MatrixType>
00097   void c_to_fortran_numbering (MatrixType& mat)
00098   {
00099     if ( !(mat.outerIndexPtr()[0]) ) 
00100     { 
00101       int i;
00102       for(i = 0; i <= mat.rows(); ++i)
00103         ++mat.outerIndexPtr()[i];
00104       for(i = 0; i < mat.nonZeros(); ++i)
00105         ++mat.innerIndexPtr()[i];
00106     }
00107   }
00108   
00109   // Convert to C-style Numbering
00110   template <typename MatrixType>
00111   void fortran_to_c_numbering (MatrixType& mat)
00112   {
00113     // Check the Numbering
00114     if ( mat.outerIndexPtr()[0] == 1 ) 
00115     { // Convert to C-style numbering
00116       int i;
00117       for(i = 0; i <= mat.rows(); ++i)
00118         --mat.outerIndexPtr()[i];
00119       for(i = 0; i < mat.nonZeros(); ++i)
00120         --mat.innerIndexPtr()[i];
00121     }
00122   }
00123 }
00124 
00125 // This is the base class to interface with PaStiX functions. 
00126 // Users should not used this class directly. 
00127 template <class Derived>
00128 class PastixBase : public SparseSolverBase<Derived>
00129 {
00130   protected:
00131     typedef SparseSolverBase<Derived> Base;
00132     using Base::derived;
00133     using Base::m_isInitialized;
00134   public:
00135     using Base::_solve_impl;
00136     
00137     typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
00138     typedef _MatrixType MatrixType;
00139     typedef typename MatrixType::Scalar Scalar;
00140     typedef typename MatrixType::RealScalar RealScalar;
00141     typedef typename MatrixType::StorageIndex StorageIndex;
00142     typedef Matrix<Scalar,Dynamic,1> Vector;
00143     typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
00144     enum {
00145       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00146       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00147     };
00148     
00149   public:
00150     
00151     PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
00152     {
00153       init();
00154     }
00155     
00156     ~PastixBase() 
00157     {
00158       clean();
00159     }
00160     
00161     template<typename Rhs,typename Dest>
00162     bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
00163     
00169     Array<StorageIndex,IPARM_SIZE,1>& iparm()
00170     {
00171       return m_iparm; 
00172     }
00173     
00178     int& iparm(int idxparam)
00179     {
00180       return m_iparm(idxparam);
00181     }
00182     
00187     Array<double,DPARM_SIZE,1>& dparm()
00188     {
00189       return m_dparm; 
00190     }
00191     
00192     
00196     double& dparm(int idxparam)
00197     {
00198       return m_dparm(idxparam);
00199     }
00200     
00201     inline Index cols() const { return m_size; }
00202     inline Index rows() const { return m_size; }
00203     
00212     ComputationInfo info() const
00213     {
00214       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00215       return m_info;
00216     }
00217     
00218   protected:
00219 
00220     // Initialize the Pastix data structure, check the matrix
00221     void init(); 
00222     
00223     // Compute the ordering and the symbolic factorization
00224     void analyzePattern(ColSpMatrix& mat);
00225     
00226     // Compute the numerical factorization
00227     void factorize(ColSpMatrix& mat);
00228     
00229     // Free all the data allocated by Pastix
00230     void clean()
00231     {
00232       eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 
00233       m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
00234       m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
00235       internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
00236                              m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
00237     }
00238     
00239     void compute(ColSpMatrix& mat);
00240     
00241     int m_initisOk; 
00242     int m_analysisIsOk;
00243     int m_factorizationIsOk;
00244     mutable ComputationInfo m_info; 
00245     mutable pastix_data_t *m_pastixdata; // Data structure for pastix
00246     mutable int m_comm; // The MPI communicator identifier
00247     mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
00248     mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
00249     mutable Matrix<StorageIndex,Dynamic,1> m_perm;  // Permutation vector
00250     mutable Matrix<StorageIndex,Dynamic,1> m_invp;  // Inverse permutation vector
00251     mutable int m_size; // Size of the matrix 
00252 }; 
00253 
00258 template <class Derived>
00259 void PastixBase<Derived>::init()
00260 {
00261   m_size = 0; 
00262   m_iparm.setZero(IPARM_SIZE);
00263   m_dparm.setZero(DPARM_SIZE);
00264   
00265   m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
00266   pastix(&m_pastixdata, MPI_COMM_WORLD,
00267          0, 0, 0, 0,
00268          0, 0, 0, 1, m_iparm.data(), m_dparm.data());
00269   
00270   m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
00271   m_iparm[IPARM_VERBOSE]             = API_VERBOSE_NOT;
00272   m_iparm[IPARM_ORDERING]            = API_ORDER_SCOTCH;
00273   m_iparm[IPARM_INCOMPLETE]          = API_NO;
00274   m_iparm[IPARM_OOC_LIMIT]           = 2000;
00275   m_iparm[IPARM_RHS_MAKING]          = API_RHS_B;
00276   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00277   
00278   m_iparm(IPARM_START_TASK) = API_TASK_INIT;
00279   m_iparm(IPARM_END_TASK) = API_TASK_INIT;
00280   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
00281                          0, 0, 0, 0, m_iparm.data(), m_dparm.data());
00282   
00283   // Check the returned error
00284   if(m_iparm(IPARM_ERROR_NUMBER)) {
00285     m_info = InvalidInput;
00286     m_initisOk = false;
00287   }
00288   else { 
00289     m_info = Success;
00290     m_initisOk = true;
00291   }
00292 }
00293 
00294 template <class Derived>
00295 void PastixBase<Derived>::compute(ColSpMatrix& mat)
00296 {
00297   eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
00298   
00299   analyzePattern(mat);  
00300   factorize(mat);
00301   
00302   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00303 }
00304 
00305 
00306 template <class Derived>
00307 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
00308 {                         
00309   eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
00310   
00311   // clean previous calls
00312   if(m_size>0)
00313     clean();
00314   
00315   m_size = internal::convert_index<int>(mat.rows());
00316   m_perm.resize(m_size);
00317   m_invp.resize(m_size);
00318   
00319   m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
00320   m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
00321   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00322                mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
00323   
00324   // Check the returned error
00325   if(m_iparm(IPARM_ERROR_NUMBER))
00326   {
00327     m_info = NumericalIssue;
00328     m_analysisIsOk = false;
00329   }
00330   else
00331   { 
00332     m_info = Success;
00333     m_analysisIsOk = true;
00334   }
00335 }
00336 
00337 template <class Derived>
00338 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
00339 {
00340 //   if(&m_cpyMat != &mat) m_cpyMat = mat;
00341   eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
00342   m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
00343   m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
00344   m_size = internal::convert_index<int>(mat.rows());
00345   
00346   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00347                mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
00348   
00349   // Check the returned error
00350   if(m_iparm(IPARM_ERROR_NUMBER))
00351   {
00352     m_info = NumericalIssue;
00353     m_factorizationIsOk = false;
00354     m_isInitialized = false;
00355   }
00356   else
00357   {
00358     m_info = Success;
00359     m_factorizationIsOk = true;
00360     m_isInitialized = true;
00361   }
00362 }
00363 
00364 /* Solve the system */
00365 template<typename Base>
00366 template<typename Rhs,typename Dest>
00367 bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
00368 {
00369   eigen_assert(m_isInitialized && "The matrix should be factorized first");
00370   EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
00371                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00372   int rhs = 1;
00373   
00374   x = b; /* on return, x is overwritten by the computed solution */
00375   
00376   for (int i = 0; i < b.cols(); i++){
00377     m_iparm[IPARM_START_TASK]          = API_TASK_SOLVE;
00378     m_iparm[IPARM_END_TASK]            = API_TASK_REFINE;
00379   
00380     internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0,
00381                            m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
00382   }
00383   
00384   // Check the returned error
00385   m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
00386   
00387   return m_iparm(IPARM_ERROR_NUMBER)==0;
00388 }
00389 
00411 template<typename _MatrixType, bool IsStrSym>
00412 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
00413 {
00414   public:
00415     typedef _MatrixType MatrixType;
00416     typedef PastixBase<PastixLU<MatrixType> > Base;
00417     typedef typename Base::ColSpMatrix ColSpMatrix;
00418     typedef typename MatrixType::StorageIndex StorageIndex;
00419     
00420   public:
00421     PastixLU() : Base()
00422     {
00423       init();
00424     }
00425     
00426     explicit PastixLU(const MatrixType& matrix):Base()
00427     {
00428       init();
00429       compute(matrix);
00430     }
00436     void compute (const MatrixType& matrix)
00437     {
00438       m_structureIsUptodate = false;
00439       ColSpMatrix temp;
00440       grabMatrix(matrix, temp);
00441       Base::compute(temp);
00442     }
00448     void analyzePattern(const MatrixType& matrix)
00449     {
00450       m_structureIsUptodate = false;
00451       ColSpMatrix temp;
00452       grabMatrix(matrix, temp);
00453       Base::analyzePattern(temp);
00454     }
00455 
00461     void factorize(const MatrixType& matrix)
00462     {
00463       ColSpMatrix temp;
00464       grabMatrix(matrix, temp);
00465       Base::factorize(temp);
00466     }
00467   protected:
00468     
00469     void init()
00470     {
00471       m_structureIsUptodate = false;
00472       m_iparm(IPARM_SYM) = API_SYM_NO;
00473       m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00474     }
00475     
00476     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
00477     {
00478       if(IsStrSym)
00479         out = matrix;
00480       else
00481       {
00482         if(!m_structureIsUptodate)
00483         {
00484           // update the transposed structure
00485           m_transposedStructure = matrix.transpose();
00486           
00487           // Set the elements of the matrix to zero 
00488           for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 
00489             for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
00490               it.valueRef() = 0.0;
00491 
00492           m_structureIsUptodate = true;
00493         }
00494         
00495         out = m_transposedStructure + matrix;
00496       }
00497       internal::c_to_fortran_numbering(out);
00498     }
00499     
00500     using Base::m_iparm;
00501     using Base::m_dparm;
00502     
00503     ColSpMatrix m_transposedStructure;
00504     bool m_structureIsUptodate;
00505 };
00506 
00523 template<typename _MatrixType, int _UpLo>
00524 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
00525 {
00526   public:
00527     typedef _MatrixType MatrixType;
00528     typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
00529     typedef typename Base::ColSpMatrix ColSpMatrix;
00530     
00531   public:
00532     enum { UpLo = _UpLo };
00533     PastixLLT() : Base()
00534     {
00535       init();
00536     }
00537     
00538     explicit PastixLLT(const MatrixType& matrix):Base()
00539     {
00540       init();
00541       compute(matrix);
00542     }
00543 
00547     void compute (const MatrixType& matrix)
00548     {
00549       ColSpMatrix temp;
00550       grabMatrix(matrix, temp);
00551       Base::compute(temp);
00552     }
00553 
00558     void analyzePattern(const MatrixType& matrix)
00559     {
00560       ColSpMatrix temp;
00561       grabMatrix(matrix, temp);
00562       Base::analyzePattern(temp);
00563     }
00567     void factorize(const MatrixType& matrix)
00568     {
00569       ColSpMatrix temp;
00570       grabMatrix(matrix, temp);
00571       Base::factorize(temp);
00572     }
00573   protected:
00574     using Base::m_iparm;
00575     
00576     void init()
00577     {
00578       m_iparm(IPARM_SYM) = API_SYM_YES;
00579       m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00580     }
00581     
00582     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
00583     {
00584       out.resize(matrix.rows(), matrix.cols());
00585       // Pastix supports only lower, column-major matrices 
00586       out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
00587       internal::c_to_fortran_numbering(out);
00588     }
00589 };
00590 
00607 template<typename _MatrixType, int _UpLo>
00608 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
00609 {
00610   public:
00611     typedef _MatrixType MatrixType;
00612     typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 
00613     typedef typename Base::ColSpMatrix ColSpMatrix;
00614     
00615   public:
00616     enum { UpLo = _UpLo };
00617     PastixLDLT():Base()
00618     {
00619       init();
00620     }
00621     
00622     explicit PastixLDLT(const MatrixType& matrix):Base()
00623     {
00624       init();
00625       compute(matrix);
00626     }
00627 
00631     void compute (const MatrixType& matrix)
00632     {
00633       ColSpMatrix temp;
00634       grabMatrix(matrix, temp);
00635       Base::compute(temp);
00636     }
00637 
00642     void analyzePattern(const MatrixType& matrix)
00643     { 
00644       ColSpMatrix temp;
00645       grabMatrix(matrix, temp);
00646       Base::analyzePattern(temp);
00647     }
00651     void factorize(const MatrixType& matrix)
00652     {
00653       ColSpMatrix temp;
00654       grabMatrix(matrix, temp);
00655       Base::factorize(temp);
00656     }
00657 
00658   protected:
00659     using Base::m_iparm;
00660     
00661     void init()
00662     {
00663       m_iparm(IPARM_SYM) = API_SYM_YES;
00664       m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00665     }
00666     
00667     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
00668     {
00669       // Pastix supports only lower, column-major matrices 
00670       out.resize(matrix.rows(), matrix.cols());
00671       out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
00672       internal::c_to_fortran_numbering(out);
00673     }
00674 };
00675 
00676 } // end namespace Eigen
00677 
00678 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends