Eigen  3.3.3
PardisoSupport.h
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 IndexType>
00044   struct pardiso_run_selector
00045   {
00046     static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
00047                       IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
00048     {
00049       IndexType 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 IndexType;
00058     static IndexType run( _MKL_DSS_HANDLE_t pt, IndexType maxfct, IndexType mnum, IndexType type, IndexType phase, IndexType n, void *a,
00059                       IndexType *ia, IndexType *ja, IndexType *perm, IndexType nrhs, IndexType *iparm, IndexType msglvl, void *b, void *x)
00060     {
00061       IndexType 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::StorageIndex StorageIndex;
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::StorageIndex StorageIndex;
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::StorageIndex StorageIndex;    
00094   };
00095 
00096 } // end namespace internal
00097 
00098 template<class Derived>
00099 class PardisoImpl : public SparseSolverBase<Derived>
00100 {
00101   protected:
00102     typedef SparseSolverBase<Derived> Base;
00103     using Base::derived;
00104     using Base::m_isInitialized;
00105     
00106     typedef internal::pardiso_traits<Derived> Traits;
00107   public:
00108     using Base::_solve_impl;
00109     
00110     typedef typename Traits::MatrixType MatrixType;
00111     typedef typename Traits::Scalar Scalar;
00112     typedef typename Traits::RealScalar RealScalar;
00113     typedef typename Traits::StorageIndex StorageIndex;
00114     typedef SparseMatrix<Scalar,RowMajor,StorageIndex> SparseMatrixType;
00115     typedef Matrix<Scalar,Dynamic,1> VectorType;
00116     typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
00117     typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
00118     typedef Array<StorageIndex,64,1,DontAlign> ParameterType;
00119     enum {
00120       ScalarIsComplex = NumTraits<Scalar>::IsComplex,
00121       ColsAtCompileTime = Dynamic,
00122       MaxColsAtCompileTime = Dynamic
00123     };
00124 
00125     PardisoImpl()
00126     {
00127       eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type");
00128       m_iparm.setZero();
00129       m_msglvl = 0; // No output
00130       m_isInitialized = false;
00131     }
00132 
00133     ~PardisoImpl()
00134     {
00135       pardisoRelease();
00136     }
00137 
00138     inline Index cols() const { return m_size; }
00139     inline Index rows() const { return m_size; }
00140   
00146     ComputationInfo info() const
00147     {
00148       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00149       return m_info;
00150     }
00151 
00155     ParameterType& pardisoParameterArray()
00156     {
00157       return m_iparm;
00158     }
00159     
00166     Derived& analyzePattern(const MatrixType& matrix);
00167     
00174     Derived& factorize(const MatrixType& matrix);
00175 
00176     Derived& compute(const MatrixType& matrix);
00177 
00178     template<typename Rhs,typename Dest>
00179     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
00180 
00181   protected:
00182     void pardisoRelease()
00183     {
00184       if(m_isInitialized) // Factorization ran at least once
00185       {
00186         internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, internal::convert_index<StorageIndex>(m_size),0, 0, 0, m_perm.data(), 0,
00187                                                           m_iparm.data(), m_msglvl, NULL, NULL);
00188         m_isInitialized = false;
00189       }
00190     }
00191 
00192     void pardisoInit(int type)
00193     {
00194       m_type = type;
00195       bool symmetric = std::abs(m_type) < 10;
00196       m_iparm[0] = 1;   // No solver default
00197       m_iparm[1] = 2;   // use Metis for the ordering
00198       m_iparm[2] = 0;   // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
00199       m_iparm[3] = 0;   // No iterative-direct algorithm
00200       m_iparm[4] = 0;   // No user fill-in reducing permutation
00201       m_iparm[5] = 0;   // Write solution into x, b is left unchanged
00202       m_iparm[6] = 0;   // Not in use
00203       m_iparm[7] = 2;   // Max numbers of iterative refinement steps
00204       m_iparm[8] = 0;   // Not in use
00205       m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
00206       m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
00207       m_iparm[11] = 0;  // Not in use
00208       m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
00209                                         // Try m_iparm[12] = 1 in case of inappropriate accuracy
00210       m_iparm[13] = 0;  // Output: Number of perturbed pivots
00211       m_iparm[14] = 0;  // Not in use
00212       m_iparm[15] = 0;  // Not in use
00213       m_iparm[16] = 0;  // Not in use
00214       m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
00215       m_iparm[18] = -1; // Output: Mflops for LU factorization
00216       m_iparm[19] = 0;  // Output: Numbers of CG Iterations
00217       
00218       m_iparm[20] = 0;  // 1x1 pivoting
00219       m_iparm[26] = 0;  // No matrix checker
00220       m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
00221       m_iparm[34] = 1;  // C indexing
00222       m_iparm[36] = 0;  // CSR
00223       m_iparm[59] = 0;  // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
00224       
00225       memset(m_pt, 0, sizeof(m_pt));
00226     }
00227 
00228   protected:
00229     // cached data to reduce reallocation, etc.
00230     
00231     void manageErrorCode(Index error) const
00232     {
00233       switch(error)
00234       {
00235         case 0:
00236           m_info = Success;
00237           break;
00238         case -4:
00239         case -7:
00240           m_info = NumericalIssue;
00241           break;
00242         default:
00243           m_info = InvalidInput;
00244       }
00245     }
00246 
00247     mutable SparseMatrixType m_matrix;
00248     mutable ComputationInfo m_info;
00249     bool m_analysisIsOk, m_factorizationIsOk;
00250     StorageIndex m_type, m_msglvl;
00251     mutable void *m_pt[64];
00252     mutable ParameterType m_iparm;
00253     mutable IntColVectorType m_perm;
00254     Index m_size;
00255     
00256 };
00257 
00258 template<class Derived>
00259 Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
00260 {
00261   m_size = a.rows();
00262   eigen_assert(a.rows() == a.cols());
00263 
00264   pardisoRelease();
00265   m_perm.setZero(m_size);
00266   derived().getMatrix(a);
00267   
00268   Index error;
00269   error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
00270                                                             m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00271                                                             m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00272   manageErrorCode(error);
00273   m_analysisIsOk = true;
00274   m_factorizationIsOk = true;
00275   m_isInitialized = true;
00276   return derived();
00277 }
00278 
00279 template<class Derived>
00280 Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
00281 {
00282   m_size = a.rows();
00283   eigen_assert(m_size == a.cols());
00284 
00285   pardisoRelease();
00286   m_perm.setZero(m_size);
00287   derived().getMatrix(a);
00288   
00289   Index error;
00290   error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size),
00291                                                             m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00292                                                             m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00293   
00294   manageErrorCode(error);
00295   m_analysisIsOk = true;
00296   m_factorizationIsOk = false;
00297   m_isInitialized = true;
00298   return derived();
00299 }
00300 
00301 template<class Derived>
00302 Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
00303 {
00304   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00305   eigen_assert(m_size == a.rows() && m_size == a.cols());
00306   
00307   derived().getMatrix(a);
00308 
00309   Index error;
00310   error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size),
00311                                                             m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00312                                                             m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
00313   
00314   manageErrorCode(error);
00315   m_factorizationIsOk = true;
00316   return derived();
00317 }
00318 
00319 template<class Derived>
00320 template<typename BDerived,typename XDerived>
00321 void PardisoImpl<Derived>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
00322 {
00323   if(m_iparm[0] == 0) // Factorization was not computed
00324   {
00325     m_info = InvalidInput;
00326     return;
00327   }
00328 
00329   //Index n = m_matrix.rows();
00330   Index nrhs = Index(b.cols());
00331   eigen_assert(m_size==b.rows());
00332   eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
00333   eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
00334   eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
00335 
00336 
00337 //  switch (transposed) {
00338 //    case SvNoTrans    : m_iparm[11] = 0 ; break;
00339 //    case SvTranspose  : m_iparm[11] = 2 ; break;
00340 //    case SvAdjoint    : m_iparm[11] = 1 ; break;
00341 //    default:
00342 //      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
00343 //      m_iparm[11] = 0;
00344 //  }
00345 
00346   Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
00347   Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
00348   
00349   // Pardiso cannot solve in-place
00350   if(rhs_ptr == x.derived().data())
00351   {
00352     tmp = b;
00353     rhs_ptr = tmp.data();
00354   }
00355   
00356   Index error;
00357   error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size),
00358                                                             m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
00359                                                             m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
00360                                                             rhs_ptr, x.derived().data());
00361 
00362   manageErrorCode(error);
00363 }
00364 
00365 
00383 template<typename MatrixType>
00384 class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
00385 {
00386   protected:
00387     typedef PardisoImpl<PardisoLU> Base;
00388     typedef typename Base::Scalar Scalar;
00389     typedef typename Base::RealScalar RealScalar;
00390     using Base::pardisoInit;
00391     using Base::m_matrix;
00392     friend class PardisoImpl< PardisoLU<MatrixType> >;
00393 
00394   public:
00395 
00396     using Base::compute;
00397     using Base::solve;
00398 
00399     PardisoLU()
00400       : Base()
00401     {
00402       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00403     }
00404 
00405     explicit PardisoLU(const MatrixType& matrix)
00406       : Base()
00407     {
00408       pardisoInit(Base::ScalarIsComplex ? 13 : 11);
00409       compute(matrix);
00410     }
00411   protected:
00412     void getMatrix(const MatrixType& matrix)
00413     {
00414       m_matrix = matrix;
00415       m_matrix.makeCompressed();
00416     }
00417 };
00418 
00438 template<typename MatrixType, int _UpLo>
00439 class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
00440 {
00441   protected:
00442     typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
00443     typedef typename Base::Scalar Scalar;
00444     typedef typename Base::RealScalar RealScalar;
00445     using Base::pardisoInit;
00446     using Base::m_matrix;
00447     friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
00448 
00449   public:
00450 
00451     typedef typename Base::StorageIndex StorageIndex;
00452     enum { UpLo = _UpLo };
00453     using Base::compute;
00454 
00455     PardisoLLT()
00456       : Base()
00457     {
00458       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00459     }
00460 
00461     explicit PardisoLLT(const MatrixType& matrix)
00462       : Base()
00463     {
00464       pardisoInit(Base::ScalarIsComplex ? 4 : 2);
00465       compute(matrix);
00466     }
00467     
00468   protected:
00469     
00470     void getMatrix(const MatrixType& matrix)
00471     {
00472       // PARDISO supports only upper, row-major matrices
00473       PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
00474       m_matrix.resize(matrix.rows(), matrix.cols());
00475       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00476       m_matrix.makeCompressed();
00477     }
00478 };
00479 
00501 template<typename MatrixType, int Options>
00502 class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
00503 {
00504   protected:
00505     typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
00506     typedef typename Base::Scalar Scalar;
00507     typedef typename Base::RealScalar RealScalar;
00508     using Base::pardisoInit;
00509     using Base::m_matrix;
00510     friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
00511 
00512   public:
00513 
00514     typedef typename Base::StorageIndex StorageIndex;
00515     using Base::compute;
00516     enum { UpLo = Options&(Upper|Lower) };
00517 
00518     PardisoLDLT()
00519       : Base()
00520     {
00521       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00522     }
00523 
00524     explicit PardisoLDLT(const MatrixType& matrix)
00525       : Base()
00526     {
00527       pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
00528       compute(matrix);
00529     }
00530     
00531     void getMatrix(const MatrixType& matrix)
00532     {
00533       // PARDISO supports only upper, row-major matrices
00534       PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
00535       m_matrix.resize(matrix.rows(), matrix.cols());
00536       m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
00537       m_matrix.makeCompressed();
00538     }
00539 };
00540 
00541 } // end namespace Eigen
00542 
00543 #endif // EIGEN_PARDISOSUPPORT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends