PaStiXSupport.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) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.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_PASTIXSUPPORT_H
00026 #define EIGEN_PASTIXSUPPORT_H
00027 
00028 namespace Eigen { 
00029 
00039 template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
00040 template<typename _MatrixType, int Options> class PastixLLT;
00041 template<typename _MatrixType, int Options> class PastixLDLT;
00042 
00043 namespace internal
00044 {
00045     
00046   template<class Pastix> struct pastix_traits;
00047 
00048   template<typename _MatrixType>
00049   struct pastix_traits< PastixLU<_MatrixType> >
00050   {
00051     typedef _MatrixType MatrixType;
00052     typedef typename _MatrixType::Scalar Scalar;
00053     typedef typename _MatrixType::RealScalar RealScalar;
00054     typedef typename _MatrixType::Index Index;
00055   };
00056 
00057   template<typename _MatrixType, int Options>
00058   struct pastix_traits< PastixLLT<_MatrixType,Options> >
00059   {
00060     typedef _MatrixType MatrixType;
00061     typedef typename _MatrixType::Scalar Scalar;
00062     typedef typename _MatrixType::RealScalar RealScalar;
00063     typedef typename _MatrixType::Index Index;
00064   };
00065 
00066   template<typename _MatrixType, int Options>
00067   struct pastix_traits< PastixLDLT<_MatrixType,Options> >
00068   {
00069     typedef _MatrixType MatrixType;
00070     typedef typename _MatrixType::Scalar Scalar;
00071     typedef typename _MatrixType::RealScalar RealScalar;
00072     typedef typename _MatrixType::Index Index;
00073   };
00074   
00075   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)
00076   {
00077     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00078     if (nbrhs == 0) x = NULL;
00079     s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
00080   }
00081   
00082   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)
00083   {
00084     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00085     if (nbrhs == 0) x = NULL;
00086     d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm); 
00087   }
00088   
00089   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)
00090   {
00091     c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm); 
00092   }
00093   
00094   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)
00095   {
00096     if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
00097     if (nbrhs == 0) x = NULL;
00098     z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); 
00099   }
00100 
00101   // Convert the matrix  to Fortran-style Numbering
00102   template <typename MatrixType>
00103   void EigenToFortranNumbering (MatrixType& mat)
00104   {
00105     if ( !(mat.outerIndexPtr()[0]) ) 
00106     { 
00107       int i;
00108       for(i = 0; i <= mat.rows(); ++i)
00109         ++mat.outerIndexPtr()[i];
00110       for(i = 0; i < mat.nonZeros(); ++i)
00111         ++mat.innerIndexPtr()[i];
00112     }
00113   }
00114   
00115   // Convert to C-style Numbering
00116   template <typename MatrixType>
00117   void EigenToCNumbering (MatrixType& mat)
00118   {
00119     // Check the Numbering
00120     if ( mat.outerIndexPtr()[0] == 1 ) 
00121     { // Convert to C-style numbering
00122       int i;
00123       for(i = 0; i <= mat.rows(); ++i)
00124         --mat.outerIndexPtr()[i];
00125       for(i = 0; i < mat.nonZeros(); ++i)
00126         --mat.innerIndexPtr()[i];
00127     }
00128   }
00129   
00130   // Symmetrize the graph of the input matrix 
00131   // In : The Input matrix to symmetrize the pattern
00132   // Out : The output matrix
00133   // StrMatTrans : The structural pattern of the transpose of In; It is 
00134   // used to optimize the future symmetrization with the same matrix pattern
00135   // WARNING It is assumed here that successive calls to this routine are done 
00136   // with matrices having the same pattern.
00137   template <typename MatrixType>
00138   void EigenSymmetrizeMatrixGraph (const MatrixType& In, MatrixType& Out, MatrixType& StrMatTrans, bool& hasTranspose)
00139   {
00140     eigen_assert(In.cols()==In.rows() && " Can only symmetrize the graph of a square matrix");
00141     if (!hasTranspose)
00142     { //First call to this routine, need to compute the structural pattern of In^T
00143       StrMatTrans = In.transpose();
00144       // Set the elements of the matrix to zero 
00145       for (int i = 0; i < StrMatTrans.rows(); i++) 
00146       {
00147         for (typename MatrixType::InnerIterator it(StrMatTrans, i); it; ++it)
00148           it.valueRef() = 0.0;
00149       }
00150       hasTranspose = true;
00151     }
00152     Out = (StrMatTrans + In).eval(); 
00153   }
00154 
00155 }
00156 
00157 // This is the base class to interface with PaStiX functions. 
00158 // Users should not used this class directly. 
00159 template <class Derived>
00160 class PastixBase
00161 {
00162   public:
00163     typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
00164     typedef _MatrixType MatrixType;
00165     typedef typename MatrixType::Scalar Scalar;
00166     typedef typename MatrixType::RealScalar RealScalar;
00167     typedef typename MatrixType::Index Index;
00168     typedef Matrix<Scalar,Dynamic,1> Vector;
00169     
00170   public:
00171     
00172     PastixBase():m_initisOk(false),m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) 
00173     {
00174       m_pastixdata = 0;
00175       m_hasTranspose = false;
00176       PastixInit();
00177     }
00178     
00179     ~PastixBase() 
00180     {
00181        PastixDestroy();
00182     }
00183     
00184     // Initialize the Pastix data structure, check the matrix
00185     void PastixInit(); 
00186     
00187     // Compute the ordering and the symbolic factorization
00188     Derived& analyzePattern (MatrixType& mat);
00189     
00190     // Compute the numerical factorization
00191     Derived& factorize (MatrixType& mat);
00192 
00197     template<typename Rhs>
00198     inline const internal::solve_retval<PastixBase, Rhs>
00199     solve(const MatrixBase<Rhs>& b) const
00200     {
00201       eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
00202       eigen_assert(rows()==b.rows()
00203                 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
00204       return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
00205     }
00206     
00207     template<typename Rhs,typename Dest>
00208     bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
00209     
00211     template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00212     void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00213     {
00214       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00215       eigen_assert(rows()==b.rows());
00216       
00217       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
00218       static const int NbColsAtOnce = 1;
00219       int rhsCols = b.cols();
00220       int size = b.rows();
00221       Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
00222       for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00223       {
00224         int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00225         tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
00226         tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
00227         dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
00228       }
00229     }
00230     
00231     Derived& derived()
00232     {
00233       return *static_cast<Derived*>(this);
00234     }
00235     const Derived& derived() const
00236     {
00237       return *static_cast<const Derived*>(this);
00238     }
00239 
00245     Array<Index,IPARM_SIZE,1>& iparm()
00246     {
00247       return m_iparm; 
00248     }
00249     
00254     int& iparm(int idxparam)
00255     {
00256       return m_iparm(idxparam);
00257     }
00258     
00263     Array<RealScalar,IPARM_SIZE,1>& dparm()
00264     {
00265       return m_dparm; 
00266     }
00267     
00268     
00273     double& dparm(int idxparam)
00274     {
00275       return m_dparm(idxparam);
00276     }
00277     
00278     inline Index cols() const { return m_size; }
00279     inline Index rows() const { return m_size; }
00280     
00289     ComputationInfo info() const
00290     {
00291       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00292       return m_info;
00293     }
00294     
00299     template<typename Rhs>
00300     inline const internal::sparse_solve_retval<PastixBase, Rhs>
00301     solve(const SparseMatrixBase<Rhs>& b) const
00302     {
00303       eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
00304       eigen_assert(rows()==b.rows()
00305                 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
00306       return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
00307     }
00308     
00309   protected:
00310     // Free all the data allocated by Pastix
00311     void PastixDestroy()
00312     {
00313       eigen_assert(m_initisOk && "The Pastix structure should be allocated first"); 
00314       m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
00315       m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
00316       internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00317                m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
00318     }
00319     
00320     Derived& compute (MatrixType& mat);
00321     
00322     int m_initisOk; 
00323     int m_analysisIsOk;
00324     int m_factorizationIsOk;
00325     bool m_isInitialized;
00326     mutable ComputationInfo m_info; 
00327     mutable pastix_data_t *m_pastixdata; // Data structure for pastix
00328     mutable SparseMatrix<Scalar, ColMajor> m_mat_null; // An input  null matrix
00329     mutable Matrix<Scalar, Dynamic,1> m_vec_null; // An input null vector
00330     mutable SparseMatrix<Scalar, ColMajor> m_StrMatTrans; // The transpose pattern of the input matrix
00331     mutable bool m_hasTranspose; // The transpose of the current matrix has already been computed
00332     mutable int m_comm; // The MPI communicator identifier
00333     mutable Matrix<Index,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
00334     mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
00335     mutable Matrix<Index,Dynamic,1> m_perm;  // Permutation vector
00336     mutable Matrix<Index,Dynamic,1> m_invp;  // Inverse permutation vector
00337     mutable int m_ordering; // ordering method to use
00338     mutable int m_amalgamation; // level of amalgamation
00339     mutable int m_size; // Size of the matrix 
00340     
00341   private:
00342     PastixBase(PastixBase& ) {}
00343     
00344 }; 
00345 
00350 template <class Derived>
00351 void PastixBase<Derived>::PastixInit()
00352 {
00353   m_size = 0; 
00354   m_iparm.resize(IPARM_SIZE);
00355   m_dparm.resize(DPARM_SIZE);
00356   
00357   m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
00358   if(m_pastixdata)
00359   { // This trick is used to reset the Pastix internal data between successive
00360     // calls with (structural) different matrices
00361     PastixDestroy();
00362     m_pastixdata = 0;
00363     m_iparm(IPARM_MODIFY_PARAMETER) = API_YES;
00364     m_hasTranspose = false; 
00365   }
00366   
00367   m_iparm(IPARM_START_TASK) = API_TASK_INIT;
00368   m_iparm(IPARM_END_TASK) = API_TASK_INIT;
00369   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00370   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00371                m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 1, m_iparm.data(), m_dparm.data());
00372   
00373   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00374   
00375   // Check the returned error
00376   if(m_iparm(IPARM_ERROR_NUMBER)) {
00377     m_info = InvalidInput;
00378     m_initisOk = false;
00379   }
00380   else { 
00381     m_info = Success;
00382     m_initisOk = true;
00383   }
00384 }
00385 
00386 template <class Derived>
00387 Derived& PastixBase<Derived>::compute(MatrixType& mat)
00388 {
00389   eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
00390   typedef typename MatrixType::Scalar Scalar;
00391   
00392   // Save the size of the current matrix 
00393   m_size = mat.rows();
00394   // Convert the matrix in fortran-style numbering 
00395   internal::EigenToFortranNumbering(mat);
00396   analyzePattern(mat);
00397   factorize(mat);
00398   m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
00399   if (m_factorizationIsOk) m_isInitialized = true;
00400   
00401   //Convert back the matrix -- Is it really necessary here 
00402   internal::EigenToCNumbering(mat);
00403   
00404   return derived();
00405 }
00406 
00407 
00408 template <class Derived>
00409 Derived& PastixBase<Derived>::analyzePattern(MatrixType& mat)
00410 {
00411   eigen_assert(m_initisOk && "PastixInit should be called first to set the default  parameters");
00412   m_size = mat.rows();
00413   m_perm.resize(m_size);
00414   m_invp.resize(m_size);
00415   
00416   // Convert the matrix in fortran-style numbering 
00417   internal::EigenToFortranNumbering(mat);
00418   
00419   m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
00420   m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
00421   
00422   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00423                mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
00424   
00425   // Check the returned error
00426   if(m_iparm(IPARM_ERROR_NUMBER)) {
00427     m_info = NumericalIssue;
00428     m_analysisIsOk = false;
00429   }
00430   else { 
00431     m_info = Success;
00432     m_analysisIsOk = true;
00433   }
00434   return derived();
00435 }
00436 
00437 template <class Derived>
00438 Derived& PastixBase<Derived>::factorize(MatrixType& mat)
00439 {
00440   eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
00441   m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
00442   m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
00443   m_size = mat.rows();
00444   
00445   // Convert the matrix in fortran-style numbering 
00446   internal::EigenToFortranNumbering(mat);
00447   
00448   internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
00449                mat.valuePtr(), m_perm.data(), m_invp.data(), m_vec_null.data(), 0, m_iparm.data(), m_dparm.data());
00450   
00451   // Check the returned error
00452   if(m_iparm(IPARM_ERROR_NUMBER)) {
00453     m_info = NumericalIssue;
00454     m_factorizationIsOk = false;
00455     m_isInitialized = false;
00456   }
00457   else { 
00458     m_info = Success;
00459     m_factorizationIsOk = true;
00460     m_isInitialized = true;
00461   }
00462   return derived();
00463 }
00464 
00465 /* Solve the system */
00466 template<typename Base>
00467 template<typename Rhs,typename Dest>
00468 bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
00469 {
00470   eigen_assert(m_isInitialized && "The matrix should be factorized first");
00471   EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
00472                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00473   int rhs = 1;
00474   
00475   x = b; /* on return, x is overwritten by the computed solution */
00476   
00477   for (int i = 0; i < b.cols(); i++){
00478     m_iparm(IPARM_START_TASK) = API_TASK_SOLVE;
00479     m_iparm(IPARM_END_TASK) = API_TASK_REFINE;
00480     m_iparm(IPARM_RHS_MAKING) = API_RHS_B;
00481     internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), m_mat_null.outerIndexPtr(), m_mat_null.innerIndexPtr(),
00482                m_mat_null.valuePtr(), m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
00483   }
00484   // Check the returned error
00485   if(m_iparm(IPARM_ERROR_NUMBER)) {
00486     m_info = NumericalIssue; 
00487     return false;
00488   }
00489   else {
00490     return true;
00491   }
00492 }
00493 
00513 template<typename _MatrixType, bool IsStrSym>
00514 class PastixLU : public PastixBase< PastixLU<_MatrixType> >
00515 {
00516   public:
00517     typedef _MatrixType MatrixType;
00518     typedef PastixBase<PastixLU<MatrixType> > Base;
00519     typedef typename MatrixType::Scalar Scalar;
00520     typedef SparseMatrix<Scalar, ColMajor> PaStiXType; 
00521     
00522   public:
00523     PastixLU():Base() {}
00524     
00525     PastixLU(const MatrixType& matrix):Base()
00526     {
00527       compute(matrix);
00528     }
00534     void compute (const MatrixType& matrix)
00535     {
00536       // Pastix supports only column-major matrices with a symmetric pattern
00537       Base::PastixInit(); 
00538       PaStiXType temp(matrix.rows(), matrix.cols());
00539       // Symmetrize the graph of the matrix
00540       if (IsStrSym)   
00541         temp = matrix;
00542       else 
00543       { 
00544         internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans, m_hasTranspose);
00545       }
00546       m_iparm[IPARM_SYM] = API_SYM_NO;
00547       m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00548       Base::compute(temp);
00549     }
00555     void analyzePattern(const MatrixType& matrix)
00556     {
00557       
00558       Base::PastixInit(); 
00559       /* Pastix supports only column-major matrices with symmetrized patterns */
00560       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00561       // Symmetrize the graph of the matrix
00562       if (IsStrSym)   
00563         temp = matrix;
00564       else 
00565       { 
00566         internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
00567       }
00568       
00569       m_iparm(IPARM_SYM) = API_SYM_NO;
00570       m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00571       Base::analyzePattern(temp);
00572     }
00573 
00579     void factorize(const MatrixType& matrix)
00580     {
00581       /* Pastix supports only column-major matrices with symmetrized patterns */
00582       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00583       // Symmetrize the graph of the matrix
00584       if (IsStrSym)   
00585         temp = matrix;
00586       else 
00587       { 
00588         internal::EigenSymmetrizeMatrixGraph<PaStiXType>(matrix, temp, m_StrMatTrans,m_hasTranspose);
00589       }
00590       m_iparm(IPARM_SYM) = API_SYM_NO;
00591       m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
00592       Base::factorize(temp);
00593     }
00594   protected:
00595     using Base::m_iparm;
00596     using Base::m_dparm;
00597     using Base::m_StrMatTrans;
00598     using Base::m_hasTranspose;
00599     
00600   private:
00601     PastixLU(PastixLU& ) {}
00602 };
00603 
00618 template<typename _MatrixType, int _UpLo>
00619 class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
00620 {
00621   public:
00622     typedef _MatrixType MatrixType;
00623     typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
00624     typedef typename MatrixType::Scalar Scalar;
00625     typedef typename MatrixType::Index Index;
00626     
00627   public:
00628     enum { UpLo = _UpLo };
00629     PastixLLT():Base() {}
00630     
00631     PastixLLT(const MatrixType& matrix):Base()
00632     {
00633       compute(matrix);
00634     }
00635 
00639     void compute (const MatrixType& matrix)
00640     {
00641       // Pastix supports only lower, column-major matrices
00642       Base::PastixInit(); // This is necessary to let PaStiX initialize its data structure between successive calls to compute
00643       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00644       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00645       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00646       m_iparm(IPARM_SYM) = API_SYM_YES;
00647       m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00648       Base::compute(temp);
00649     }
00650 
00655     void analyzePattern(const MatrixType& matrix)
00656     {
00657       Base::PastixInit(); 
00658       // Pastix supports only lower, column-major matrices
00659       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00660       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00661       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00662       m_iparm(IPARM_SYM) = API_SYM_YES;
00663       m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00664       Base::analyzePattern(temp);
00665     }
00669     void factorize(const MatrixType& matrix)
00670     {
00671       // Pastix supports only lower, column-major matrices 
00672       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00673       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00674       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00675       m_iparm(IPARM_SYM) = API_SYM_YES;
00676       m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
00677       Base::factorize(temp);
00678     }
00679   protected:
00680     using Base::m_iparm;
00681     
00682   private:
00683     PastixLLT(PastixLLT& ) {}
00684 };
00685 
00700 template<typename _MatrixType, int _UpLo>
00701 class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
00702 {
00703 public:
00704     typedef _MatrixType MatrixType;
00705     typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base; 
00706     typedef typename MatrixType::Scalar Scalar;
00707     typedef typename MatrixType::Index Index;
00708     
00709   public:
00710     enum { UpLo = _UpLo };
00711     PastixLDLT():Base() {}
00712     
00713     PastixLDLT(const MatrixType& matrix):Base()
00714     {
00715       compute(matrix);
00716     }
00717 
00721     void compute (const MatrixType& matrix)
00722     {
00723       Base::PastixInit();
00724       // Pastix supports only lower, column-major matrices
00725       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00726       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00727       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00728       m_iparm(IPARM_SYM) = API_SYM_YES;
00729       m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00730       Base::compute(temp);
00731     }
00732 
00737     void analyzePattern(const MatrixType& matrix)
00738     { 
00739       Base::PastixInit();
00740       // Pastix supports only lower, column-major matrices
00741       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00742       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00743       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00744     
00745       m_iparm(IPARM_SYM) = API_SYM_YES;
00746       m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00747       Base::analyzePattern(temp);
00748     }
00752     void factorize(const MatrixType& matrix)
00753     {
00754       // Pastix supports only lower, column-major matrices
00755       SparseMatrix<Scalar, ColMajor> temp(matrix.rows(), matrix.cols());
00756       PermutationMatrix<Dynamic,Dynamic,Index> pnull;
00757       temp.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>().twistedBy(pnull);
00758 
00759       m_iparm(IPARM_SYM) = API_SYM_YES;
00760       m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
00761       Base::factorize(temp);
00762     }
00763 
00764   protected:
00765     using Base::m_iparm;
00766     
00767   private:
00768     PastixLDLT(PastixLDLT& ) {}
00769 };
00770 
00771 namespace internal {
00772 
00773 template<typename _MatrixType, typename Rhs>
00774 struct solve_retval<PastixBase<_MatrixType>, Rhs>
00775   : solve_retval_base<PastixBase<_MatrixType>, Rhs>
00776 {
00777   typedef PastixBase<_MatrixType> Dec;
00778   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00779 
00780   template<typename Dest> void evalTo(Dest& dst) const
00781   {
00782     dec()._solve(rhs(),dst);
00783   }
00784 };
00785 
00786 template<typename _MatrixType, typename Rhs>
00787 struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
00788   : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
00789 {
00790   typedef PastixBase<_MatrixType> Dec;
00791   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00792 
00793   template<typename Dest> void evalTo(Dest& dst) const
00794   {
00795     dec()._solve_sparse(rhs(),dst);
00796   }
00797 };
00798 
00799 } // end namespace internal
00800 
00801 } // end namespace Eigen
00802 
00803 #endif