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)
78 , n_(Teuchos::as<int_t>(this->globalNumRows_))
79 , perm_(this->globalNumRows_)
81 , is_contiguous_(true)
86 PMKL::_INTEGER_t iparm_temp[64];
87 PMKL::_INTEGER_t mtype_temp =
mtype_;
88 PMKL::pardisoinit(
pt_, &mtype_temp, iparm_temp);
90 for(
int i = 0; i < 64; ++i ){
95 if constexpr ( std::is_same_v<solver_magnitude_type, PMKL::_REAL_t> ) {
103 #ifdef HAVE_AMESOS2_DEBUG
109 template <
class Matrix,
class Vector>
116 void *bdummy, *xdummy;
120 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
121 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
122 nzvals_view_.data(), rowptr_view_.data(),
123 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
124 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
127 check_pardiso_mkl_error(Amesos2::CLEAN, error);
131 template<
class Matrix,
class Vector>
142 template <
class Matrix,
class Vector>
149 #ifdef HAVE_AMESOS2_TIMERS
150 Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
154 void *bdummy, *xdummy;
156 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
157 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
158 nzvals_view_.data(), rowptr_view_.data(),
159 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
160 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
163 check_pardiso_mkl_error(Amesos2::SYMBFACT, error);
168 this->setNnzLU(iparm_[17]);
174 template <
class Matrix,
class Vector>
181 #ifdef HAVE_AMESOS2_TIMERS
182 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
186 void *bdummy, *xdummy;
188 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
189 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
190 nzvals_view_.data(), rowptr_view_.data(),
191 colind_view_.data(), perm_.getRawPtr(), &nrhs_, iparm_,
192 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
195 check_pardiso_mkl_error(Amesos2::NUMFACT, error);
201 template <
class Matrix,
class Vector>
211 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
212 nrhs_ = as<int_t>(X->getGlobalNumVectors());
214 const size_t val_store_size = as<size_t>(ld_rhs * nrhs_);
215 xvals_.resize(val_store_size);
216 bvals_.resize(val_store_size);
219 #ifdef HAVE_AMESOS2_TIMERS
220 Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
221 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
226 solver_scalar_type>::do_get(B, bvals_(),
228 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
229 this->rowIndexBase_);
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_view_.data()),
246 const_cast<int_t*>(rowptr_view_.data()),
247 const_cast<int_t*>(colind_view_.data()),
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_(),
268 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED);
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();
295 if( parameterList->isParameter(
"IPARM(2)") )
297 RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry(
"IPARM(2)").validator();
298 parameterList->getEntry(
"IPARM(2)").setValidator(fillin_validator);
299 iparm_[1] = getIntegralValue<int>(*parameterList,
"IPARM(2)");
303 if( parameterList->isParameter(
"IPARM(4)") )
305 RCP<const ParameterEntryValidator> prec_validator = valid_params->getEntry(
"IPARM(4)").validator();
306 parameterList->getEntry(
"IPARM(4)").setValidator(prec_validator);
307 iparm_[3] = getIntegralValue<int>(*parameterList,
"IPARM(4)");
311 if( parameterList->isParameter(
"IPARM(8)") )
313 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IPARM(8)").validator();
314 parameterList->getEntry(
"IPARM(8)").setValidator(refine_validator);
315 iparm_[7] = getIntegralValue<int>(*parameterList,
"IPARM(8)");
319 if( parameterList->isParameter(
"IPARM(10)") )
321 RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry(
"IPARM(10)").validator();
322 parameterList->getEntry(
"IPARM(10)").setValidator(pivot_perturb_validator);
323 iparm_[9] = getIntegralValue<int>(*parameterList,
"IPARM(10)");
328 iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
331 if( parameterList->isParameter(
"IPARM(12)") )
333 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"IPARM(12)").validator();
334 parameterList->getEntry(
"IPARM(12)").setValidator(trans_validator);
335 iparm_[11] = getIntegralValue<int>(*parameterList,
"IPARM(12)");
339 if( parameterList->isParameter(
"IPARM(13)") )
341 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"IPARM(13)").validator();
342 parameterList->getEntry(
"IPARM(13)").setValidator(trans_validator);
343 iparm_[12] = getIntegralValue<int>(*parameterList,
"IPARM(13)");
347 if( parameterList->isParameter(
"IPARM(18)") )
349 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry(
"IPARM(18)").validator();
350 parameterList->getEntry(
"IPARM(18)").setValidator(report_validator);
351 iparm_[17] = getIntegralValue<int>(*parameterList,
"IPARM(18)");
355 if( parameterList->isParameter(
"IPARM(24)") )
357 RCP<const ParameterEntryValidator> par_fact_validator = valid_params->getEntry(
"IPARM(24)").validator();
358 parameterList->getEntry(
"IPARM(24)").setValidator(par_fact_validator);
359 iparm_[23] = getIntegralValue<int>(*parameterList,
"IPARM(24)");
363 if( parameterList->isParameter(
"IPARM(25)") )
365 RCP<const ParameterEntryValidator> par_fbsolve_validator = valid_params->getEntry(
"IPARM(25)").validator();
366 parameterList->getEntry(
"IPARM(25)").setValidator(par_fbsolve_validator);
367 iparm_[24] = getIntegralValue<int>(*parameterList,
"IPARM(25)");
371 if( parameterList->isParameter(
"IPARM(60)") )
373 RCP<const ParameterEntryValidator> ooc_validator = valid_params->getEntry(
"IPARM(60)").validator();
374 parameterList->getEntry(
"IPARM(60)").setValidator(ooc_validator);
375 iparm_[59] = getIntegralValue<int>(*parameterList,
"IPARM(60)");
378 if( parameterList->isParameter(
"IsContiguous") ){
379 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
404 template <
class Matrix,
class Vector>
405 Teuchos::RCP<const Teuchos::ParameterList>
411 using Teuchos::tuple;
412 using Teuchos::toString;
413 using Teuchos::EnhancedNumberValidator;
414 using Teuchos::setStringToIntegralParameter;
415 using Teuchos::anyNumberParameterEntryValidator;
417 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
419 if( is_null(valid_params) ){
420 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
424 PMKL::_INTEGER_t mtype_temp = mtype_;
425 PMKL::_INTEGER_t iparm_temp[64];
426 PMKL::pardisoinit(pt_dummy,
427 const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
428 const_cast<PMKL::_INTEGER_t*>(iparm_temp));
430 setStringToIntegralParameter<int>(
"IPARM(2)", toString(iparm_temp[1]),
431 "Fill-in reducing ordering for the input matrix",
432 tuple<string>(
"0",
"2",
"3"),
433 tuple<string>(
"The minimum degree algorithm",
434 "Nested dissection algorithm from METIS",
435 "OpenMP parallel nested dissection algorithm"),
439 Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
440 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
441 iparm_4_validator->setMin(0);
442 pl->set(
"IPARM(4)" , as<int>(iparm_temp[3]) ,
"Preconditioned CGS/CG",
445 setStringToIntegralParameter<int>(
"IPARM(12)", toString(iparm_temp[11]),
446 "Solve with transposed or conjugate transposed matrix A",
447 tuple<string>(
"0",
"1",
"2"),
448 tuple<string>(
"Non-transposed",
449 "Conjugate-transposed",
454 setStringToIntegralParameter<int>(
"IPARM(13)", toString(iparm_temp[12]),
455 "Use weighted matching",
456 tuple<string>(
"0",
"1"),
457 tuple<string>(
"No matching",
"Use matching"),
461 setStringToIntegralParameter<int>(
"IPARM(24)", toString(iparm_temp[23]),
462 "Parallel factorization control",
463 tuple<string>(
"0",
"1"),
464 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
465 "PARDISO uses the new two-level factorization algorithm"),
469 setStringToIntegralParameter<int>(
"IPARM(25)", toString(iparm_temp[24]),
470 "Parallel forward/backward solve control",
471 tuple<string>(
"0",
"1"),
472 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
473 "PARDISO uses the sequential forward and backward solve"),
477 setStringToIntegralParameter<int>(
"IPARM(60)", toString(iparm_temp[59]),
478 "PARDISO mode (OOC mode)",
479 tuple<string>(
"0",
"2"),
480 tuple<string>(
"In-core PARDISO",
481 "Out-of-core PARDISO. The OOC PARDISO can solve very "
482 "large problems by holding the matrix factors in files "
483 "on the disk. Hence the amount of RAM required by OOC "
484 "PARDISO is significantly reduced."),
488 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
489 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
491 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
492 accept_int.allowInt(
true );
494 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
495 anyNumberParameterEntryValidator(preferred_int, accept_int));
497 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
498 anyNumberParameterEntryValidator(preferred_int, accept_int));
500 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
501 anyNumberParameterEntryValidator(preferred_int, accept_int));
503 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
513 template <
class Matrix,
class Vector>
517 #ifdef HAVE_AMESOS2_TIMERS
518 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
522 if( current_phase == PREORDERING )
return(
false );
525 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
526 Kokkos::resize(colind_view_, this->globalNumNonZeros_);
527 Kokkos::resize(rowptr_view_, this->globalNumRows_ + 1);
532 #ifdef HAVE_AMESOS2_TIMERS
533 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
538 host_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
539 this->matrixA_.ptr(),
540 nzvals_view_, colind_view_, rowptr_view_, nnz_ret,
541 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
543 this->rowIndexBase_);
550 template <
class Matrix,
class Vector>
556 Teuchos::broadcast(*(this->getComm()), 0, &error_i);
558 if( error == 0 )
return;
560 std::string errmsg =
"Other error";
563 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
566 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
569 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
573 "PardisoMKL reported error: 'Zero pivot, numerical "
574 "factorization or iterative refinement problem'";
577 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
580 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
583 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
586 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
589 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
592 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
595 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
599 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
603 template <
class Matrix,
class Vector>
616 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
617 std::invalid_argument,
618 "Cannot set a real Pardiso matrix type with scalar type complex" );
621 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
622 std::invalid_argument,
623 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
626 TEUCHOS_TEST_FOR_EXCEPTION(
true,
627 std::invalid_argument,
628 "Symmetric matrices are not yet supported by the Amesos2 interface" );
634 template <
class Matrix,
class Vector>
637 template <
class Matrix,
class Vector>
638 const typename PardisoMKL<Matrix,Vector>::int_t
641 template <
class Matrix,
class Vector>
642 const typename PardisoMKL<Matrix,Vector>::int_t
645 template <
class Matrix,
class Vector>
646 const typename PardisoMKL<Matrix,Vector>::int_t
652 #endif // AMESOS2_PARDISOMKL_DEF_HPP
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
~PardisoMKL()
Destructor.
Definition: Amesos2_PardisoMKL_def.hpp:110
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
void set_pardiso_mkl_matrix_type(int_t mtype=0)
Definition: Amesos2_PardisoMKL_def.hpp:605
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:267
void * pt_[64]
PardisoMKL internal data address pointer.
Definition: Amesos2_PardisoMKL_decl.hpp:289
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:176
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using PardisoMKL.
Definition: Amesos2_PardisoMKL_def.hpp:144
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:515
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:659
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_PardisoMKL_def.hpp:286
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:203
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:552
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:373
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_PardisoMKL_def.hpp:406
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_PardisoMKL_def.hpp:277
int_t mtype_
The matrix type. We deal only with unsymmetrix matrices.
Definition: Amesos2_PardisoMKL_decl.hpp:291
int_t iparm_[64]
Definition: Amesos2_PardisoMKL_decl.hpp:301
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_PardisoMKL_def.hpp:133