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 
00209   protected:
00210     void pardisoRelease()
00211     {
00212       if(m_initialized) // Factorization ran at least once
00213       {
00214         internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
00215                                                    m_iparm.data(), m_msglvl, 0, 0);
00216       }
00217     }
00218 
00219     void pardisoInit(int type)
00220     {
00221       m_type = type;
00222       bool symmetric = abs(m_type) < 10;
00223       m_iparm[0] = 1;   // No solver default
00224       m_iparm[1] = 3;   // use Metis for the ordering
00225       m_iparm[2] = 1;   // Numbers of processors, value of OMP_NUM_THREADS
00226       m_iparm[3] = 0;   // No iterative-direct algorithm
00227       m_iparm[4] = 0;   // No user fill-in reducing permutation
00228       m_iparm[5] = 0;   // Write solution into x
00229       m_iparm[6] = 0;   // Not in use
00230       m_iparm[7] = 2;   // Max numbers of iterative refinement steps
00231       m_iparm[8] = 0;   // Not in use
00232       m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
00233       m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
00234       m_iparm[11] = 0;  // Not in use
00235       m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
00236                                         // Try m_iparm[12] = 1 in case of inappropriate accuracy
00237       m_iparm[13] = 0;  // Output: Number of perturbed pivots
00238       m_iparm[14] = 0;  // Not in use
00239       m_iparm[15] = 0;  // Not in use
00240       m_iparm[16] = 0;  // Not in use
00241       m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
00242       m_iparm[18] = -1; // Output: Mflops for LU factorization
00243       m_iparm[19] = 0;  // Output: Numbers of CG Iterations
00244       
00245       m_iparm[20] = 0;  // 1x1 pivoting
00246       m_iparm[26] = 0;  // No matrix checker
00247       m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
00248       m_iparm[34] = 1;  // C indexing
00249       m_iparm[59] = 1;  // Automatic switch between In-Core and Out-of-Core modes
00250     }
00251 
00252   protected:
00253     // cached data to reduce reallocation, etc.
00254     
00255     void manageErrorCode(Index error)
00256     {
00257       switch(error)
00258       {
00259         case 0:
00260           m_info = Success;
00261           break;
00262         case -4:
00263         case -7:
00264           m_info = NumericalIssue;
00265           break;
00266         default:
00267           m_info = InvalidInput;
00268       }
00269     }
00270 
00271     mutable SparseMatrixType m_matrix;
00272     ComputationInfo m_info;
00273     bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
00274     Index m_type, m_msglvl;
00275     mutable void *m_pt[64];
00276     mutable ParameterType m_iparm;
00277     mutable IntColVectorType m_perm;
00278     Index m_size;
00279     
00280   private:
00281     PardisoImpl(PardisoImpl &) {}
00282 };
00283 
00284 template<class Derived>
00285 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
00286 {
00287   m_size = a.rows();
00288   eigen_assert(a.rows() == a.cols());
00289 
00290   pardisoRelease();
00291   memset(m_pt, 0, sizeof(m_pt));
00292   m_perm.setZero(m_size);
00293   derived().getMatrix(a);
00294   
00295   Index error;
00296   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
00297                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00298                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00299 
00300   manageErrorCode(error);
00301   m_analysisIsOk = true;
00302   m_factorizationIsOk = true;
00303   m_initialized = true;
00304   return derived();
00305 }
00306 
00307 template<class Derived>
00308 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
00309 {
00310   m_size = a.rows();
00311   eigen_assert(m_size == 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, 11, 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 = false;
00326   m_initialized = true;
00327   return derived();
00328 }
00329 
00330 template<class Derived>
00331 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
00332 {
00333   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00334   eigen_assert(m_size == a.rows() && m_size == a.cols());
00335   
00336   derived().getMatrix(a);
00337 
00338   Index error;  
00339   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
00340                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00341                                                      m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00342   
00343   manageErrorCode(error);
00344   m_factorizationIsOk = true;
00345   return derived();
00346 }
00347 
00348 template<class Base>
00349 template<typename BDerived,typename XDerived>
00350 bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
00351 {
00352   if(m_iparm[0] == 0) // Factorization was not computed
00353     return false;
00354 
00355   //Index n = m_matrix.rows();
00356   Index nrhs = Index(b.cols());
00357   eigen_assert(m_size==b.rows());
00358   eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
00359   eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
00360   eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
00361 
00362 
00363 //  switch (transposed) {
00364 //    case SvNoTrans    : m_iparm[11] = 0 ; break;
00365 //    case SvTranspose  : m_iparm[11] = 2 ; break;
00366 //    case SvAdjoint    : m_iparm[11] = 1 ; break;
00367 //    default:
00368 //      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
00369 //      m_iparm[11] = 0;
00370 //  }
00371 
00372   Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
00373   Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
00374   
00375   // Pardiso cannot solve in-place
00376   if(rhs_ptr == x.derived().data())
00377   {
00378     tmp = b;
00379     rhs_ptr = tmp.data();
00380   }
00381   
00382   Index error;
00383   error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
00384                                                      m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00385                                                      m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
00386                                                      rhs_ptr, x.derived().data());
00387 
00388   return error==0;
00389 }
00390 
00391 
00404 template<typename MatrixType>
00405 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
00406 {
00407   protected:
00408     typedef PardisoImpl< PardisoLU<MatrixType> > Base;
00409     typedef typename Base::Scalar Scalar;
00410     typedef typename Base::RealScalar RealScalar;
00411     using Base::pardisoInit;
00412     using Base::m_matrix;
00413     friend class PardisoImpl< PardisoLU<MatrixType> >;
00414 
00415   public:
00416 
00417     using Base::compute;
00418     using Base::solve;
00419 
00420     PardisoLU()
00421       : Base()
00422     {
00423       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00424     }
00425 
00426     PardisoLU(const MatrixType& matrix)
00427       : Base()
00428     {
00429       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00430       compute(matrix);
00431     }
00432   protected:
00433     void getMatrix(const MatrixType& matrix)
00434     {
00435       m_matrix = matrix;
00436     }
00437     
00438   private:
00439     PardisoLU(PardisoLU& ) {}
00440 };
00441 
00456 template<typename MatrixType, int _UpLo>
00457 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
00458 {
00459   protected:
00460     typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
00461     typedef typename Base::Scalar Scalar;
00462     typedef typename Base::Index Index;
00463     typedef typename Base::RealScalar RealScalar;
00464     using Base::pardisoInit;
00465     using Base::m_matrix;
00466     friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
00467 
00468   public:
00469 
00470     enum { UpLo = _UpLo };
00471     using Base::compute;
00472     using Base::solve;
00473 
00474     PardisoLLT()
00475       : Base()
00476     {
00477       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00478     }
00479 
00480     PardisoLLT(const MatrixType& matrix)
00481       : Base()
00482     {
00483       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00484       compute(matrix);
00485     }
00486     
00487   protected:
00488     
00489     void getMatrix(const MatrixType& matrix)
00490     {
00491       // PARDISO supports only upper, row-major matrices
00492       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00493       m_matrix.resize(matrix.rows(), matrix.cols());
00494       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00495     }
00496     
00497   private:
00498     PardisoLLT(PardisoLLT& ) {}
00499 };
00500 
00517 template<typename MatrixType, int Options>
00518 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
00519 {
00520   protected:
00521     typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
00522     typedef typename Base::Scalar Scalar;
00523     typedef typename Base::Index Index;
00524     typedef typename Base::RealScalar RealScalar;
00525     using Base::pardisoInit;
00526     using Base::m_matrix;
00527     friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
00528 
00529   public:
00530 
00531     using Base::compute;
00532     using Base::solve;
00533     enum { UpLo = Options&(Upper|Lower) };
00534 
00535     PardisoLDLT()
00536       : Base()
00537     {
00538       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00539     }
00540 
00541     PardisoLDLT(const MatrixType& matrix)
00542       : Base()
00543     {
00544       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00545       compute(matrix);
00546     }
00547     
00548     void getMatrix(const MatrixType& matrix)
00549     {
00550       // PARDISO supports only upper, row-major matrices
00551       PermutationMatrix<Dynamic,Dynamic,Index> p_null;
00552       m_matrix.resize(matrix.rows(), matrix.cols());
00553       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00554     }
00555     
00556   private:
00557     PardisoLDLT(PardisoLDLT& ) {}
00558 };
00559 
00560 namespace internal {
00561   
00562 template<typename _Derived, typename Rhs>
00563 struct solve_retval<PardisoImpl<_Derived>, Rhs>
00564   : solve_retval_base<PardisoImpl<_Derived>, Rhs>
00565 {
00566   typedef PardisoImpl<_Derived> Dec;
00567   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00568 
00569   template<typename Dest> void evalTo(Dest& dst) const
00570   {
00571     dec()._solve(rhs(),dst);
00572   }
00573 };
00574 
00575 template<typename Derived, typename Rhs>
00576 struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
00577   : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
00578 {
00579   typedef PardisoImpl<Derived> Dec;
00580   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00581 
00582   template<typename Dest> void evalTo(Dest& dst) const
00583   {
00584     this->defaultEvalTo(dst);
00585   }
00586 };
00587 
00588 } // end namespace internal
00589 
00590 } // end namespace Eigen
00591 
00592 #endif // EIGEN_PARDISOSUPPORT_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:59:39