Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_PardisoMKL_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package 
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00044 
00053 #ifndef AMESOS2_PARDISOMKL_DEF_HPP
00054 #define AMESOS2_PARDISOMKL_DEF_HPP
00055 
00056 #include <map>
00057 
00058 #include <Teuchos_Tuple.hpp>
00059 #include <Teuchos_toString.hpp>
00060 #include <Teuchos_StandardParameterEntryValidators.hpp>
00061 
00062 #include "Amesos2_SolverCore_def.hpp"
00063 #include "Amesos2_PardisoMKL_decl.hpp"
00064 
00065 
00066 namespace Amesos2 {
00067 
00068   namespace PMKL {
00069 #   include <mkl.h>
00070 #   include <mkl_pardiso.h>
00071   }
00072 
00073   template <class Matrix, class Vector>
00074   PardisoMKL<Matrix,Vector>::PardisoMKL(Teuchos::RCP<const Matrix> A,
00075                                         Teuchos::RCP<Vector>       X,
00076                                         Teuchos::RCP<const Vector> B)
00077     : SolverCore<Amesos2::PardisoMKL,Matrix,Vector>(A, X, B) // instantiate superclass
00078     , nzvals_()
00079     , colind_()
00080     , rowptr_()
00081     , n_(Teuchos::as<int_t>(this->globalNumRows_))
00082     , perm_(this->globalNumRows_)
00083     , nrhs_(0)
00084   {
00085     // set the default matrix type
00086     set_pardiso_mkl_matrix_type();
00087 
00088     PMKL::_INTEGER_t iparm_temp[64];
00089     PMKL::_INTEGER_t mtype_temp = mtype_;
00090     PMKL::pardisoinit(pt_, &mtype_temp, iparm_temp);
00091 
00092     for( int i = 0; i < 64; ++i ){
00093       iparm_[i] = iparm_temp[i];
00094     }
00095 
00096     // set single or double precision
00097     if( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
00098       iparm_[27] = 1;           // single-precision
00099     } else {
00100       iparm_[27] = 0;           // double-precision
00101     }
00102 
00103     // Reset some of the default parameters
00104     iparm_[34] = 1;             // Use zero-based indexing
00105 #ifdef HAVE_AMESOS2_DEBUG
00106     iparm_[26] = 1;             // turn the Pardiso matrix checker on
00107 #endif
00108   }
00109 
00110 
00111   template <class Matrix, class Vector>
00112   PardisoMKL<Matrix,Vector>::~PardisoMKL( )
00113   {
00114     /*
00115      * Free any memory allocated by the PardisoMKL library functions
00116      */
00117     int_t error = 0;
00118     void *bdummy, *xdummy;
00119 
00120     if( this->root_ ){
00121       int_t phase = -1;         // release all internal solver memory
00122       function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
00123                              const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
00124                              nzvals_.getRawPtr(), rowptr_.getRawPtr(),
00125                              colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
00126                              const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &error );
00127     }
00128 
00129     check_pardiso_mkl_error(Amesos2::CLEAN, error);
00130   }
00131 
00132 
00133   template<class Matrix, class Vector>
00134   int
00135   PardisoMKL<Matrix,Vector>::preOrdering_impl()
00136   {
00137     // preOrdering done in PardisoMKL during "Analysis" (aka symbolic
00138     // factorization) phase
00139 
00140     return(0);
00141   }
00142 
00143 
00144   template <class Matrix, class Vector>
00145   int
00146   PardisoMKL<Matrix,Vector>::symbolicFactorization_impl()
00147   {
00148     int_t error = 0;
00149 
00150     if( this->root_ ){
00151 #ifdef HAVE_AMESOS2_TIMERS
00152       Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
00153 #endif
00154 
00155       int_t phase = 11;
00156       void *bdummy, *xdummy;
00157 
00158       function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
00159                              const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
00160                              nzvals_.getRawPtr(), rowptr_.getRawPtr(),
00161                              colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
00162                              const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &error );
00163     }
00164 
00165     check_pardiso_mkl_error(Amesos2::SYMBFACT, error);
00166 
00167     // Pardiso only lets you retrieve the total number of factor
00168     // non-zeros, not for each individually.  We should document how
00169     // such a situation is reported.
00170     this->setNnzLU(iparm_[17]);
00171 
00172     return(0);
00173   }
00174 
00175 
00176   template <class Matrix, class Vector>
00177   int
00178   PardisoMKL<Matrix,Vector>::numericFactorization_impl()
00179   {
00180     int_t error = 0;
00181 
00182     if( this->root_ ){
00183 #ifdef HAVE_AMESOS2_TIMERS
00184       Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
00185 #endif
00186 
00187       int_t phase = 22;
00188       void *bdummy, *xdummy;
00189 
00190       function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
00191                              const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
00192                              nzvals_.getRawPtr(), rowptr_.getRawPtr(),
00193                              colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
00194                              const_cast<int_t*>(&msglvl_), &bdummy, &xdummy, &error );
00195     }
00196 
00197     check_pardiso_mkl_error(Amesos2::NUMFACT, error);
00198 
00199     return( 0 );
00200   }
00201 
00202 
00203   template <class Matrix, class Vector>
00204   int
00205   PardisoMKL<Matrix,Vector>::solve_impl(const Teuchos::Ptr<MultiVecAdapter<Vector> >       X,
00206                                         const Teuchos::Ptr<const MultiVecAdapter<Vector> > B) const
00207   {
00208     using Teuchos::as;
00209 
00210     int_t error = 0;
00211 
00212     // Get B data
00213     const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
00214     nrhs_ = as<int_t>(X->getGlobalNumVectors());
00215 
00216     const size_t val_store_size = as<size_t>(ld_rhs * nrhs_);
00217     xvals_.resize(val_store_size);
00218     bvals_.resize(val_store_size);
00219 
00220     {                             // Get values from RHS B
00221 #ifdef HAVE_AMESOS2_TIMERS
00222       Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
00223       Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
00224 #endif
00225       Util::get_1d_copy_helper<
00226         MultiVecAdapter<Vector>,
00227         solver_scalar_type>::do_get(B, bvals_(),
00228                                     as<size_t>(ld_rhs),
00229                                     ROOTED);
00230     }
00231 
00232     if( this->root_ ){
00233 #ifdef HAVE_AMESOS2_TIMERS
00234       Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
00235 #endif
00236 
00237       const int_t phase = 33;
00238 
00239       function_map::pardiso( pt_,
00240                              const_cast<int_t*>(&maxfct_),
00241                              const_cast<int_t*>(&mnum_),
00242                              const_cast<int_t*>(&mtype_),
00243                              const_cast<int_t*>(&phase),
00244                              const_cast<int_t*>(&n_),
00245                              const_cast<solver_scalar_type*>(nzvals_.getRawPtr()),
00246                              const_cast<int_t*>(rowptr_.getRawPtr()),
00247                              const_cast<int_t*>(colind_.getRawPtr()),
00248                              const_cast<int_t*>(perm_.getRawPtr()),
00249                              &nrhs_,
00250                              const_cast<int_t*>(iparm_),
00251                              const_cast<int_t*>(&msglvl_),
00252                              as<void*>(bvals_.getRawPtr()),
00253                              as<void*>(xvals_.getRawPtr()), &error );
00254     }
00255 
00256     check_pardiso_mkl_error(Amesos2::SOLVE, error);
00257 
00258     /* Export X from root to the global space */
00259     {
00260 #ifdef HAVE_AMESOS2_TIMERS
00261       Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
00262 #endif
00263 
00264       Util::put_1d_data_helper<
00265       MultiVecAdapter<Vector>,
00266         solver_scalar_type>::do_put(X, xvals_(),
00267                                     as<size_t>(ld_rhs),
00268                                     ROOTED);
00269   }
00270 
00271     return( 0 );
00272 }
00273 
00274 
00275   template <class Matrix, class Vector>
00276   bool
00277   PardisoMKL<Matrix,Vector>::matrixShapeOK_impl() const
00278   {
00279     // PardisoMKL supports square matrices
00280     return( this->globalNumRows_ == this->globalNumCols_ );
00281   }
00282 
00283 
00284 template <class Matrix, class Vector>
00285 void
00286 PardisoMKL<Matrix,Vector>::setParameters_impl(const Teuchos::RCP<Teuchos::ParameterList> & parameterList )
00287 {
00288   iparm_[1]  = validators[2]->getIntegralValue(*parameterList, "IPARM(2)",
00289                  validators[2]->getDefaultParameterName());
00290   iparm_[3]  = parameterList->get<int>("IPARM(4)" , (int)iparm_[3]);
00291   iparm_[7]  = parameterList->get<int>("IPARM(8)" , (int)iparm_[7]);
00292   iparm_[9]  = parameterList->get<int>("IPARM(10)", (int)iparm_[9]);
00293   iparm_[17] = parameterList->get<int>("IPARM(18)", (int)iparm_[17]);
00294   iparm_[23] = validators[24]->getIntegralValue(*parameterList, "IPARM(24)",
00295             validators[24]->getDefaultParameterName());
00296   iparm_[24] = validators[25]->getIntegralValue(*parameterList, "IPARM(25)",
00297             validators[25]->getDefaultParameterName());
00298   iparm_[59] = validators[60]->getIntegralValue(*parameterList, "IPARM(60)",
00299             validators[60]->getDefaultParameterName());
00300 }
00301 
00302 
00303 /*
00304  * TODO: It would be nice if the parameters could be expressed as
00305  * either all string or as all integers.  I see no way of doing this
00306  * at present with the standard validators.  However, we could create
00307  * our own validators or kindly ask the Teuchos team to add some
00308  * features for use.
00309  *
00310  * The issue is that with the current validators we cannot specify
00311  * arbitrary sets of numbers that are the only allowed parameters.
00312  * For example the IPARM(2) parameter can take only the values 0, 2,
00313  * and 3.  The EnhancedNumberValidator can take a min value, and max
00314  * value, and a step size, but with those options there is no way to
00315  * specify the needed set.
00316  *
00317  * Another missing feature is the ability to give docstrings for such
00318  * numbers.  For example IPARM(25) can take on the values 0 and 1.
00319  * This would be easy enough to accomplish with just a number
00320  * validator, but then have no way to document the effect of each
00321  * value.
00322  */
00323 template <class Matrix, class Vector>
00324 Teuchos::RCP<const Teuchos::ParameterList>
00325 PardisoMKL<Matrix,Vector>::getValidParameters_impl() const
00326 {
00327   using std::string;
00328   using Teuchos::as;
00329   using Teuchos::RCP;
00330   using Teuchos::tuple;
00331   using Teuchos::toString;
00332   using Teuchos::EnhancedNumberValidator;
00333   using Teuchos::setStringToIntegralParameter;
00334   using Teuchos::anyNumberParameterEntryValidator;
00335   using Teuchos::stringToIntegralParameterEntryValidator;
00336   typedef Teuchos::StringToIntegralParameterEntryValidator<int> STIPEV;
00337   Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
00338     Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
00339 
00340   static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
00341 
00342   if( is_null(valid_params) ){
00343     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00344 
00345     // Use pardisoinit to get some default values;
00346     void *pt_dummy[64];
00347     PMKL::_INTEGER_t mtype_temp = mtype_;
00348     PMKL::_INTEGER_t iparm_temp[64];
00349     PMKL::pardisoinit(pt_dummy,
00350                       const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
00351                       const_cast<PMKL::_INTEGER_t*>(iparm_temp));
00352 
00353     // Initialize our parameter validators, saving the string to int validators for later
00354     RCP<STIPEV> iparm_2_validator
00355       = stringToIntegralParameterEntryValidator<int>(tuple<string>("0", "2", "3"),
00356                  tuple<string>("The minimum degree algorithm",
00357                    "Nested dissection algorithm from METIS",
00358                    "OpenMP parallel nested dissection algorithm"),
00359                  tuple<int>(0, 2, 3),
00360                  toString(iparm_temp[1]));
00361     validators.insert( std::pair<int,RCP<STIPEV> >(2, iparm_2_validator) );
00362     
00363     Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
00364       = Teuchos::rcp( new EnhancedNumberValidator<int>() );
00365     iparm_4_validator->setMin(0);
00366 
00367     RCP<STIPEV> iparm_24_validator
00368       = stringToIntegralParameterEntryValidator<int>(tuple<string>("0", "1"),
00369                  tuple<string>("PARDISO uses the previous algorithm for factorization",
00370                    "PARDISO uses the new two-level factorization algorithm"),
00371                  tuple<int>(0, 1),
00372                  toString(iparm_temp[23]));
00373     validators.insert( std::pair<int,RCP<STIPEV> >(24, iparm_24_validator) );
00374 
00375     RCP<STIPEV> iparm_25_validator
00376       = stringToIntegralParameterEntryValidator<int>(tuple<string>("0", "1"),
00377                  tuple<string>("PARDISO uses the parallel algorithm for the solve step",
00378                    "PARDISO uses the sequential forward and backward solve"),
00379                  tuple<int>(0, 1),
00380                  toString(iparm_temp[24]));
00381     validators.insert( std::pair<int,RCP<STIPEV> >(25, iparm_25_validator) );
00382 
00383     RCP<STIPEV> iparm_60_validator
00384       = stringToIntegralParameterEntryValidator<int>(tuple<string>("0", "2"),
00385                  tuple<string>("In-core PARDISO",
00386                    "Out-of-core PARDISO.  The OOC PARDISO can solve very "
00387                    "large problems by holding the matrix factors in files "
00388                    "on the disk. Hence the amount of RAM required by OOC "
00389                    "PARDISO is significantly reduced."),
00390                  tuple<int>(0, 2),
00391                  toString(iparm_temp[59]));
00392     validators.insert( std::pair<int,RCP<STIPEV> >(60, iparm_60_validator) );
00393 
00394     Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int( false );
00395     accept_int.allowInt( true );
00396 
00397     pl->set("IPARM(2)" , validators[2]->getDefaultParameterName(),
00398       "Fill-in reducing ordering for the input matrix", validators[2]);
00399 
00400     pl->set("IPARM(4)" , as<int>(iparm_temp[3]) , "Preconditioned CGS/CG",
00401             iparm_4_validator);
00402 
00403     pl->set("IPARM(8)" , as<int>(iparm_temp[8]) , "Iterative refinement step",
00404             anyNumberParameterEntryValidator(preferred_int, accept_int));
00405 
00406     pl->set("IPARM(10)", as<int>(iparm_temp[9]) , "Pivoting perturbation",
00407             anyNumberParameterEntryValidator(preferred_int, accept_int));
00408 
00409     pl->set("IPARM(18)", as<int>(iparm_temp[17]), "Report the number of non-zero elements in the factors",
00410             anyNumberParameterEntryValidator(preferred_int, accept_int));
00411 
00412     pl->set("IPARM(24)", validators[24]->getDefaultParameterName(),
00413       "Parallel factorization control", validators[24]);
00414     
00415     pl->set("IPARM(25)", validators[25]->getDefaultParameterName(),
00416       "Parallel forward/backward solve control", validators[25]);
00417 
00418     pl->set("IPARM(60)", validators[60]->getDefaultParameterName(),
00419       "PARDISO mode (OOC mode)", validators[60]);
00420 
00421     valid_params = pl;
00422   }
00423 
00424   return valid_params;
00425 }
00426 
00427 
00428 
00429 template <class Matrix, class Vector>
00430 bool
00431 PardisoMKL<Matrix,Vector>::loadA_impl(EPhase current_phase)
00432 {
00433 #ifdef HAVE_AMESOS2_TIMERS
00434   Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
00435 #endif
00436 
00437   // PardisoMKL does not need matrix data in the pre-ordering phase
00438   if( current_phase == PREORDERING ) return( false );
00439 
00440   if( this->root_ ){
00441     nzvals_.resize(this->globalNumNonZeros_);
00442     colind_.resize(this->globalNumNonZeros_);
00443     rowptr_.resize(this->globalNumRows_ + 1);
00444   }
00445 
00446   int_t nnz_ret = 0;
00447   {
00448 #ifdef HAVE_AMESOS2_TIMERS
00449     Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
00450 #endif
00451 
00452     Util::get_crs_helper<
00453     MatrixAdapter<Matrix>,
00454       solver_scalar_type,
00455       int_t,int_t>::do_get(this->matrixA_.ptr(),
00456                            nzvals_(), colind_(), rowptr_(),
00457                            nnz_ret, ROOTED, SORTED_INDICES);
00458 }
00459 
00460   return( true );
00461 }
00462 
00463 
00464 template <class Matrix, class Vector>
00465 void
00466 PardisoMKL<Matrix,Vector>::check_pardiso_mkl_error(EPhase phase,
00467                                                    int_t error) const
00468 {
00469   int error_i = error;
00470   Teuchos::broadcast(*(this->getComm()), 0, &error_i); // We only care about root's value
00471 
00472   if( error == 0 ) return;      // No error
00473 
00474   std::string errmsg = "Other error";
00475   switch( error ){
00476   case -1:
00477     errmsg = "PardisoMKL reported error: 'Input inconsistent'";
00478     break;
00479   case -2:
00480     errmsg = "PardisoMKL reported error: 'Not enough memory'";
00481     break;
00482   case -3:
00483     errmsg = "PardisoMKL reported error: 'Reordering problem'";
00484     break;
00485   case -4:
00486     errmsg =
00487       "PardisoMKL reported error: 'Zero pivot, numerical "
00488       "factorization or iterative refinement problem'";
00489     break;
00490   case -5:
00491     errmsg = "PardisoMKL reported error: 'Unclassified (internal) error'";
00492     break;
00493   case -6:
00494     errmsg = "PardisoMKL reported error: 'Reordering failed'";
00495     break;
00496   case -7:
00497     errmsg = "PardisoMKL reported error: 'Diagonal matrix is singular'";
00498     break;
00499   case -8:
00500     errmsg = "PardisoMKL reported error: '32-bit integer overflow problem'";
00501     break;
00502   case -9:
00503     errmsg = "PardisoMKL reported error: 'Not enough memory for OOC'";
00504     break;
00505   case -10:
00506     errmsg = "PardisoMKL reported error: 'Problems with opening OOC temporary files'";
00507     break;
00508   case -11:
00509     errmsg = "PardisoMKL reported error: 'Read/write problem with OOC data file'";
00510     break;
00511   }
00512 
00513   TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error, errmsg );
00514 }
00515 
00516 
00517 template <class Matrix, class Vector>
00518 void
00519 PardisoMKL<Matrix,Vector>::set_pardiso_mkl_matrix_type(int_t mtype)
00520 {
00521   if( mtype == 0 ){
00522     if( complex_ ){
00523       mtype_ = 13;              // complex, unsymmetric
00524     } else {
00525       mtype_ = 11;              // real, unsymmetric
00526     }
00527   } else {
00528     switch( mtype ){
00529     case 11:
00530       TEUCHOS_TEST_FOR_EXCEPTION( complex_,
00531                           std::invalid_argument,
00532                           "Cannot set a real Pardiso matrix type with scalar type complex" );
00533       mtype_ = 11; break;
00534     case 13:
00535       TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
00536                           std::invalid_argument,
00537                           "Cannot set a complex Pardiso matrix type with non-complex scalars" );
00538       mtype_ = 13; break;
00539     default:
00540       TEUCHOS_TEST_FOR_EXCEPTION( true,
00541                           std::invalid_argument,
00542                           "Symmetric matrices are not yet supported by the Amesos2 interface" );
00543     }
00544   }
00545 }
00546 
00547 
00548 template <class Matrix, class Vector>
00549 const char* PardisoMKL<Matrix,Vector>::name = "PARDISOMKL";
00550 
00551 template <class Matrix, class Vector>
00552 const typename PardisoMKL<Matrix,Vector>::int_t
00553 PardisoMKL<Matrix,Vector>::msglvl_ = 0;
00554 
00555 template <class Matrix, class Vector>
00556 const typename PardisoMKL<Matrix,Vector>::int_t
00557 PardisoMKL<Matrix,Vector>::maxfct_ = 1;
00558 
00559 template <class Matrix, class Vector>
00560 const typename PardisoMKL<Matrix,Vector>::int_t
00561 PardisoMKL<Matrix,Vector>::mnum_ = 1;
00562 
00563 
00564 } // end namespace Amesos
00565 
00566 #endif  // AMESOS2_PARDISOMKL_DEF_HPP