PardisoSupport.h
Go to the documentation of this file.
00001 /*
00002  Copyright (c) 2011, Intel Corporation. All rights reserved.
00003 
00004  Redistribution and use in source and binary forms, with or without modification,
00005  are permitted provided that the following conditions are met:
00006 
00007  * Redistributions of source code must retain the above copyright notice, this
00008    list of conditions and the following disclaimer.
00009  * Redistributions in binary form must reproduce the above copyright notice,
00010    this list of conditions and the following disclaimer in the documentation
00011    and/or other materials provided with the distribution.
00012  * Neither the name of Intel Corporation nor the names of its contributors may
00013    be used to endorse or promote products derived from this software without
00014    specific prior written permission.
00015 
00016  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
00017  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
00018  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
00019  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
00020  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00021  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00022  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
00023  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
00024  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00025  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00026 
00027  ********************************************************************************
00028  *   Content : Eigen bindings to Intel(R) MKL PARDISO
00029  ********************************************************************************
00030 */
00031 
00032 #ifndef EIGEN_PARDISOSUPPORT_H
00033 #define EIGEN_PARDISOSUPPORT_H
00034 
00035 namespace Eigen { 
00036 
00037 template<typename _MatrixType> class PardisoLU;
00038 template<typename _MatrixType, int Options=Upper> class PardisoLLT;
00039 template<typename _MatrixType, int Options=Upper> class PardisoLDLT;
00040 
00041 namespace internal
00042 {
00043   template<typename Index>
00044   struct pardiso_run_selector
00045   {
00046     static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
00047                       Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
00048     {
00049       Index error = 0;
00050       ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
00051       return error;
00052     }
00053   };
00054   template<>
00055   struct pardiso_run_selector<long long int>
00056   {
00057     typedef long long int Index;
00058     static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
00059                       Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
00060     {
00061       Index error = 0;
00062       ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
00063       return error;
00064     }
00065   };
00066 
00067   template<class Pardiso> struct pardiso_traits;
00068 
00069   template<typename _MatrixType>
00070   struct pardiso_traits< PardisoLU<_MatrixType> >
00071   {
00072     typedef _MatrixType MatrixType;
00073     typedef typename _MatrixType::Scalar Scalar;
00074     typedef typename _MatrixType::RealScalar RealScalar;
00075     typedef typename _MatrixType::Index Index;
00076   };
00077 
00078   template<typename _MatrixType, int Options>
00079   struct pardiso_traits< PardisoLLT<_MatrixType, Options> >
00080   {
00081     typedef _MatrixType MatrixType;
00082     typedef typename _MatrixType::Scalar Scalar;
00083     typedef typename _MatrixType::RealScalar RealScalar;
00084     typedef typename _MatrixType::Index Index;
00085   };
00086 
00087   template<typename _MatrixType, int Options>
00088   struct pardiso_traits< PardisoLDLT<_MatrixType, Options> >
00089   {
00090     typedef _MatrixType MatrixType;
00091     typedef typename _MatrixType::Scalar Scalar;
00092     typedef typename _MatrixType::RealScalar RealScalar;
00093     typedef typename _MatrixType::Index Index;    
00094   };
00095 
00096 }
00097 
00098 template<class Derived>
00099 class PardisoImpl
00100 {
00101     typedef internal::pardiso_traits<Derived> Traits;
00102   public:
00103     typedef typename Traits::MatrixType MatrixType;
00104     typedef typename Traits::Scalar Scalar;
00105     typedef typename Traits::RealScalar RealScalar;
00106     typedef typename Traits::Index Index;
00107     typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType;
00108     typedef Matrix<Scalar,Dynamic,1> VectorType;
00109     typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
00110     typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
00111     enum {
00112       ScalarIsComplex = NumTraits<Scalar>::IsComplex
00113     };
00114 
00115     PardisoImpl()
00116     {
00117       eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
00118       m_iparm.setZero();
00119       m_msglvl = 0; // No output
00120       m_initialized = false;
00121     }
00122 
00123     ~PardisoImpl()
00124     {
00125       pardisoRelease();
00126     }
00127 
00128     inline Index cols() const { return m_size; }
00129     inline Index rows() const { return m_size; }
00130   
00136     ComputationInfo info() const
00137     {
00138       eigen_assert(m_initialized && "Decomposition is not initialized.");
00139       return m_info;
00140     }
00141 
00145     Array<Index,64,1>& pardisoParameterArray()
00146     {
00147       return m_iparm;
00148     }
00149     
00156     Derived& analyzePattern(const MatrixType& matrix);
00157     
00164     Derived& factorize(const MatrixType& matrix);
00165 
00166     Derived& compute(const MatrixType& matrix);
00167     
00172     template<typename Rhs>
00173     inline const internal::solve_retval<PardisoImpl, Rhs>
00174     solve(const MatrixBase<Rhs>& b) const
00175     {
00176       eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00177       eigen_assert(rows()==b.rows()
00178                 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00179       return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00180     }
00181 
00186     template<typename Rhs>
00187     inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
00188     solve(const SparseMatrixBase<Rhs>& b) const
00189     {
00190       eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00191       eigen_assert(rows()==b.rows()
00192                 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00193       return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00194     }
00195 
00196     Derived& derived()
00197     {
00198       return *static_cast<Derived*>(this);
00199     }
00200     const Derived& derived() const
00201     {
00202       return *static_cast<const Derived*>(this);
00203     }
00204 
00205     template<typename BDerived, typename XDerived>
00206     bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
00207 
00209     template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00210     void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00211     {
00212       eigen_assert(m_size==b.rows());
00213 
00214       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
00215       static const int NbColsAtOnce = 4;
00216       int rhsCols = b.cols();
00217       int size = b.rows();
00218       // Pardiso cannot solve in-place,
00219       // so we need two temporaries
00220       Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols);
00221       Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols);
00222       for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00223       {
00224         int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00225         tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
00226         tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
00227         dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
00228       }
00229     }
00230 
00231   protected:
00232     void pardisoRelease()
00233     {
00234       if(m_initialized) // Factorization ran at least once
00235       {
00236         internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
00237                                                    m_iparm.data(), m_msglvl, 0, 0);
00238       }
00239     }
00240 
00241     void pardisoInit(int type)
00242     {
00243       m_type = type;
00244       bool symmetric = abs(m_type) < 10;
00245       m_iparm[0] = 1;   // No solver default
00246       m_iparm[1] = 3;   // use Metis for the ordering
00247       m_iparm[2] = 1;   // Numbers of processors, value of OMP_NUM_THREADS
00248       m_iparm[3] = 0;   // No iterative-direct algorithm
00249       m_iparm[4] = 0;   // No user fill-in reducing permutation
00250       m_iparm[5] = 0;   // Write solution into x
00251       m_iparm[6] = 0;   // Not in use
00252       m_iparm[7] = 2;   // Max numbers of iterative refinement steps
00253       m_iparm[8] = 0;   // Not in use
00254       m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
00255       m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
00256       m_iparm[11] = 0;  // Not in use
00257       m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
00258                                         // Try m_iparm[12] = 1 in case of inappropriate accuracy
00259       m_iparm[13] = 0;  // Output: Number of perturbed pivots
00260       m_iparm[14] = 0;  // Not in use
00261       m_iparm[15] = 0;  // Not in use
00262       m_iparm[16] = 0;  // Not in use
00263       m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
00264       m_iparm[18] = -1; // Output: Mflops for LU factorization
00265       m_iparm[19] = 0;  // Output: Numbers of CG Iterations
00266       
00267       m_iparm[20] = 0;  // 1x1 pivoting
00268       m_iparm[26] = 0;  // No matrix checker
00269       m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
00270       m_iparm[34] = 1;  // C indexing
00271       m_iparm[59] = 1;  // Automatic switch between In-Core and Out-of-Core modes
00272     }
00273 
00274   protected:
00275     // cached data to reduce reallocation, etc.
00276     
00277     void manageErrorCode(Index error)
00278     {
00279       switch(error)
00280       {
00281         case 0:
00282           m_info = Success;
00283           break;
00284         case -4:
00285         case -7:
00286           m_info = NumericalIssue;
00287           break;
00288         default:
00289           m_info = InvalidInput;
00290       }
00291     }
00292 
00293     mutable SparseMatrixType m_matrix;
00294     ComputationInfo m_info;
00295     bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
00296     Index m_type, m_msglvl;
00297     mutable void *m_pt[64];
00298     mutable Array<Index,64,1> m_iparm;
00299     mutable IntColVectorType m_perm;
00300     Index m_size;
00301     
00302   private:
00303     PardisoImpl(PardisoImpl &) {}
00304 };
00305 
00306 template<class Derived>
00307 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
00308 {
00309   m_size = a.rows();
00310   eigen_assert(a.rows() == a.cols());
00311 
00312   pardisoRelease();
00313   memset(m_pt, 0, sizeof(m_pt));
00314   m_perm.setZero(m_size);
00315   derived().getMatrix(a);
00316   
00317   Index error;
00318   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
00319                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00320                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00321 
00322   manageErrorCode(error);
00323   m_analysisIsOk = true;
00324   m_factorizationIsOk = true;
00325   m_initialized = true;
00326   return derived();
00327 }
00328 
00329 template<class Derived>
00330 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
00331 {
00332   m_size = a.rows();
00333   eigen_assert(m_size == a.cols());
00334 
00335   pardisoRelease();
00336   memset(m_pt, 0, sizeof(m_pt));
00337   m_perm.setZero(m_size);
00338   derived().getMatrix(a);
00339   
00340   Index error;
00341   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
00342                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00343                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00344   
00345   manageErrorCode(error);
00346   m_analysisIsOk = true;
00347   m_factorizationIsOk = false;
00348   m_initialized = true;
00349   return derived();
00350 }
00351 
00352 template<class Derived>
00353 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
00354 {
00355   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00356   eigen_assert(m_size == a.rows() && m_size == a.cols());
00357   
00358   derived().getMatrix(a);
00359 
00360   Index error;  
00361   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
00362                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00363                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00364   
00365   manageErrorCode(error);
00366   m_factorizationIsOk = true;
00367   return derived();
00368 }
00369 
00370 template<class Base>
00371 template<typename BDerived,typename XDerived>
00372 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
00373 {
00374   if(m_iparm[0] == 0) // Factorization was not computed
00375     return false;
00376 
00377   //Index n = m_matrix.rows();
00378   Index nrhs = Index(b.cols());
00379   eigen_assert(m_size==b.rows());
00380   eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
00381   eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
00382   eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
00383 
00384 
00385 //  switch (transposed) {
00386 //    case SvNoTrans    : m_iparm[11] = 0 ; break;
00387 //    case SvTranspose  : m_iparm[11] = 2 ; break;
00388 //    case SvAdjoint    : m_iparm[11] = 1 ; break;
00389 //    default:
00390 //      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
00391 //      m_iparm[11] = 0;
00392 //  }
00393 
00394   Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
00395   Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
00396   
00397   // Pardiso cannot solve in-place
00398   if(rhs_ptr == x.derived().data())
00399   {
00400     tmp = b;
00401     rhs_ptr = tmp.data();
00402   }
00403   
00404   Index error;
00405   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
00406                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00407                                                      m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
00408                                                      rhs_ptr, x.derived().data());
00409 
00410   return error==0;
00411 }
00412 
00413 
00426 template<typename MatrixType>
00427 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
00428 {
00429   protected:
00430     typedef PardisoImpl< PardisoLU<MatrixType> > Base;
00431     typedef typename Base::Scalar Scalar;
00432     typedef typename Base::RealScalar RealScalar;
00433     using Base::pardisoInit;
00434     using Base::m_matrix;
00435     friend class PardisoImpl< PardisoLU<MatrixType> >;
00436 
00437   public:
00438 
00439     using Base::compute;
00440     using Base::solve;
00441 
00442     PardisoLU()
00443       : Base()
00444     {
00445       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00446     }
00447 
00448     PardisoLU(const MatrixType& matrix)
00449       : Base()
00450     {
00451       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00452       compute(matrix);
00453     }
00454   protected:
00455     void getMatrix(const MatrixType& matrix)
00456     {
00457       m_matrix = matrix;
00458     }
00459     
00460   private:
00461     PardisoLU(PardisoLU& ) {}
00462 };
00463 
00478 template<typename MatrixType, int _UpLo>
00479 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
00480 {
00481   protected:
00482     typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
00483     typedef typename Base::Scalar Scalar;
00484     typedef typename Base::Index Index;
00485     typedef typename Base::RealScalar RealScalar;
00486     using Base::pardisoInit;
00487     using Base::m_matrix;
00488     friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
00489 
00490   public:
00491 
00492     enum { UpLo = _UpLo };
00493     using Base::compute;
00494     using Base::solve;
00495 
00496     PardisoLLT()
00497       : Base()
00498     {
00499       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00500     }
00501 
00502     PardisoLLT(const MatrixType& matrix)
00503       : Base()
00504     {
00505       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00506       compute(matrix);
00507     }
00508     
00509   protected:
00510     
00511     void getMatrix(const MatrixType& matrix)
00512     {
00513       // PARDISO supports only upper, row-major matrices
00514       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00515       m_matrix.resize(matrix.rows(), matrix.cols());
00516       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00517     }
00518     
00519   private:
00520     PardisoLLT(PardisoLLT& ) {}
00521 };
00522 
00539 template<typename MatrixType, int Options>
00540 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
00541 {
00542   protected:
00543     typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
00544     typedef typename Base::Scalar Scalar;
00545     typedef typename Base::Index Index;
00546     typedef typename Base::RealScalar RealScalar;
00547     using Base::pardisoInit;
00548     using Base::m_matrix;
00549     friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
00550 
00551   public:
00552 
00553     using Base::compute;
00554     using Base::solve;
00555     enum { UpLo = Options&(Upper|Lower) };
00556 
00557     PardisoLDLT()
00558       : Base()
00559     {
00560       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00561     }
00562 
00563     PardisoLDLT(const MatrixType& matrix)
00564       : Base()
00565     {
00566       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00567       compute(matrix);
00568     }
00569     
00570     void getMatrix(const MatrixType& matrix)
00571     {
00572       // PARDISO supports only upper, row-major matrices
00573       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00574       m_matrix.resize(matrix.rows(), matrix.cols());
00575       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00576     }
00577     
00578   private:
00579     PardisoLDLT(PardisoLDLT& ) {}
00580 };
00581 
00582 namespace internal {
00583   
00584 template<typename _Derived, typename Rhs>
00585 struct solve_retval<PardisoImpl<_Derived>, Rhs>
00586   : solve_retval_base<PardisoImpl<_Derived>, Rhs>
00587 {
00588   typedef PardisoImpl<_Derived> Dec;
00589   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00590 
00591   template<typename Dest> void evalTo(Dest& dst) const
00592   {
00593     dec()._solve(rhs(),dst);
00594   }
00595 };
00596 
00597 template<typename Derived, typename Rhs>
00598 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
00599   : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
00600 {
00601   typedef PardisoImpl<Derived> Dec;
00602   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00603 
00604   template<typename Dest> void evalTo(Dest& dst) const
00605   {
00606     dec().derived()._solve_sparse(rhs(),dst);
00607   }
00608 };
00609 
00610 } // end namespace internal
00611 
00612 } // end namespace Eigen
00613 
00614 #endif // EIGEN_PARDISOSUPPORT_H