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();
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)");
361 if( parameterList->isParameter(
"IsContiguous") ){
362 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
387 template <
class Matrix,
class Vector>
388 Teuchos::RCP<const Teuchos::ParameterList>
394 using Teuchos::tuple;
395 using Teuchos::toString;
396 using Teuchos::EnhancedNumberValidator;
397 using Teuchos::setStringToIntegralParameter;
398 using Teuchos::anyNumberParameterEntryValidator;
400 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
402 if( is_null(valid_params) ){
403 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
407 PMKL::_INTEGER_t mtype_temp = mtype_;
408 PMKL::_INTEGER_t iparm_temp[64];
409 PMKL::pardisoinit(pt_dummy,
410 const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
411 const_cast<PMKL::_INTEGER_t*>(iparm_temp));
413 setStringToIntegralParameter<int>(
"IPARM(2)", toString(iparm_temp[1]),
414 "Fill-in reducing ordering for the input matrix",
415 tuple<string>(
"0",
"2",
"3"),
416 tuple<string>(
"The minimum degree algorithm",
417 "Nested dissection algorithm from METIS",
418 "OpenMP parallel nested dissection algorithm"),
422 Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
423 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
424 iparm_4_validator->setMin(0);
425 pl->set(
"IPARM(4)" , as<int>(iparm_temp[3]) ,
"Preconditioned CGS/CG",
428 setStringToIntegralParameter<int>(
"IPARM(12)", toString(iparm_temp[11]),
429 "Solve with transposed or conjugate transposed matrix A",
430 tuple<string>(
"0",
"1",
"2"),
431 tuple<string>(
"Non-transposed",
432 "Conjugate-transposed",
437 setStringToIntegralParameter<int>(
"IPARM(24)", toString(iparm_temp[23]),
438 "Parallel factorization control",
439 tuple<string>(
"0",
"1"),
440 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
441 "PARDISO uses the new two-level factorization algorithm"),
445 setStringToIntegralParameter<int>(
"IPARM(25)", toString(iparm_temp[24]),
446 "Parallel forward/backward solve control",
447 tuple<string>(
"0",
"1"),
448 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
449 "PARDISO uses the sequential forward and backward solve"),
453 setStringToIntegralParameter<int>(
"IPARM(60)", toString(iparm_temp[59]),
454 "PARDISO mode (OOC mode)",
455 tuple<string>(
"0",
"2"),
456 tuple<string>(
"In-core PARDISO",
457 "Out-of-core PARDISO. The OOC PARDISO can solve very "
458 "large problems by holding the matrix factors in files "
459 "on the disk. Hence the amount of RAM required by OOC "
460 "PARDISO is significantly reduced."),
464 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
465 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
467 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
468 accept_int.allowInt(
true );
470 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
471 anyNumberParameterEntryValidator(preferred_int, accept_int));
473 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
474 anyNumberParameterEntryValidator(preferred_int, accept_int));
476 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
477 anyNumberParameterEntryValidator(preferred_int, accept_int));
479 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
489 template <
class Matrix,
class Vector>
493 #ifdef HAVE_AMESOS2_TIMERS
494 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
498 if( current_phase == PREORDERING )
return(
false );
501 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
502 Kokkos::resize(colind_view_, this->globalNumNonZeros_);
503 Kokkos::resize(rowptr_view_, this->globalNumRows_ + 1);
508 #ifdef HAVE_AMESOS2_TIMERS
509 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
514 host_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
515 this->matrixA_.ptr(),
516 nzvals_view_, colind_view_, rowptr_view_, nnz_ret,
517 (is_contiguous_ ==
true) ? ROOTED : CONTIGUOUS_AND_ROOTED,
519 this->rowIndexBase_);
526 template <
class Matrix,
class Vector>
532 Teuchos::broadcast(*(this->getComm()), 0, &error_i);
534 if( error == 0 )
return;
536 std::string errmsg =
"Other error";
539 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
542 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
545 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
549 "PardisoMKL reported error: 'Zero pivot, numerical "
550 "factorization or iterative refinement problem'";
553 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
556 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
559 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
562 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
565 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
568 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
571 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
575 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
579 template <
class Matrix,
class Vector>
592 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
593 std::invalid_argument,
594 "Cannot set a real Pardiso matrix type with scalar type complex" );
597 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
598 std::invalid_argument,
599 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
602 TEUCHOS_TEST_FOR_EXCEPTION(
true,
603 std::invalid_argument,
604 "Symmetric matrices are not yet supported by the Amesos2 interface" );
610 template <
class Matrix,
class Vector>
613 template <
class Matrix,
class Vector>
614 const typename PardisoMKL<Matrix,Vector>::int_t
617 template <
class Matrix,
class Vector>
618 const typename PardisoMKL<Matrix,Vector>::int_t
621 template <
class Matrix,
class Vector>
622 const typename PardisoMKL<Matrix,Vector>::int_t
628 #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:581
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:491
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:528
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:389
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