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( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
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_ );
224 if ( is_contiguous_ ==
true ) {
227 solver_scalar_type>::do_get(B, bvals_(),
229 ROOTED, this->rowIndexBase_);
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 ) {
275 solver_scalar_type>::do_put(X, xvals_(),
282 solver_scalar_type>::do_put(X, xvals_(),
292 template <
class Matrix,
class Vector>
297 return( this->globalNumRows_ == this->globalNumCols_ );
301 template <
class Matrix,
class Vector>
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");
404template <
class Matrix,
class Vector>
405Teuchos::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(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");
506template <
class Matrix,
class Vector>
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 ) {
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_,
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_,
551template <
class Matrix,
class Vector>
557 Teuchos::broadcast(*(this->getComm()), 0, &error_i);
559 if( error == 0 )
return;
561 std::string errmsg =
"Other error";
564 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
567 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
570 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
574 "PardisoMKL reported error: 'Zero pivot, numerical "
575 "factorization or iterative refinement problem'";
578 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
581 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
584 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
587 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
590 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
593 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
596 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
600 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
604template <
class Matrix,
class Vector>
617 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
618 std::invalid_argument,
619 "Cannot set a real Pardiso matrix type with scalar type complex" );
622 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
623 std::invalid_argument,
624 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
627 TEUCHOS_TEST_FOR_EXCEPTION(
true,
628 std::invalid_argument,
629 "Symmetric matrices are not yet supported by the Amesos2 interface" );
635template <
class Matrix,
class Vector>
638template <
class Matrix,
class Vector>
639const typename PardisoMKL<Matrix,Vector>::int_t
642template <
class Matrix,
class Vector>
643const typename PardisoMKL<Matrix,Vector>::int_t
646template <
class Matrix,
class Vector>
647const typename PardisoMKL<Matrix,Vector>::int_t
A template class that does nothing useful besides show developers what, in general,...
@ ROOTED
Definition: Amesos2_TypeDecl.hpp:127
@ CONTIGUOUS_AND_ROOTED
Definition: Amesos2_TypeDecl.hpp:128
@ SORTED_INDICES
Definition: Amesos2_TypeDecl.hpp:142
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Amesos2 interface to the PardisoMKL package.
Definition: Amesos2_PardisoMKL_decl.hpp:84
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_PardisoMKL_def.hpp:406
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_PardisoMKL_def.hpp:303
~PardisoMKL()
Destructor.
Definition: Amesos2_PardisoMKL_def.hpp:110
int_t mtype_
The matrix type. We deal only with unsymmetrix matrices.
Definition: Amesos2_PardisoMKL_decl.hpp:291
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
int_t iparm_[64]
Definition: Amesos2_PardisoMKL_decl.hpp:301
int numericFactorization_impl()
PardisoMKL specific numeric factorization.
Definition: Amesos2_PardisoMKL_def.hpp:176
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_PardisoMKL_def.hpp:294
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_PardisoMKL_def.hpp:508
void set_pardiso_mkl_matrix_type(int_t mtype=0)
Definition: Amesos2_PardisoMKL_def.hpp:606
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
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_PardisoMKL_def.hpp:133
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:553
void * pt_[64]
PardisoMKL internal data address pointer.
Definition: Amesos2_PardisoMKL_decl.hpp:289
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using PardisoMKL.
Definition: Amesos2_PardisoMKL_def.hpp:144
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition: Amesos2_SolverCore_decl.hpp:106
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:267
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:663
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:373