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     typedef Array<Index,64,1,DontAlign> ParameterType;
00112     enum {
00113       ScalarIsComplex = NumTraits<Scalar>::IsComplex
00114     };
00115 
00116     PardisoImpl()
00117     {
00118       eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
00119       m_iparm.setZero();
00120       m_msglvl = 0; // No output
00121       m_initialized = false;
00122     }
00123 
00124     ~PardisoImpl()
00125     {
00126       pardisoRelease();
00127     }
00128 
00129     inline Index cols() const { return m_size; }
00130     inline Index rows() const { return m_size; }
00131   
00137     ComputationInfo info() const
00138     {
00139       eigen_assert(m_initialized && "Decomposition is not initialized.");
00140       return m_info;
00141     }
00142 
00146     ParameterType& pardisoParameterArray()
00147     {
00148       return m_iparm;
00149     }
00150     
00157     Derived& analyzePattern(const MatrixType& matrix);
00158     
00165     Derived& factorize(const MatrixType& matrix);
00166 
00167     Derived& compute(const MatrixType& matrix);
00168     
00173     template<typename Rhs>
00174     inline const internal::solve_retval<PardisoImpl, Rhs>
00175     solve(const MatrixBase<Rhs>& b) const
00176     {
00177       eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00178       eigen_assert(rows()==b.rows()
00179                 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00180       return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00181     }
00182 
00187     template<typename Rhs>
00188     inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
00189     solve(const SparseMatrixBase<Rhs>& b) const
00190     {
00191       eigen_assert(m_initialized && "Pardiso solver is not initialized.");
00192       eigen_assert(rows()==b.rows()
00193                 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
00194       return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
00195     }
00196 
00197     Derived& derived()
00198     {
00199       return *static_cast<Derived*>(this);
00200     }
00201     const Derived& derived() const
00202     {
00203       return *static_cast<const Derived*>(this);
00204     }
00205 
00206     template<typename BDerived, typename XDerived>
00207     bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
00208 
00210     template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00211     void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00212     {
00213       eigen_assert(m_size==b.rows());
00214 
00215       // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
00216       static const int NbColsAtOnce = 4;
00217       int rhsCols = b.cols();
00218       int size = b.rows();
00219       // Pardiso cannot solve in-place,
00220       // so we need two temporaries
00221       Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_rhs(size,rhsCols);
00222       Eigen::Matrix<DestScalar,Dynamic,Dynamic,ColMajor> tmp_res(size,rhsCols);
00223       for(int k=0; k<rhsCols; k+=NbColsAtOnce)
00224       {
00225         int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
00226         tmp_rhs.leftCols(actualCols) = b.middleCols(k,actualCols);
00227         tmp_res.leftCols(actualCols) = derived().solve(tmp_rhs.leftCols(actualCols));
00228         dest.middleCols(k,actualCols) = tmp_res.leftCols(actualCols).sparseView();
00229       }
00230     }
00231 
00232   protected:
00233     void pardisoRelease()
00234     {
00235       if(m_initialized) // Factorization ran at least once
00236       {
00237         internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
00238                                                    m_iparm.data(), m_msglvl, 0, 0);
00239       }
00240     }
00241 
00242     void pardisoInit(int type)
00243     {
00244       m_type = type;
00245       bool symmetric = abs(m_type) < 10;
00246       m_iparm[0] = 1;   // No solver default
00247       m_iparm[1] = 3;   // use Metis for the ordering
00248       m_iparm[2] = 1;   // Numbers of processors, value of OMP_NUM_THREADS
00249       m_iparm[3] = 0;   // No iterative-direct algorithm
00250       m_iparm[4] = 0;   // No user fill-in reducing permutation
00251       m_iparm[5] = 0;   // Write solution into x
00252       m_iparm[6] = 0;   // Not in use
00253       m_iparm[7] = 2;   // Max numbers of iterative refinement steps
00254       m_iparm[8] = 0;   // Not in use
00255       m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
00256       m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
00257       m_iparm[11] = 0;  // Not in use
00258       m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
00259                                         // Try m_iparm[12] = 1 in case of inappropriate accuracy
00260       m_iparm[13] = 0;  // Output: Number of perturbed pivots
00261       m_iparm[14] = 0;  // Not in use
00262       m_iparm[15] = 0;  // Not in use
00263       m_iparm[16] = 0;  // Not in use
00264       m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
00265       m_iparm[18] = -1; // Output: Mflops for LU factorization
00266       m_iparm[19] = 0;  // Output: Numbers of CG Iterations
00267       
00268       m_iparm[20] = 0;  // 1x1 pivoting
00269       m_iparm[26] = 0;  // No matrix checker
00270       m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
00271       m_iparm[34] = 1;  // C indexing
00272       m_iparm[59] = 1;  // Automatic switch between In-Core and Out-of-Core modes
00273     }
00274 
00275   protected:
00276     // cached data to reduce reallocation, etc.
00277     
00278     void manageErrorCode(Index error)
00279     {
00280       switch(error)
00281       {
00282         case 0:
00283           m_info = Success;
00284           break;
00285         case -4:
00286         case -7:
00287           m_info = NumericalIssue;
00288           break;
00289         default:
00290           m_info = InvalidInput;
00291       }
00292     }
00293 
00294     mutable SparseMatrixType m_matrix;
00295     ComputationInfo m_info;
00296     bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
00297     Index m_type, m_msglvl;
00298     mutable void *m_pt[64];
00299     mutable ParameterType m_iparm;
00300     mutable IntColVectorType m_perm;
00301     Index m_size;
00302     
00303   private:
00304     PardisoImpl(PardisoImpl &) {}
00305 };
00306 
00307 template<class Derived>
00308 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
00309 {
00310   m_size = a.rows();
00311   eigen_assert(a.rows() == a.cols());
00312 
00313   pardisoRelease();
00314   memset(m_pt, 0, sizeof(m_pt));
00315   m_perm.setZero(m_size);
00316   derived().getMatrix(a);
00317   
00318   Index error;
00319   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
00320                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00321                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00322 
00323   manageErrorCode(error);
00324   m_analysisIsOk = true;
00325   m_factorizationIsOk = true;
00326   m_initialized = true;
00327   return derived();
00328 }
00329 
00330 template<class Derived>
00331 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
00332 {
00333   m_size = a.rows();
00334   eigen_assert(m_size == a.cols());
00335 
00336   pardisoRelease();
00337   memset(m_pt, 0, sizeof(m_pt));
00338   m_perm.setZero(m_size);
00339   derived().getMatrix(a);
00340   
00341   Index error;
00342   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
00343                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00344                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00345   
00346   manageErrorCode(error);
00347   m_analysisIsOk = true;
00348   m_factorizationIsOk = false;
00349   m_initialized = true;
00350   return derived();
00351 }
00352 
00353 template<class Derived>
00354 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
00355 {
00356   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00357   eigen_assert(m_size == a.rows() && m_size == a.cols());
00358   
00359   derived().getMatrix(a);
00360 
00361   Index error;  
00362   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
00363                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00364                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00365   
00366   manageErrorCode(error);
00367   m_factorizationIsOk = true;
00368   return derived();
00369 }
00370 
00371 template<class Base>
00372 template<typename BDerived,typename XDerived>
00373 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
00374 {
00375   if(m_iparm[0] == 0) // Factorization was not computed
00376     return false;
00377 
00378   //Index n = m_matrix.rows();
00379   Index nrhs = Index(b.cols());
00380   eigen_assert(m_size==b.rows());
00381   eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
00382   eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
00383   eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
00384 
00385 
00386 //  switch (transposed) {
00387 //    case SvNoTrans    : m_iparm[11] = 0 ; break;
00388 //    case SvTranspose  : m_iparm[11] = 2 ; break;
00389 //    case SvAdjoint    : m_iparm[11] = 1 ; break;
00390 //    default:
00391 //      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
00392 //      m_iparm[11] = 0;
00393 //  }
00394 
00395   Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
00396   Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
00397   
00398   // Pardiso cannot solve in-place
00399   if(rhs_ptr == x.derived().data())
00400   {
00401     tmp = b;
00402     rhs_ptr = tmp.data();
00403   }
00404   
00405   Index error;
00406   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
00407                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00408                                                      m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
00409                                                      rhs_ptr, x.derived().data());
00410 
00411   return error==0;
00412 }
00413 
00414 
00427 template<typename MatrixType>
00428 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
00429 {
00430   protected:
00431     typedef PardisoImpl< PardisoLU<MatrixType> > Base;
00432     typedef typename Base::Scalar Scalar;
00433     typedef typename Base::RealScalar RealScalar;
00434     using Base::pardisoInit;
00435     using Base::m_matrix;
00436     friend class PardisoImpl< PardisoLU<MatrixType> >;
00437 
00438   public:
00439 
00440     using Base::compute;
00441     using Base::solve;
00442 
00443     PardisoLU()
00444       : Base()
00445     {
00446       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00447     }
00448 
00449     PardisoLU(const MatrixType& matrix)
00450       : Base()
00451     {
00452       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00453       compute(matrix);
00454     }
00455   protected:
00456     void getMatrix(const MatrixType& matrix)
00457     {
00458       m_matrix = matrix;
00459     }
00460     
00461   private:
00462     PardisoLU(PardisoLU& ) {}
00463 };
00464 
00479 template<typename MatrixType, int _UpLo>
00480 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
00481 {
00482   protected:
00483     typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
00484     typedef typename Base::Scalar Scalar;
00485     typedef typename Base::Index Index;
00486     typedef typename Base::RealScalar RealScalar;
00487     using Base::pardisoInit;
00488     using Base::m_matrix;
00489     friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
00490 
00491   public:
00492 
00493     enum { UpLo = _UpLo };
00494     using Base::compute;
00495     using Base::solve;
00496 
00497     PardisoLLT()
00498       : Base()
00499     {
00500       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00501     }
00502 
00503     PardisoLLT(const MatrixType& matrix)
00504       : Base()
00505     {
00506       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00507       compute(matrix);
00508     }
00509     
00510   protected:
00511     
00512     void getMatrix(const MatrixType& matrix)
00513     {
00514       // PARDISO supports only upper, row-major matrices
00515       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00516       m_matrix.resize(matrix.rows(), matrix.cols());
00517       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00518     }
00519     
00520   private:
00521     PardisoLLT(PardisoLLT& ) {}
00522 };
00523 
00540 template<typename MatrixType, int Options>
00541 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
00542 {
00543   protected:
00544     typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
00545     typedef typename Base::Scalar Scalar;
00546     typedef typename Base::Index Index;
00547     typedef typename Base::RealScalar RealScalar;
00548     using Base::pardisoInit;
00549     using Base::m_matrix;
00550     friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
00551 
00552   public:
00553 
00554     using Base::compute;
00555     using Base::solve;
00556     enum { UpLo = Options&(Upper|Lower) };
00557 
00558     PardisoLDLT()
00559       : Base()
00560     {
00561       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00562     }
00563 
00564     PardisoLDLT(const MatrixType& matrix)
00565       : Base()
00566     {
00567       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00568       compute(matrix);
00569     }
00570     
00571     void getMatrix(const MatrixType& matrix)
00572     {
00573       // PARDISO supports only upper, row-major matrices
00574       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00575       m_matrix.resize(matrix.rows(), matrix.cols());
00576       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00577     }
00578     
00579   private:
00580     PardisoLDLT(PardisoLDLT& ) {}
00581 };
00582 
00583 namespace internal {
00584   
00585 template<typename _Derived, typename Rhs>
00586 struct solve_retval<PardisoImpl<_Derived>, Rhs>
00587   : solve_retval_base<PardisoImpl<_Derived>, Rhs>
00588 {
00589   typedef PardisoImpl<_Derived> Dec;
00590   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00591 
00592   template<typename Dest> void evalTo(Dest& dst) const
00593   {
00594     dec()._solve(rhs(),dst);
00595   }
00596 };
00597 
00598 template<typename Derived, typename Rhs>
00599 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
00600   : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
00601 {
00602   typedef PardisoImpl<Derived> Dec;
00603   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00604 
00605   template<typename Dest> void evalTo(Dest& dst) const
00606   {
00607     dec().derived()._solve_sparse(rhs(),dst);
00608   }
00609 };
00610 
00611 } // end namespace internal
00612 
00613 } // end namespace Eigen
00614 
00615 #endif // EIGEN_PARDISOSUPPORT_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:37