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_)
84 , is_contiguous_(true)
89 PMKL::_INTEGER_t iparm_temp[64];
90 PMKL::_INTEGER_t mtype_temp =
mtype_;
91 PMKL::pardisoinit(
pt_, &mtype_temp, iparm_temp);
93 for(
int i = 0; i < 64; ++i ){
98 if( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
106 #ifdef HAVE_AMESOS2_DEBUG
112 template <
class Matrix,
class Vector>
119 void *bdummy, *xdummy;
123 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
124 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
125 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
126 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
127 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
130 check_pardiso_mkl_error(Amesos2::CLEAN, error);
134 template<
class Matrix,
class Vector>
145 template <
class Matrix,
class Vector>
152 #ifdef HAVE_AMESOS2_TIMERS
153 Teuchos::TimeMonitor symbFactTimer( this->timers_.symFactTime_ );
157 void *bdummy, *xdummy;
159 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
160 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
161 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
162 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
163 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
166 check_pardiso_mkl_error(Amesos2::SYMBFACT, error);
171 this->setNnzLU(iparm_[17]);
177 template <
class Matrix,
class Vector>
184 #ifdef HAVE_AMESOS2_TIMERS
185 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
189 void *bdummy, *xdummy;
191 function_map::pardiso( pt_, const_cast<int_t*>(&maxfct_),
192 const_cast<int_t*>(&mnum_), &mtype_, &phase, &n_,
193 nzvals_.getRawPtr(), rowptr_.getRawPtr(),
194 colind_.getRawPtr(), perm_.getRawPtr(), &nrhs_, iparm_,
195 const_cast<int_t*
>(&msglvl_), &bdummy, &xdummy, &error );
198 check_pardiso_mkl_error(Amesos2::NUMFACT, error);
204 template <
class Matrix,
class Vector>
214 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
215 nrhs_ = as<int_t>(X->getGlobalNumVectors());
217 const size_t val_store_size = as<size_t>(ld_rhs * nrhs_);
218 xvals_.resize(val_store_size);
219 bvals_.resize(val_store_size);
222 #ifdef HAVE_AMESOS2_TIMERS
223 Teuchos::TimeMonitor mvConvTimer( this->timers_.vecConvTime_ );
224 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
227 if ( is_contiguous_ ==
true ) {
230 solver_scalar_type>::do_get(B, bvals_(),
232 ROOTED, this->rowIndexBase_);
237 solver_scalar_type>::do_get(B, bvals_(),
239 CONTIGUOUS_AND_ROOTED, this->rowIndexBase_);
244 #ifdef HAVE_AMESOS2_TIMERS
245 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
248 const int_t phase = 33;
250 function_map::pardiso( pt_,
251 const_cast<int_t*>(&maxfct_),
252 const_cast<int_t*>(&mnum_),
253 const_cast<int_t*>(&mtype_),
254 const_cast<int_t*>(&phase),
255 const_cast<int_t*>(&n_),
256 const_cast<solver_scalar_type*>(nzvals_.getRawPtr()),
257 const_cast<int_t*>(rowptr_.getRawPtr()),
258 const_cast<int_t*>(colind_.getRawPtr()),
259 const_cast<int_t*>(perm_.getRawPtr()),
261 const_cast<int_t*>(iparm_),
262 const_cast<int_t*
>(&msglvl_),
263 as<void*>(bvals_.getRawPtr()),
264 as<void*>(xvals_.getRawPtr()), &error );
267 check_pardiso_mkl_error(Amesos2::SOLVE, error);
271 #ifdef HAVE_AMESOS2_TIMERS
272 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
275 if ( is_contiguous_ ==
true ) {
278 solver_scalar_type>::do_put(X, xvals_(),
285 solver_scalar_type>::do_put(X, xvals_(),
287 CONTIGUOUS_AND_ROOTED);
295 template <
class Matrix,
class Vector>
300 return( this->globalNumRows_ == this->globalNumCols_ );
304 template <
class Matrix,
class Vector>
309 using Teuchos::getIntegralValue;
310 using Teuchos::ParameterEntryValidator;
312 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
314 if( parameterList->isParameter(
"IPARM(2)") )
316 RCP<const ParameterEntryValidator> fillin_validator = valid_params->getEntry(
"IPARM(2)").validator();
317 parameterList->getEntry(
"IPARM(2)").setValidator(fillin_validator);
318 iparm_[1] = getIntegralValue<int>(*parameterList,
"IPARM(2)");
321 if( parameterList->isParameter(
"IPARM(4)") )
323 RCP<const ParameterEntryValidator> prec_validator = valid_params->getEntry(
"IPARM(4)").validator();
324 parameterList->getEntry(
"IPARM(4)").setValidator(prec_validator);
325 iparm_[3] = getIntegralValue<int>(*parameterList,
"IPARM(4)");
328 if( parameterList->isParameter(
"IPARM(8)") )
330 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IPARM(8)").validator();
331 parameterList->getEntry(
"IPARM(8)").setValidator(refine_validator);
332 iparm_[7] = getIntegralValue<int>(*parameterList,
"IPARM(8)");
335 if( parameterList->isParameter(
"IPARM(10)") )
337 RCP<const ParameterEntryValidator> pivot_perturb_validator = valid_params->getEntry(
"IPARM(10)").validator();
338 parameterList->getEntry(
"IPARM(10)").setValidator(pivot_perturb_validator);
339 iparm_[9] = getIntegralValue<int>(*parameterList,
"IPARM(10)");
344 iparm_[11] = this->control_.useTranspose_ ? 2 : 0;
346 if( parameterList->isParameter(
"IPARM(12)") )
348 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"IPARM(12)").validator();
349 parameterList->getEntry(
"IPARM(12)").setValidator(trans_validator);
350 iparm_[11] = getIntegralValue<int>(*parameterList,
"IPARM(12)");
353 if( parameterList->isParameter(
"IPARM(18)") )
355 RCP<const ParameterEntryValidator> report_validator = valid_params->getEntry(
"IPARM(18)").validator();
356 parameterList->getEntry(
"IPARM(18)").setValidator(report_validator);
357 iparm_[17] = getIntegralValue<int>(*parameterList,
"IPARM(18)");
360 if( parameterList->isParameter(
"IPARM(24)") )
362 RCP<const ParameterEntryValidator> par_fact_validator = valid_params->getEntry(
"IPARM(24)").validator();
363 parameterList->getEntry(
"IPARM(24)").setValidator(par_fact_validator);
364 iparm_[23] = getIntegralValue<int>(*parameterList,
"IPARM(24)");
367 if( parameterList->isParameter(
"IPARM(25)") )
369 RCP<const ParameterEntryValidator> par_fbsolve_validator = valid_params->getEntry(
"IPARM(25)").validator();
370 parameterList->getEntry(
"IPARM(25)").setValidator(par_fbsolve_validator);
371 iparm_[24] = getIntegralValue<int>(*parameterList,
"IPARM(25)");
374 if( parameterList->isParameter(
"IPARM(60)") )
376 RCP<const ParameterEntryValidator> ooc_validator = valid_params->getEntry(
"IPARM(60)").validator();
377 parameterList->getEntry(
"IPARM(60)").setValidator(ooc_validator);
378 iparm_[59] = getIntegralValue<int>(*parameterList,
"IPARM(60)");
381 if( parameterList->isParameter(
"IsContiguous") ){
382 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
407 template <
class Matrix,
class Vector>
408 Teuchos::RCP<const Teuchos::ParameterList>
414 using Teuchos::tuple;
415 using Teuchos::toString;
416 using Teuchos::EnhancedNumberValidator;
417 using Teuchos::setStringToIntegralParameter;
418 using Teuchos::anyNumberParameterEntryValidator;
420 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
422 if( is_null(valid_params) ){
423 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
427 PMKL::_INTEGER_t mtype_temp = mtype_;
428 PMKL::_INTEGER_t iparm_temp[64];
429 PMKL::pardisoinit(pt_dummy,
430 const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
431 const_cast<PMKL::_INTEGER_t*>(iparm_temp));
433 setStringToIntegralParameter<int>(
"IPARM(2)", toString(iparm_temp[1]),
434 "Fill-in reducing ordering for the input matrix",
435 tuple<string>(
"0",
"2",
"3"),
436 tuple<string>(
"The minimum degree algorithm",
437 "Nested dissection algorithm from METIS",
438 "OpenMP parallel nested dissection algorithm"),
442 Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
443 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
444 iparm_4_validator->setMin(0);
445 pl->set(
"IPARM(4)" , as<int>(iparm_temp[3]) ,
"Preconditioned CGS/CG",
448 setStringToIntegralParameter<int>(
"IPARM(12)", toString(iparm_temp[11]),
449 "Solve with transposed or conjugate transposed matrix A",
450 tuple<string>(
"0",
"1",
"2"),
451 tuple<string>(
"Non-transposed",
452 "Conjugate-transposed",
457 setStringToIntegralParameter<int>(
"IPARM(24)", toString(iparm_temp[23]),
458 "Parallel factorization control",
459 tuple<string>(
"0",
"1"),
460 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
461 "PARDISO uses the new two-level factorization algorithm"),
465 setStringToIntegralParameter<int>(
"IPARM(25)", toString(iparm_temp[24]),
466 "Parallel forward/backward solve control",
467 tuple<string>(
"0",
"1"),
468 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
469 "PARDISO uses the sequential forward and backward solve"),
473 setStringToIntegralParameter<int>(
"IPARM(60)", toString(iparm_temp[59]),
474 "PARDISO mode (OOC mode)",
475 tuple<string>(
"0",
"2"),
476 tuple<string>(
"In-core PARDISO",
477 "Out-of-core PARDISO. The OOC PARDISO can solve very "
478 "large problems by holding the matrix factors in files "
479 "on the disk. Hence the amount of RAM required by OOC "
480 "PARDISO is significantly reduced."),
484 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
485 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
487 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
488 accept_int.allowInt(
true );
490 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
491 anyNumberParameterEntryValidator(preferred_int, accept_int));
493 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
494 anyNumberParameterEntryValidator(preferred_int, accept_int));
496 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
497 anyNumberParameterEntryValidator(preferred_int, accept_int));
499 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
509 template <
class Matrix,
class Vector>
513 #ifdef HAVE_AMESOS2_TIMERS
514 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
518 if( current_phase == PREORDERING )
return(
false );
521 nzvals_.resize(this->globalNumNonZeros_);
522 colind_.resize(this->globalNumNonZeros_);
523 rowptr_.resize(this->globalNumRows_ + 1);
528 #ifdef HAVE_AMESOS2_TIMERS
529 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
532 if ( is_contiguous_ ==
true ) {
536 int_t,int_t>::do_get(this->matrixA_.ptr(),
537 nzvals_(), colind_(), rowptr_(),
538 nnz_ret, ROOTED, SORTED_INDICES, this->rowIndexBase_);
544 int_t,int_t>::do_get(this->matrixA_.ptr(),
545 nzvals_(), colind_(), rowptr_(),
546 nnz_ret, CONTIGUOUS_AND_ROOTED, SORTED_INDICES, this->rowIndexBase_);
554 template <
class Matrix,
class Vector>
560 Teuchos::broadcast(*(this->getComm()), 0, &error_i);
562 if( error == 0 )
return;
564 std::string errmsg =
"Other error";
567 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
570 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
573 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
577 "PardisoMKL reported error: 'Zero pivot, numerical "
578 "factorization or iterative refinement problem'";
581 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
584 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
587 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
590 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
593 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
596 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
599 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
603 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
607 template <
class Matrix,
class Vector>
620 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
621 std::invalid_argument,
622 "Cannot set a real Pardiso matrix type with scalar type complex" );
625 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
626 std::invalid_argument,
627 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
630 TEUCHOS_TEST_FOR_EXCEPTION(
true,
631 std::invalid_argument,
632 "Symmetric matrices are not yet supported by the Amesos2 interface" );
638 template <
class Matrix,
class Vector>
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
649 template <
class Matrix,
class Vector>
650 const typename PardisoMKL<Matrix,Vector>::int_t
656 #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:113
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:665
void set_pardiso_mkl_matrix_type(int_t mtype=0)
Definition: Amesos2_PardisoMKL_def.hpp:609
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:266
void * pt_[64]
PardisoMKL internal data address pointer.
Definition: Amesos2_PardisoMKL_decl.hpp:285
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:179
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using PardisoMKL.
Definition: Amesos2_PardisoMKL_def.hpp:147
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:511
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_PardisoMKL_def.hpp:306
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:206
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:556
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:322
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:409
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_PardisoMKL_def.hpp:297
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:136