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( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
103#ifdef HAVE_AMESOS2_DEBUG
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);
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]);
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);
204 const Teuchos::Ptr<
const MultiVecAdapter<Vector> > B)
const
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_ );
224 if ( is_contiguous_ ==
true ) {
226 MultiVecAdapter<Vector>,
227 solver_scalar_type>::do_get(B, bvals_(),
229 ROOTED, this->rowIndexBase_);
233 MultiVecAdapter<Vector>,
234 solver_scalar_type>::do_get(B, bvals_(),
241#ifdef HAVE_AMESOS2_TIMERS
242 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
245 const int_t phase = 33;
247 function_map::pardiso( pt_,
248 const_cast<int_t*
>(&maxfct_),
249 const_cast<int_t*
>(&mnum_),
250 const_cast<int_t*
>(&mtype_),
251 const_cast<int_t*
>(&phase),
252 const_cast<int_t*
>(&n_),
253 const_cast<solver_scalar_type*
>(nzvals_view_.data()),
254 const_cast<int_t*
>(rowptr_view_.data()),
255 const_cast<int_t*
>(colind_view_.data()),
256 const_cast<int_t*
>(perm_.getRawPtr()),
258 const_cast<int_t*
>(iparm_),
259 const_cast<int_t*
>(&msglvl_),
260 as<void*>(bvals_.getRawPtr()),
261 as<void*>(xvals_.getRawPtr()), &error );
264 check_pardiso_mkl_error(Amesos2::SOLVE, error);
268#ifdef HAVE_AMESOS2_TIMERS
269 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
272 if ( is_contiguous_ ==
true ) {
274 MultiVecAdapter<Vector>,
275 solver_scalar_type>::do_put(X, xvals_(),
281 MultiVecAdapter<Vector>,
282 solver_scalar_type>::do_put(X, xvals_(),
306 using Teuchos::getIntegralValue;
307 using Teuchos::ParameterEntryValidator;
309 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
311 if( parameterList->isParameter(
"IPARM(2)") )
313 RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry(
"IPARM(2)").validator();
314 parameterList->getEntry(
"IPARM(2)").setValidator(fillin_validator);
315 iparm_[1] = getIntegralValue<int>(*parameterList,
"IPARM(2)");
318 if( parameterList->isParameter(
"IPARM(4)") )
320 RCP<const ParameterEntryValidator> prec_validator = valid_params->getEntry(
"IPARM(4)").validator();
321 parameterList->getEntry(
"IPARM(4)").setValidator(prec_validator);
322 iparm_[3] = getIntegralValue<int>(*parameterList,
"IPARM(4)");
325 if( parameterList->isParameter(
"IPARM(8)") )
327 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IPARM(8)").validator();
328 parameterList->getEntry(
"IPARM(8)").setValidator(refine_validator);
329 iparm_[7] = getIntegralValue<int>(*parameterList,
"IPARM(8)");
332 if( parameterList->isParameter(
"IPARM(10)") )
334 RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry(
"IPARM(10)").validator();
335 parameterList->getEntry(
"IPARM(10)").setValidator(pivot_perturb_validator);
336 iparm_[9] = getIntegralValue<int>(*parameterList,
"IPARM(10)");
341 iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
343 if( parameterList->isParameter(
"IPARM(12)") )
345 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"IPARM(12)").validator();
346 parameterList->getEntry(
"IPARM(12)").setValidator(trans_validator);
347 iparm_[11] = getIntegralValue<int>(*parameterList,
"IPARM(12)");
350 if( parameterList->isParameter(
"IPARM(18)") )
352 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry(
"IPARM(18)").validator();
353 parameterList->getEntry(
"IPARM(18)").setValidator(report_validator);
354 iparm_[17] = getIntegralValue<int>(*parameterList,
"IPARM(18)");
357 if( parameterList->isParameter(
"IPARM(24)") )
359 RCP<const ParameterEntryValidator> par_fact_validator = valid_params->getEntry(
"IPARM(24)").validator();
360 parameterList->getEntry(
"IPARM(24)").setValidator(par_fact_validator);
361 iparm_[23] = getIntegralValue<int>(*parameterList,
"IPARM(24)");
364 if( parameterList->isParameter(
"IPARM(25)") )
366 RCP<const ParameterEntryValidator> par_fbsolve_validator = valid_params->getEntry(
"IPARM(25)").validator();
367 parameterList->getEntry(
"IPARM(25)").setValidator(par_fbsolve_validator);
368 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");
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(24)", toString(iparm_temp[23]),
455 "Parallel factorization control",
456 tuple<string>(
"0",
"1"),
457 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
458 "PARDISO uses the new two-level factorization algorithm"),
462 setStringToIntegralParameter<int>(
"IPARM(25)", toString(iparm_temp[24]),
463 "Parallel forward/backward solve control",
464 tuple<string>(
"0",
"1"),
465 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
466 "PARDISO uses the sequential forward and backward solve"),
470 setStringToIntegralParameter<int>(
"IPARM(60)", toString(iparm_temp[59]),
471 "PARDISO mode (OOC mode)",
472 tuple<string>(
"0",
"2"),
473 tuple<string>(
"In-core PARDISO",
474 "Out-of-core PARDISO. The OOC PARDISO can solve very "
475 "large problems by holding the matrix factors in files "
476 "on the disk. Hence the amount of RAM required by OOC "
477 "PARDISO is significantly reduced."),
481 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
482 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
484 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
485 accept_int.allowInt(
true );
487 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
488 anyNumberParameterEntryValidator(preferred_int, accept_int));
490 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
491 anyNumberParameterEntryValidator(preferred_int, accept_int));
493 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
494 anyNumberParameterEntryValidator(preferred_int, accept_int));
496 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
510#ifdef HAVE_AMESOS2_TIMERS
511 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
515 if( current_phase == PREORDERING )
return(
false );
518 Kokkos::resize(nzvals_view_, this->globalNumNonZeros_);
519 Kokkos::resize(colind_view_, this->globalNumNonZeros_);
520 Kokkos::resize(rowptr_view_, this->globalNumRows_ + 1);
525#ifdef HAVE_AMESOS2_TIMERS
526 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
529 if ( is_contiguous_ ==
true ) {
531 MatrixAdapter<Matrix>,
532 host_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
533 this->matrixA_.ptr(),
534 nzvals_view_, colind_view_, rowptr_view_,
539 MatrixAdapter<Matrix>,
540 host_value_type_array, host_ordinal_type_array, host_size_type_array>::do_get(
541 this->matrixA_.ptr(),
542 nzvals_view_, colind_view_, rowptr_view_,
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