53 #ifndef AMESOS2_PARDISOMKL_DEF_HPP 54 #define AMESOS2_PARDISOMKL_DEF_HPP 58 #include <Teuchos_Tuple.hpp> 59 #include <Teuchos_toString.hpp> 60 #include <Teuchos_StandardParameterEntryValidators.hpp> 70 # include <mkl_pardiso.h> 73 template <
class Matrix,
class Vector>
75 Teuchos::RCP<Vector> X,
76 Teuchos::RCP<const Vector> B)
81 , n_(
Teuchos::as<int_t>(this->globalNumRows_))
82 , perm_(this->globalNumRows_)
88 PMKL::_INTEGER_t iparm_temp[64];
89 PMKL::_INTEGER_t mtype_temp =
mtype_;
90 PMKL::pardisoinit(
pt_, &mtype_temp, iparm_temp);
92 for(
int i = 0; i < 64; ++i ){
97 if( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
105 #ifdef HAVE_AMESOS2_DEBUG 111 template <
class Matrix,
class Vector>
118 void *bdummy, *xdummy;
122 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
123 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
124 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
125 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
126 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
129 check_pardiso_mkl_error(Amesos2::CLEAN, error);
133 template<
class Matrix,
class Vector>
144 template <
class Matrix,
class Vector>
151 #ifdef HAVE_AMESOS2_TIMERS 152 Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
156 void *bdummy, *xdummy;
158 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
159 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
160 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
161 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
162 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
165 check_pardiso_mkl_error(Amesos2::SYMBFACT, error);
170 this->setNnzLU(iparm_[17]);
176 template <
class Matrix,
class Vector>
183 #ifdef HAVE_AMESOS2_TIMERS 184 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
188 void *bdummy, *xdummy;
190 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
191 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
192 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
193 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
194 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
197 check_pardiso_mkl_error(Amesos2::NUMFACT, error);
203 template <
class Matrix,
class Vector>
213 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
214 nrhs_ = as<int_t>(X->getGlobalNumVectors());
216 const size_t val_store_size = as<size_t>(ld_rhs * nrhs_);
217 xvals_.resize(val_store_size);
218 bvals_.resize(val_store_size);
221 #ifdef HAVE_AMESOS2_TIMERS 222 Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
223 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
227 solver_scalar_type>::do_get(B, bvals_(),
233 #ifdef HAVE_AMESOS2_TIMERS 234 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
237 const int_t phase = 33;
239 function_map::pardiso( pt_,
240 const_cast<int_t*>(&maxfct_),
241 const_cast<int_t*>(&mnum_),
242 const_cast<int_t*>(&mtype_),
243 const_cast<int_t*>(&phase),
244 const_cast<int_t*>(&n_),
245 const_cast<solver_scalar_type*>(nzvals_.getRawPtr()),
246 const_cast<int_t*>(rowptr_.getRawPtr()),
247 const_cast<int_t*>(colind_.getRawPtr()),
248 const_cast<int_t*>(perm_.getRawPtr()),
250 const_cast<int_t*>(iparm_),
251 const_cast<int_t*
>(&msglvl_),
252 as<void*>(bvals_.getRawPtr()),
253 as<void*>(xvals_.getRawPtr()), &error );
256 check_pardiso_mkl_error(Amesos2::SOLVE, error);
260 #ifdef HAVE_AMESOS2_TIMERS 261 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
266 solver_scalar_type>::do_put(X, xvals_(),
275 template <
class Matrix,
class Vector>
280 return( this->globalNumRows_ == this->globalNumCols_ );
284 template <
class Matrix,
class Vector>
289 using Teuchos::getIntegralValue;
290 using Teuchos::ParameterEntryValidator;
292 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
294 if( parameterList->isParameter(
"IPARM(2)") )
296 RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry(
"IPARM(2)").validator();
297 parameterList->getEntry(
"IPARM(2)").setValidator(fillin_validator);
298 iparm_[1] = getIntegralValue<int>(*parameterList,
"IPARM(2)");
301 if( parameterList->isParameter(
"IPARM(4)") )
303 RCP<const ParameterEntryValidator> prec_validator = valid_params->getEntry(
"IPARM(4)").validator();
304 parameterList->getEntry(
"IPARM(4)").setValidator(prec_validator);
305 iparm_[3] = getIntegralValue<int>(*parameterList,
"IPARM(4)");
308 if( parameterList->isParameter(
"IPARM(8)") )
310 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IPARM(8)").validator();
311 parameterList->getEntry(
"IPARM(8)").setValidator(refine_validator);
312 iparm_[7] = getIntegralValue<int>(*parameterList,
"IPARM(8)");
315 if( parameterList->isParameter(
"IPARM(10)") )
317 RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry(
"IPARM(10)").validator();
318 parameterList->getEntry(
"IPARM(10)").setValidator(pivot_perturb_validator);
319 iparm_[9] = getIntegralValue<int>(*parameterList,
"IPARM(10)");
324 iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
326 if( parameterList->isParameter(
"IPARM(12)") )
328 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"IPARM(12)").validator();
329 parameterList->getEntry(
"IPARM(12)").setValidator(trans_validator);
330 iparm_[11] = getIntegralValue<int>(*parameterList,
"IPARM(12)");
333 if( parameterList->isParameter(
"IPARM(18)") )
335 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry(
"IPARM(18)").validator();
336 parameterList->getEntry(
"IPARM(18)").setValidator(report_validator);
337 iparm_[17] = getIntegralValue<int>(*parameterList,
"IPARM(18)");
340 if( parameterList->isParameter(
"IPARM(24)") )
342 RCP<const ParameterEntryValidator> par_fact_validator = valid_params->getEntry(
"IPARM(24)").validator();
343 parameterList->getEntry(
"IPARM(24)").setValidator(par_fact_validator);
344 iparm_[23] = getIntegralValue<int>(*parameterList,
"IPARM(24)");
347 if( parameterList->isParameter(
"IPARM(25)") )
349 RCP<const ParameterEntryValidator> par_fbsolve_validator = valid_params->getEntry(
"IPARM(25)").validator();
350 parameterList->getEntry(
"IPARM(25)").setValidator(par_fbsolve_validator);
351 iparm_[24] = getIntegralValue<int>(*parameterList,
"IPARM(25)");
354 if( parameterList->isParameter(
"IPARM(60)") )
356 RCP<const ParameterEntryValidator> ooc_validator = valid_params->getEntry(
"IPARM(60)").validator();
357 parameterList->getEntry(
"IPARM(60)").setValidator(ooc_validator);
358 iparm_[59] = getIntegralValue<int>(*parameterList,
"IPARM(60)");
383 template <
class Matrix,
class Vector>
384 Teuchos::RCP<const Teuchos::ParameterList>
390 using Teuchos::tuple;
391 using Teuchos::toString;
392 using Teuchos::EnhancedNumberValidator;
393 using Teuchos::setStringToIntegralParameter;
394 using Teuchos::anyNumberParameterEntryValidator;
396 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
398 if( is_null(valid_params) ){
399 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
403 PMKL::_INTEGER_t mtype_temp = mtype_;
404 PMKL::_INTEGER_t iparm_temp[64];
405 PMKL::pardisoinit(pt_dummy,
406 const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
407 const_cast<PMKL::_INTEGER_t*>(iparm_temp));
409 setStringToIntegralParameter<int>(
"IPARM(2)", toString(iparm_temp[1]),
410 "Fill-in reducing ordering for the input matrix",
411 tuple<string>(
"0",
"2",
"3"),
412 tuple<string>(
"The minimum degree algorithm",
413 "Nested dissection algorithm from METIS",
414 "OpenMP parallel nested dissection algorithm"),
418 Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
419 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
420 iparm_4_validator->setMin(0);
421 pl->set(
"IPARM(4)" , as<int>(iparm_temp[3]) ,
"Preconditioned CGS/CG",
424 setStringToIntegralParameter<int>(
"IPARM(12)", toString(iparm_temp[11]),
425 "Solve with transposed or conjugate transposed matrix A",
426 tuple<string>(
"0",
"1",
"2"),
427 tuple<string>(
"Non-transposed",
428 "Conjugate-transposed",
433 setStringToIntegralParameter<int>(
"IPARM(24)", toString(iparm_temp[23]),
434 "Parallel factorization control",
435 tuple<string>(
"0",
"1"),
436 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
437 "PARDISO uses the new two-level factorization algorithm"),
441 setStringToIntegralParameter<int>(
"IPARM(25)", toString(iparm_temp[24]),
442 "Parallel forward/backward solve control",
443 tuple<string>(
"0",
"1"),
444 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
445 "PARDISO uses the sequential forward and backward solve"),
449 setStringToIntegralParameter<int>(
"IPARM(60)", toString(iparm_temp[59]),
450 "PARDISO mode (OOC mode)",
451 tuple<string>(
"0",
"2"),
452 tuple<string>(
"In-core PARDISO",
453 "Out-of-core PARDISO. The OOC PARDISO can solve very " 454 "large problems by holding the matrix factors in files " 455 "on the disk. Hence the amount of RAM required by OOC " 456 "PARDISO is significantly reduced."),
460 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
461 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
463 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
464 accept_int.allowInt(
true );
466 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
467 anyNumberParameterEntryValidator(preferred_int, accept_int));
469 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
470 anyNumberParameterEntryValidator(preferred_int, accept_int));
472 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
473 anyNumberParameterEntryValidator(preferred_int, accept_int));
483 template <
class Matrix,
class Vector>
487 #ifdef HAVE_AMESOS2_TIMERS 488 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
492 if( current_phase == PREORDERING )
return(
false );
495 nzvals_.resize(this->globalNumNonZeros_);
496 colind_.resize(this->globalNumNonZeros_);
497 rowptr_.resize(this->globalNumRows_ + 1);
502 #ifdef HAVE_AMESOS2_TIMERS 503 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
509 int_t,int_t>::do_get(this->matrixA_.ptr(),
510 nzvals_(), colind_(), rowptr_(),
518 template <
class Matrix,
class Vector>
524 Teuchos::broadcast(*(this->getComm()), 0, &error_i);
526 if( error == 0 )
return;
528 std::string errmsg =
"Other error";
531 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
534 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
537 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
541 "PardisoMKL reported error: 'Zero pivot, numerical " 542 "factorization or iterative refinement problem'";
545 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
548 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
551 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
554 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
557 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
560 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
563 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
567 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
571 template <
class Matrix,
class Vector>
584 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
585 std::invalid_argument,
586 "Cannot set a real Pardiso matrix type with scalar type complex" );
589 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
590 std::invalid_argument,
591 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
594 TEUCHOS_TEST_FOR_EXCEPTION(
true,
595 std::invalid_argument,
596 "Symmetric matrices are not yet supported by the Amesos2 interface" );
602 template <
class Matrix,
class Vector>
605 template <
class Matrix,
class Vector>
606 const typename PardisoMKL<Matrix,Vector>::int_t
609 template <
class Matrix,
class Vector>
610 const typename PardisoMKL<Matrix,Vector>::int_t
613 template <
class Matrix,
class Vector>
614 const typename PardisoMKL<Matrix,Vector>::int_t
620 #endif // AMESOS2_PARDISOMKL_DEF_HPP Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Definition: Amesos2_TypeDecl.hpp:141
~PardisoMKL()
Destructor.
Definition: Amesos2_PardisoMKL_def.hpp:112
PardisoMKL(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_PardisoMKL_def.hpp:74
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:591
void check_pardiso_mkl_error(EPhase phase, int_t error) const
Throws an appropriate runtime error in the event that error < 0 .
Definition: Amesos2_PardisoMKL_def.hpp:520
void set_pardiso_mkl_matrix_type(int_t mtype=0)
Definition: Amesos2_PardisoMKL_def.hpp:573
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
PardisoMKL specific solve.
Definition: Amesos2_PardisoMKL_def.hpp:205
void * pt_[64]
PardisoMKL internal data address pointer.
Definition: Amesos2_PardisoMKL_decl.hpp:285
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A template class that does nothing useful besides show developers what, in general, needs to be done to add a new solver interface to the Amesos2 collection.
Amesos2 interface to the PardisoMKL package.
Definition: Amesos2_PardisoMKL_decl.hpp:83
int numericFactorization_impl()
PardisoMKL specific numeric factorization.
Definition: Amesos2_PardisoMKL_def.hpp:178
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_PardisoMKL_def.hpp:277
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using PardisoMKL.
Definition: Amesos2_PardisoMKL_def.hpp:146
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_PardisoMKL_def.hpp:485
Definition: Amesos2_Cholmod_TypeMap.hpp:92
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_PardisoMKL_def.hpp:286
Definition: Amesos2_TypeDecl.hpp:127
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_PardisoMKL_def.hpp:385
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
int_t mtype_
The matrix type. We deal only with unsymmetrix matrices.
Definition: Amesos2_PardisoMKL_decl.hpp:287
int_t iparm_[64]
Definition: Amesos2_PardisoMKL_decl.hpp:297
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_PardisoMKL_def.hpp:135