Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/Eigen/src/PaStiXSupport/PaStiXSupport.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_PASTIXSUPPORT_H
0011 #define EIGEN_PASTIXSUPPORT_H
0012 
0013 namespace Eigen { 
0014 
0015 #if defined(DCOMPLEX)
0016   #define PASTIX_COMPLEX  COMPLEX
0017   #define PASTIX_DCOMPLEX DCOMPLEX
0018 #else
0019   #define PASTIX_COMPLEX  std::complex<float>
0020   #define PASTIX_DCOMPLEX std::complex<double>
0021 #endif
0022 
0023 /** \ingroup PaStiXSupport_Module
0024   * \brief Interface to the PaStix solver
0025   * 
0026   * This class is used to solve the linear systems A.X = B via the PaStix library. 
0027   * The matrix can be either real or complex, symmetric or not.
0028   *
0029   * \sa TutorialSparseDirectSolvers
0030   */
0031 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
0032 template<typename _MatrixType, int Options> class PastixLLT;
0033 template<typename _MatrixType, int Options> class PastixLDLT;
0034 
0035 namespace internal
0036 {
0037     
0038   template<class Pastix> struct pastix_traits;
0039 
0040   template<typename _MatrixType>
0041   struct pastix_traits< PastixLU<_MatrixType> >
0042   {
0043     typedef _MatrixType MatrixType;
0044     typedef typename _MatrixType::Scalar Scalar;
0045     typedef typename _MatrixType::RealScalar RealScalar;
0046     typedef typename _MatrixType::StorageIndex StorageIndex;
0047   };
0048 
0049   template<typename _MatrixType, int Options>
0050   struct pastix_traits< PastixLLT<_MatrixType,Options> >
0051   {
0052     typedef _MatrixType MatrixType;
0053     typedef typename _MatrixType::Scalar Scalar;
0054     typedef typename _MatrixType::RealScalar RealScalar;
0055     typedef typename _MatrixType::StorageIndex StorageIndex;
0056   };
0057 
0058   template<typename _MatrixType, int Options>
0059   struct pastix_traits< PastixLDLT<_MatrixType,Options> >
0060   {
0061     typedef _MatrixType MatrixType;
0062     typedef typename _MatrixType::Scalar Scalar;
0063     typedef typename _MatrixType::RealScalar RealScalar;
0064     typedef typename _MatrixType::StorageIndex StorageIndex;
0065   };
0066   
0067   inline 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)
0068   {
0069     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
0070     if (nbrhs == 0) {x = NULL; nbrhs=1;}
0071     s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
0072   }
0073   
0074   inline 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)
0075   {
0076     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
0077     if (nbrhs == 0) {x = NULL; nbrhs=1;}
0078     d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
0079   }
0080   
0081   inline 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)
0082   {
0083     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
0084     if (nbrhs == 0) {x = NULL; nbrhs=1;}
0085     c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); 
0086   }
0087   
0088   inline 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)
0089   {
0090     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
0091     if (nbrhs == 0) {x = NULL; nbrhs=1;}
0092     z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); 
0093   }
0094 
0095   // Convert the matrix  to Fortran-style Numbering
0096   template <typename MatrixType>
0097   void c_to_fortran_numbering (MatrixType& mat)
0098   {
0099     if ( !(mat.outerIndexPtr()[0]) ) 
0100     { 
0101       int i;
0102       for(i = 0; i <= mat.rows(); ++i)
0103         ++mat.outerIndexPtr()[i];
0104       for(i = 0; i < mat.nonZeros(); ++i)
0105         ++mat.innerIndexPtr()[i];
0106     }
0107   }
0108   
0109   // Convert to C-style Numbering
0110   template <typename MatrixType>
0111   void fortran_to_c_numbering (MatrixType& mat)
0112   {
0113     // Check the Numbering
0114     if ( mat.outerIndexPtr()[0] == 1 ) 
0115     { // Convert to C-style numbering
0116       int i;
0117       for(i = 0; i <= mat.rows(); ++i)
0118         --mat.outerIndexPtr()[i];
0119       for(i = 0; i < mat.nonZeros(); ++i)
0120         --mat.innerIndexPtr()[i];
0121     }
0122   }
0123 }
0124 
0125 // This is the base class to interface with PaStiX functions. 
0126 // Users should not used this class directly. 
0127 template <class Derived>
0128 class PastixBase : public SparseSolverBase<Derived>
0129 {
0130   protected:
0131     typedef SparseSolverBase<Derived> Base;
0132     using Base::derived;
0133     using Base::m_isInitialized;
0134   public:
0135     using Base::_solve_impl;
0136     
0137     typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
0138     typedef _MatrixType MatrixType;
0139     typedef typename MatrixType::Scalar Scalar;
0140     typedef typename MatrixType::RealScalar RealScalar;
0141     typedef typename MatrixType::StorageIndex StorageIndex;
0142     typedef Matrix<Scalar,Dynamic,1> Vector;
0143     typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
0144     enum {
0145       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0146       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0147     };
0148     
0149   public:
0150     
0151     PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0)
0152     {
0153       init();
0154     }
0155     
0156     ~PastixBase() 
0157     {
0158       clean();
0159     }
0160     
0161     template<typename Rhs,typename Dest>
0162     bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
0163     
0164     /** Returns a reference to the integer vector IPARM of PaStiX parameters
0165       * to modify the default parameters. 
0166       * The statistics related to the different phases of factorization and solve are saved here as well
0167       * \sa analyzePattern() factorize()
0168       */
0169     Array<StorageIndex,IPARM_SIZE,1>& iparm()
0170     {
0171       return m_iparm; 
0172     }
0173     
0174     /** Return a reference to a particular index parameter of the IPARM vector 
0175      * \sa iparm()
0176      */
0177     
0178     int& iparm(int idxparam)
0179     {
0180       return m_iparm(idxparam);
0181     }
0182     
0183      /** Returns a reference to the double vector DPARM of PaStiX parameters 
0184       * The statistics related to the different phases of factorization and solve are saved here as well
0185       * \sa analyzePattern() factorize()
0186       */
0187     Array<double,DPARM_SIZE,1>& dparm()
0188     {
0189       return m_dparm; 
0190     }
0191     
0192     
0193     /** Return a reference to a particular index parameter of the DPARM vector 
0194      * \sa dparm()
0195      */
0196     double& dparm(int idxparam)
0197     {
0198       return m_dparm(idxparam);
0199     }
0200     
0201     inline Index cols() const { return m_size; }
0202     inline Index rows() const { return m_size; }
0203     
0204      /** \brief Reports whether previous computation was successful.
0205       *
0206       * \returns \c Success if computation was successful,
0207       *          \c NumericalIssue if the PaStiX reports a problem
0208       *          \c InvalidInput if the input matrix is invalid
0209       *
0210       * \sa iparm()          
0211       */
0212     ComputationInfo info() const
0213     {
0214       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0215       return m_info;
0216     }
0217     
0218   protected:
0219 
0220     // Initialize the Pastix data structure, check the matrix
0221     void init(); 
0222     
0223     // Compute the ordering and the symbolic factorization
0224     void analyzePattern(ColSpMatrix& mat);
0225     
0226     // Compute the numerical factorization
0227     void factorize(ColSpMatrix& mat);
0228     
0229     // Free all the data allocated by Pastix
0230     void clean()
0231     {
0232       eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 
0233       m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
0234       m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
0235       internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
0236                              m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
0237     }
0238     
0239     void compute(ColSpMatrix& mat);
0240     
0241     int m_initisOk; 
0242     int m_analysisIsOk;
0243     int m_factorizationIsOk;
0244     mutable ComputationInfo m_info; 
0245     mutable pastix_data_t *m_pastixdata; // Data structure for pastix
0246     mutable int m_comm; // The MPI communicator identifier
0247     mutable Array<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
0248     mutable Array<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
0249     mutable Matrix<StorageIndex,Dynamic,1> m_perm;  // Permutation vector
0250     mutable Matrix<StorageIndex,Dynamic,1> m_invp;  // Inverse permutation vector
0251     mutable int m_size; // Size of the matrix 
0252 }; 
0253 
0254  /** Initialize the PaStiX data structure. 
0255    *A first call to this function fills iparm and dparm with the default PaStiX parameters
0256    * \sa iparm() dparm()
0257    */
0258 template <class Derived>
0259 void PastixBase<Derived>::init()
0260 {
0261   m_size = 0; 
0262   m_iparm.setZero(IPARM_SIZE);
0263   m_dparm.setZero(DPARM_SIZE);
0264   
0265   m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
0266   pastix(&m_pastixdata, MPI_COMM_WORLD,
0267          0, 0, 0, 0,
0268          0, 0, 0, 1, m_iparm.data(), m_dparm.data());
0269   
0270   m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
0271   m_iparm[IPARM_VERBOSE]             = API_VERBOSE_NOT;
0272   m_iparm[IPARM_ORDERING]            = API_ORDER_SCOTCH;
0273   m_iparm[IPARM_INCOMPLETE]          = API_NO;
0274   m_iparm[IPARM_OOC_LIMIT]           = 2000;
0275   m_iparm[IPARM_RHS_MAKING]          = API_RHS_B;
0276   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
0277   
0278   m_iparm(IPARM_START_TASK) = API_TASK_INIT;
0279   m_iparm(IPARM_END_TASK) = API_TASK_INIT;
0280   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
0281                          0, 0, 0, 0, m_iparm.data(), m_dparm.data());
0282   
0283   // Check the returned error
0284   if(m_iparm(IPARM_ERROR_NUMBER)) {
0285     m_info = InvalidInput;
0286     m_initisOk = false;
0287   }
0288   else { 
0289     m_info = Success;
0290     m_initisOk = true;
0291   }
0292 }
0293 
0294 template <class Derived>
0295 void PastixBase<Derived>::compute(ColSpMatrix& mat)
0296 {
0297   eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
0298   
0299   analyzePattern(mat);  
0300   factorize(mat);
0301   
0302   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
0303 }
0304 
0305 
0306 template <class Derived>
0307 void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
0308 {                         
0309   eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
0310   
0311   // clean previous calls
0312   if(m_size>0)
0313     clean();
0314   
0315   m_size = internal::convert_index<int>(mat.rows());
0316   m_perm.resize(m_size);
0317   m_invp.resize(m_size);
0318   
0319   m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
0320   m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
0321   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
0322                mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
0323   
0324   // Check the returned error
0325   if(m_iparm(IPARM_ERROR_NUMBER))
0326   {
0327     m_info = NumericalIssue;
0328     m_analysisIsOk = false;
0329   }
0330   else
0331   { 
0332     m_info = Success;
0333     m_analysisIsOk = true;
0334   }
0335 }
0336 
0337 template <class Derived>
0338 void PastixBase<Derived>::factorize(ColSpMatrix& mat)
0339 {
0340 //   if(&m_cpyMat != &mat) m_cpyMat = mat;
0341   eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
0342   m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
0343   m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
0344   m_size = internal::convert_index<int>(mat.rows());
0345   
0346   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
0347                mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
0348   
0349   // Check the returned error
0350   if(m_iparm(IPARM_ERROR_NUMBER))
0351   {
0352     m_info = NumericalIssue;
0353     m_factorizationIsOk = false;
0354     m_isInitialized = false;
0355   }
0356   else
0357   {
0358     m_info = Success;
0359     m_factorizationIsOk = true;
0360     m_isInitialized = true;
0361   }
0362 }
0363 
0364 /* Solve the system */
0365 template<typename Base>
0366 template<typename Rhs,typename Dest>
0367 bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
0368 {
0369   eigen_assert(m_isInitialized && "The matrix should be factorized first");
0370   EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
0371                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
0372   int rhs = 1;
0373   
0374   x = b; /* on return, x is overwritten by the computed solution */
0375   
0376   for (int i = 0; i < b.cols(); i++){
0377     m_iparm[IPARM_START_TASK]          = API_TASK_SOLVE;
0378     m_iparm[IPARM_END_TASK]            = API_TASK_REFINE;
0379   
0380     internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, internal::convert_index<int>(x.rows()), 0, 0, 0,
0381                            m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
0382   }
0383   
0384   // Check the returned error
0385   m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
0386   
0387   return m_iparm(IPARM_ERROR_NUMBER)==0;
0388 }
0389 
0390 /** \ingroup PaStiXSupport_Module
0391   * \class PastixLU
0392   * \brief Sparse direct LU solver based on PaStiX library
0393   * 
0394   * This class is used to solve the linear systems A.X = B with a supernodal LU 
0395   * factorization in the PaStiX library. The matrix A should be squared and nonsingular
0396   * PaStiX requires that the matrix A has a symmetric structural pattern. 
0397   * This interface can symmetrize the input matrix otherwise. 
0398   * The vectors or matrices X and B can be either dense or sparse.
0399   * 
0400   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0401   * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false
0402   * NOTE : Note that if the analysis and factorization phase are called separately, 
0403   * the input matrix will be symmetrized at each call, hence it is advised to 
0404   * symmetrize the matrix in a end-user program and set \p IsStrSym to true
0405   *
0406   * \implsparsesolverconcept
0407   *
0408   * \sa \ref TutorialSparseSolverConcept, class SparseLU
0409   * 
0410   */
0411 template<typename _MatrixType, bool IsStrSym>
0412 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
0413 {
0414   public:
0415     typedef _MatrixType MatrixType;
0416     typedef PastixBase<PastixLU<MatrixType> > Base;
0417     typedef typename Base::ColSpMatrix ColSpMatrix;
0418     typedef typename MatrixType::StorageIndex StorageIndex;
0419     
0420   public:
0421     PastixLU() : Base()
0422     {
0423       init();
0424     }
0425     
0426     explicit PastixLU(const MatrixType& matrix):Base()
0427     {
0428       init();
0429       compute(matrix);
0430     }
0431     /** Compute the LU supernodal factorization of \p matrix. 
0432       * iparm and dparm can be used to tune the PaStiX parameters. 
0433       * see the PaStiX user's manual
0434       * \sa analyzePattern() factorize()
0435       */
0436     void compute (const MatrixType& matrix)
0437     {
0438       m_structureIsUptodate = false;
0439       ColSpMatrix temp;
0440       grabMatrix(matrix, temp);
0441       Base::compute(temp);
0442     }
0443     /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern. 
0444       * Several ordering methods can be used at this step. See the PaStiX user's manual. 
0445       * The result of this operation can be used with successive matrices having the same pattern as \p matrix
0446       * \sa factorize()
0447       */
0448     void analyzePattern(const MatrixType& matrix)
0449     {
0450       m_structureIsUptodate = false;
0451       ColSpMatrix temp;
0452       grabMatrix(matrix, temp);
0453       Base::analyzePattern(temp);
0454     }
0455 
0456     /** Compute the LU supernodal factorization of \p matrix
0457       * WARNING The matrix \p matrix should have the same structural pattern 
0458       * as the same used in the analysis phase.
0459       * \sa analyzePattern()
0460       */ 
0461     void factorize(const MatrixType& matrix)
0462     {
0463       ColSpMatrix temp;
0464       grabMatrix(matrix, temp);
0465       Base::factorize(temp);
0466     }
0467   protected:
0468     
0469     void init()
0470     {
0471       m_structureIsUptodate = false;
0472       m_iparm(IPARM_SYM) = API_SYM_NO;
0473       m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
0474     }
0475     
0476     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
0477     {
0478       if(IsStrSym)
0479         out = matrix;
0480       else
0481       {
0482         if(!m_structureIsUptodate)
0483         {
0484           // update the transposed structure
0485           m_transposedStructure = matrix.transpose();
0486           
0487           // Set the elements of the matrix to zero 
0488           for (Index j=0; j<m_transposedStructure.outerSize(); ++j) 
0489             for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
0490               it.valueRef() = 0.0;
0491 
0492           m_structureIsUptodate = true;
0493         }
0494         
0495         out = m_transposedStructure + matrix;
0496       }
0497       internal::c_to_fortran_numbering(out);
0498     }
0499     
0500     using Base::m_iparm;
0501     using Base::m_dparm;
0502     
0503     ColSpMatrix m_transposedStructure;
0504     bool m_structureIsUptodate;
0505 };
0506 
0507 /** \ingroup PaStiXSupport_Module
0508   * \class PastixLLT
0509   * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library
0510   * 
0511   * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization
0512   * available in the PaStiX library. The matrix A should be symmetric and positive definite
0513   * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX
0514   * The vectors or matrices X and B can be either dense or sparse
0515   * 
0516   * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0517   * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX
0518   *
0519   * \implsparsesolverconcept
0520   *
0521   * \sa \ref TutorialSparseSolverConcept, class SimplicialLLT
0522   */
0523 template<typename _MatrixType, int _UpLo>
0524 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
0525 {
0526   public:
0527     typedef _MatrixType MatrixType;
0528     typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
0529     typedef typename Base::ColSpMatrix ColSpMatrix;
0530     
0531   public:
0532     enum { UpLo = _UpLo };
0533     PastixLLT() : Base()
0534     {
0535       init();
0536     }
0537     
0538     explicit PastixLLT(const MatrixType& matrix):Base()
0539     {
0540       init();
0541       compute(matrix);
0542     }
0543 
0544     /** Compute the L factor of the LL^T supernodal factorization of \p matrix 
0545       * \sa analyzePattern() factorize()
0546       */
0547     void compute (const MatrixType& matrix)
0548     {
0549       ColSpMatrix temp;
0550       grabMatrix(matrix, temp);
0551       Base::compute(temp);
0552     }
0553 
0554      /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern
0555       * The result of this operation can be used with successive matrices having the same pattern as \p matrix
0556       * \sa factorize()
0557       */
0558     void analyzePattern(const MatrixType& matrix)
0559     {
0560       ColSpMatrix temp;
0561       grabMatrix(matrix, temp);
0562       Base::analyzePattern(temp);
0563     }
0564       /** Compute the LL^T supernodal numerical factorization of \p matrix 
0565         * \sa analyzePattern()
0566         */
0567     void factorize(const MatrixType& matrix)
0568     {
0569       ColSpMatrix temp;
0570       grabMatrix(matrix, temp);
0571       Base::factorize(temp);
0572     }
0573   protected:
0574     using Base::m_iparm;
0575     
0576     void init()
0577     {
0578       m_iparm(IPARM_SYM) = API_SYM_YES;
0579       m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
0580     }
0581     
0582     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
0583     {
0584       out.resize(matrix.rows(), matrix.cols());
0585       // Pastix supports only lower, column-major matrices 
0586       out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
0587       internal::c_to_fortran_numbering(out);
0588     }
0589 };
0590 
0591 /** \ingroup PaStiXSupport_Module
0592   * \class PastixLDLT
0593   * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library
0594   * 
0595   * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization
0596   * available in the PaStiX library. The matrix A should be symmetric and positive definite
0597   * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX
0598   * The vectors or matrices X and B can be either dense or sparse
0599   * 
0600   * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
0601   * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX
0602   *
0603   * \implsparsesolverconcept
0604   *
0605   * \sa \ref TutorialSparseSolverConcept, class SimplicialLDLT
0606   */
0607 template<typename _MatrixType, int _UpLo>
0608 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
0609 {
0610   public:
0611     typedef _MatrixType MatrixType;
0612     typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 
0613     typedef typename Base::ColSpMatrix ColSpMatrix;
0614     
0615   public:
0616     enum { UpLo = _UpLo };
0617     PastixLDLT():Base()
0618     {
0619       init();
0620     }
0621     
0622     explicit PastixLDLT(const MatrixType& matrix):Base()
0623     {
0624       init();
0625       compute(matrix);
0626     }
0627 
0628     /** Compute the L and D factors of the LDL^T factorization of \p matrix 
0629       * \sa analyzePattern() factorize()
0630       */
0631     void compute (const MatrixType& matrix)
0632     {
0633       ColSpMatrix temp;
0634       grabMatrix(matrix, temp);
0635       Base::compute(temp);
0636     }
0637 
0638     /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern
0639       * The result of this operation can be used with successive matrices having the same pattern as \p matrix
0640       * \sa factorize()
0641       */
0642     void analyzePattern(const MatrixType& matrix)
0643     { 
0644       ColSpMatrix temp;
0645       grabMatrix(matrix, temp);
0646       Base::analyzePattern(temp);
0647     }
0648     /** Compute the LDL^T supernodal numerical factorization of \p matrix 
0649       * 
0650       */
0651     void factorize(const MatrixType& matrix)
0652     {
0653       ColSpMatrix temp;
0654       grabMatrix(matrix, temp);
0655       Base::factorize(temp);
0656     }
0657 
0658   protected:
0659     using Base::m_iparm;
0660     
0661     void init()
0662     {
0663       m_iparm(IPARM_SYM) = API_SYM_YES;
0664       m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
0665     }
0666     
0667     void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
0668     {
0669       // Pastix supports only lower, column-major matrices 
0670       out.resize(matrix.rows(), matrix.cols());
0671       out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
0672       internal::c_to_fortran_numbering(out);
0673     }
0674 };
0675 
0676 } // end namespace Eigen
0677 
0678 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!