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_)
88 PMKL::_INTEGER_t iparm_temp[64];
89 PMKL::_INTEGER_t mtype_temp =
mtype_;
90 PMKL::pardisoinit(
pt_, &mtype_temp, iparm_temp);
92 for(
int i = 0; i < 64; ++i ){
97 if( Meta::is_same<solver_magnitude_type, PMKL::_REAL_t>::value ){
105 #ifdef HAVE_AMESOS2_DEBUG 111 template <
class Matrix,
class Vector>
118 void *bdummy, *xdummy;
122 function_map::pardiso(
pt_, const_cast<int_t*>(&maxfct_),
123 const_cast<int_t*>(&mnum_), &
mtype_, &phase, &
n_,
126 const_cast<int_t*
>(&
msglvl_), &bdummy, &xdummy, &error );
133 template<
class Matrix,
class Vector>
144 template <
class Matrix,
class Vector>
151 #ifdef HAVE_AMESOS2_TIMERS 152 Teuchos::TimeMonitor symbFactTimer( this->
timers_.symFactTime_ );
156 void *bdummy, *xdummy;
158 function_map::pardiso(
pt_, const_cast<int_t*>(&maxfct_),
159 const_cast<int_t*>(&mnum_), &
mtype_, &phase, &
n_,
162 const_cast<int_t*
>(&
msglvl_), &bdummy, &xdummy, &error );
176 template <
class Matrix,
class Vector>
183 #ifdef HAVE_AMESOS2_TIMERS 184 Teuchos::TimeMonitor numFactTimer( this->
timers_.numFactTime_ );
188 void *bdummy, *xdummy;
190 function_map::pardiso(
pt_, const_cast<int_t*>(&maxfct_),
191 const_cast<int_t*>(&mnum_), &
mtype_, &phase, &
n_,
194 const_cast<int_t*
>(&
msglvl_), &bdummy, &xdummy, &error );
203 template <
class Matrix,
class Vector>
213 const global_size_type ld_rhs = this->
root_ ? X->getGlobalLength() : 0;
214 nrhs_ = as<int_t>(X->getGlobalNumVectors());
216 const size_t val_store_size = as<size_t>(ld_rhs *
nrhs_);
217 xvals_.resize(val_store_size);
218 bvals_.resize(val_store_size);
221 #ifdef HAVE_AMESOS2_TIMERS 222 Teuchos::TimeMonitor mvConvTimer( this->
timers_.vecConvTime_ );
223 Teuchos::TimeMonitor redistTimer( this->
timers_.vecRedistTime_ );
227 solver_scalar_type>::do_get(B,
bvals_(),
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_.getRawPtr()),
246 const_cast<int_t*>(
rowptr_.getRawPtr()),
247 const_cast<int_t*>(
colind_.getRawPtr()),
248 const_cast<int_t*>(
perm_.getRawPtr()),
250 const_cast<int_t*>(
iparm_),
252 as<void*>(
bvals_.getRawPtr()),
253 as<void*>(
xvals_.getRawPtr()), &error );
260 #ifdef HAVE_AMESOS2_TIMERS 261 Teuchos::TimeMonitor redistTimer(this->
timers_.vecRedistTime_);
266 solver_scalar_type>::do_put(X,
xvals_(),
275 template <
class Matrix,
class Vector>
284 template <
class Matrix,
class Vector>
288 iparm_[1] = validators[2]->getIntegralValue(*parameterList,
"IPARM(2)",
289 validators[2]->getDefaultParameterName());
290 iparm_[3] = parameterList->get<
int>(
"IPARM(4)" , (int)
iparm_[3]);
291 iparm_[7] = parameterList->get<
int>(
"IPARM(8)" , (int)
iparm_[7]);
292 iparm_[9] = parameterList->get<
int>(
"IPARM(10)", (int)
iparm_[9]);
293 iparm_[17] = parameterList->get<
int>(
"IPARM(18)", (int)
iparm_[17]);
294 iparm_[23] = validators[24]->getIntegralValue(*parameterList,
"IPARM(24)",
295 validators[24]->getDefaultParameterName());
296 iparm_[24] = validators[25]->getIntegralValue(*parameterList,
"IPARM(25)",
297 validators[25]->getDefaultParameterName());
298 iparm_[59] = validators[60]->getIntegralValue(*parameterList,
"IPARM(60)",
299 validators[60]->getDefaultParameterName());
323 template <
class Matrix,
class Vector>
324 Teuchos::RCP<const Teuchos::ParameterList>
330 using Teuchos::tuple;
331 using Teuchos::toString;
332 using Teuchos::EnhancedNumberValidator;
333 using Teuchos::setStringToIntegralParameter;
334 using Teuchos::anyNumberParameterEntryValidator;
335 using Teuchos::stringToIntegralParameterEntryValidator;
336 typedef Teuchos::StringToIntegralParameterEntryValidator<int> STIPEV;
337 Teuchos::AnyNumberParameterEntryValidator::EPreferredType preferred_int =
338 Teuchos::AnyNumberParameterEntryValidator::PREFER_INT;
340 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
342 if( is_null(valid_params) ){
343 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
347 PMKL::_INTEGER_t mtype_temp =
mtype_;
348 PMKL::_INTEGER_t iparm_temp[64];
349 PMKL::pardisoinit(pt_dummy,
350 const_cast<PMKL::_INTEGER_t*>(&mtype_temp),
351 const_cast<PMKL::_INTEGER_t*>(iparm_temp));
354 RCP<STIPEV> iparm_2_validator
355 = stringToIntegralParameterEntryValidator<int>(tuple<string>(
"0",
"2",
"3"),
356 tuple<string>(
"The minimum degree algorithm",
357 "Nested dissection algorithm from METIS",
358 "OpenMP parallel nested dissection algorithm"),
360 toString(iparm_temp[1]));
361 validators.insert( std::pair<
int,RCP<STIPEV> >(2, iparm_2_validator) );
363 Teuchos::RCP<EnhancedNumberValidator<int> > iparm_4_validator
364 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
365 iparm_4_validator->setMin(0);
367 RCP<STIPEV> iparm_24_validator
368 = stringToIntegralParameterEntryValidator<int>(tuple<string>(
"0",
"1"),
369 tuple<string>(
"PARDISO uses the previous algorithm for factorization",
370 "PARDISO uses the new two-level factorization algorithm"),
372 toString(iparm_temp[23]));
373 validators.insert( std::pair<
int,RCP<STIPEV> >(24, iparm_24_validator) );
375 RCP<STIPEV> iparm_25_validator
376 = stringToIntegralParameterEntryValidator<int>(tuple<string>(
"0",
"1"),
377 tuple<string>(
"PARDISO uses the parallel algorithm for the solve step",
378 "PARDISO uses the sequential forward and backward solve"),
380 toString(iparm_temp[24]));
381 validators.insert( std::pair<
int,RCP<STIPEV> >(25, iparm_25_validator) );
383 RCP<STIPEV> iparm_60_validator
384 = stringToIntegralParameterEntryValidator<int>(tuple<string>(
"0",
"2"),
385 tuple<string>(
"In-core PARDISO",
386 "Out-of-core PARDISO. The OOC PARDISO can solve very " 387 "large problems by holding the matrix factors in files " 388 "on the disk. Hence the amount of RAM required by OOC " 389 "PARDISO is significantly reduced."),
391 toString(iparm_temp[59]));
392 validators.insert( std::pair<
int,RCP<STIPEV> >(60, iparm_60_validator) );
394 Teuchos::AnyNumberParameterEntryValidator::AcceptedTypes accept_int(
false );
395 accept_int.allowInt(
true );
397 pl->set(
"IPARM(2)" , validators[2]->getDefaultParameterName(),
398 "Fill-in reducing ordering for the input matrix", validators[2]);
400 pl->set(
"IPARM(4)" , as<int>(iparm_temp[3]) ,
"Preconditioned CGS/CG",
403 pl->set(
"IPARM(8)" , as<int>(iparm_temp[8]) ,
"Iterative refinement step",
404 anyNumberParameterEntryValidator(preferred_int, accept_int));
406 pl->set(
"IPARM(10)", as<int>(iparm_temp[9]) ,
"Pivoting perturbation",
407 anyNumberParameterEntryValidator(preferred_int, accept_int));
409 pl->set(
"IPARM(18)", as<int>(iparm_temp[17]),
"Report the number of non-zero elements in the factors",
410 anyNumberParameterEntryValidator(preferred_int, accept_int));
412 pl->set(
"IPARM(24)", validators[24]->getDefaultParameterName(),
413 "Parallel factorization control", validators[24]);
415 pl->set(
"IPARM(25)", validators[25]->getDefaultParameterName(),
416 "Parallel forward/backward solve control", validators[25]);
418 pl->set(
"IPARM(60)", validators[60]->getDefaultParameterName(),
419 "PARDISO mode (OOC mode)", validators[60]);
429 template <
class Matrix,
class Vector>
433 #ifdef HAVE_AMESOS2_TIMERS 434 Teuchos::TimeMonitor convTimer(this->
timers_.mtxConvTime_);
438 if( current_phase == PREORDERING )
return(
false );
448 #ifdef HAVE_AMESOS2_TIMERS 449 Teuchos::TimeMonitor mtxRedistTimer( this->
timers_.mtxRedistTime_ );
455 int_t,int_t>::do_get(this->
matrixA_.ptr(),
464 template <
class Matrix,
class Vector>
470 Teuchos::broadcast(*(this->
getComm()), 0, &error_i);
472 if( error == 0 )
return;
474 std::string errmsg =
"Other error";
477 errmsg =
"PardisoMKL reported error: 'Input inconsistent'";
480 errmsg =
"PardisoMKL reported error: 'Not enough memory'";
483 errmsg =
"PardisoMKL reported error: 'Reordering problem'";
487 "PardisoMKL reported error: 'Zero pivot, numerical " 488 "factorization or iterative refinement problem'";
491 errmsg =
"PardisoMKL reported error: 'Unclassified (internal) error'";
494 errmsg =
"PardisoMKL reported error: 'Reordering failed'";
497 errmsg =
"PardisoMKL reported error: 'Diagonal matrix is singular'";
500 errmsg =
"PardisoMKL reported error: '32-bit integer overflow problem'";
503 errmsg =
"PardisoMKL reported error: 'Not enough memory for OOC'";
506 errmsg =
"PardisoMKL reported error: 'Problems with opening OOC temporary files'";
509 errmsg =
"PardisoMKL reported error: 'Read/write problem with OOC data file'";
513 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, errmsg );
517 template <
class Matrix,
class Vector>
530 TEUCHOS_TEST_FOR_EXCEPTION( complex_,
531 std::invalid_argument,
532 "Cannot set a real Pardiso matrix type with scalar type complex" );
535 TEUCHOS_TEST_FOR_EXCEPTION( !complex_,
536 std::invalid_argument,
537 "Cannot set a complex Pardiso matrix type with non-complex scalars" );
540 TEUCHOS_TEST_FOR_EXCEPTION(
true,
541 std::invalid_argument,
542 "Symmetric matrices are not yet supported by the Amesos2 interface" );
548 template <
class Matrix,
class Vector>
551 template <
class Matrix,
class Vector>
552 const typename PardisoMKL<Matrix,Vector>::int_t
555 template <
class Matrix,
class Vector>
556 const typename PardisoMKL<Matrix,Vector>::int_t
559 template <
class Matrix,
class Vector>
560 const typename PardisoMKL<Matrix,Vector>::int_t
566 #endif // AMESOS2_PARDISOMKL_DEF_HPP Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Definition: Amesos2_TypeDecl.hpp:141
~PardisoMKL()
Destructor.
Definition: Amesos2_PardisoMKL_def.hpp:112
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
Teuchos::Array< solver_scalar_type > nzvals_
Stores the values of the nonzero entries for PardisoMKL.
Definition: Amesos2_PardisoMKL_decl.hpp:273
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:591
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
int_t n_
Number of equations in the sparse linear system.
Definition: Amesos2_PardisoMKL_decl.hpp:288
bool root_
If true, then this is the root processor.
Definition: Amesos2_SolverCore_decl.hpp:507
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
int_t nrhs_
number of righthand-side vectors
Definition: Amesos2_PardisoMKL_decl.hpp:292
void set_pardiso_mkl_matrix_type(int_t mtype=0)
Definition: Amesos2_PardisoMKL_def.hpp:519
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
static const int_t msglvl_
The messaging level. Set to 1 if you wish for Pardiso MKL to print statistical info.
Definition: Amesos2_PardisoMKL_decl.hpp:299
Teuchos::Array< solver_scalar_type > bvals_
Persisting, contiguous, 1D store for B.
Definition: Amesos2_PardisoMKL_decl.hpp:281
void * pt_[64]
PardisoMKL internal data address pointer.
Definition: Amesos2_PardisoMKL_decl.hpp:284
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Teuchos::Array< int_t > perm_
Permutation vector.
Definition: Amesos2_PardisoMKL_decl.hpp:290
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:178
Teuchos::Array< int_t > rowptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_PardisoMKL_decl.hpp:277
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition: Amesos2_SolverCore_decl.hpp:363
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using PardisoMKL.
Definition: Amesos2_PardisoMKL_def.hpp:146
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:431
Definition: Amesos2_Cholmod_TypeMap.hpp:92
Teuchos::Array< solver_scalar_type > xvals_
Persisting, contiguous, 1D store for X.
Definition: Amesos2_PardisoMKL_decl.hpp:279
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_PardisoMKL_def.hpp:286
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
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:205
void setNnzLU(size_t nnz)
Set the number of non-zero values in the and factors.
Definition: Amesos2_SolverCore_decl.hpp:452
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:466
Teuchos::Array< int_t > colind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_PardisoMKL_decl.hpp:275
Definition: Amesos2_TypeDecl.hpp:127
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:296
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:175
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_PardisoMKL_def.hpp:325
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:286
int_t iparm_[64]
Definition: Amesos2_PardisoMKL_decl.hpp:296
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_PardisoMKL_def.hpp:135