53 #ifndef AMESOS2_SUPERLU_DEF_HPP 54 #define AMESOS2_SUPERLU_DEF_HPP 56 #include <Teuchos_Tuple.hpp> 57 #include <Teuchos_ParameterList.hpp> 58 #include <Teuchos_StandardParameterEntryValidators.hpp> 66 template <
class Matrix,
class Vector>
68 Teuchos::RCP<const Matrix> A,
69 Teuchos::RCP<Vector> X,
70 Teuchos::RCP<const Vector> B )
81 SLU::set_default_options(&(data_.options));
83 data_.options.PrintStat = SLU::NO;
85 SLU::StatInit(&(data_.stat));
93 data_.relax = SLU::sp_ienv(2);
94 data_.panel_size = SLU::sp_ienv(1);
100 data_.X.Store = NULL;
101 data_.B.Store = NULL;
107 template <
class Matrix,
class Vector>
115 SLU::StatFree( &(data_.stat) ) ;
118 if ( data_.A.Store != NULL ){
119 SLU::Destroy_SuperMatrix_Store( &(data_.A) );
123 if ( data_.L.Store != NULL ){
124 SLU::Destroy_SuperNode_Matrix( &(data_.L) );
125 SLU::Destroy_CompCol_Matrix( &(data_.U) );
129 template <
class Matrix,
class Vector>
133 std::ostringstream oss;
134 oss <<
"SuperLU solver interface";
136 oss <<
", \"ILUTP\" : {";
137 oss <<
"drop tol = " << data_.options.ILU_DropTol;
138 oss <<
", fill factor = " << data_.options.ILU_FillFactor;
139 oss <<
", fill tol = " << data_.options.ILU_FillTol;
140 switch(data_.options.ILU_MILU) {
152 oss <<
", regular ILU";
154 switch(data_.options.ILU_Norm) {
163 oss <<
", infinity-norm";
167 oss <<
", direct solve";
183 template<
class Matrix,
class Vector>
195 int permc_spec = data_.options.ColPerm;
196 if ( permc_spec != SLU::MY_PERMC && this->root_ ){
197 #ifdef HAVE_AMESOS2_TIMERS 198 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
201 SLU::get_perm_c(permc_spec, &(data_.A), data_.perm_c.getRawPtr());
208 template <
class Matrix,
class Vector>
223 same_symbolic_ =
false;
224 data_.options.Fact = SLU::DOFACT;
230 template <
class Matrix,
class Vector>
239 if ( !same_symbolic_ && data_.L.Store != NULL ){
240 SLU::Destroy_SuperNode_Matrix( &(data_.L) );
241 SLU::Destroy_CompCol_Matrix( &(data_.U) );
242 data_.L.Store = NULL;
243 data_.U.Store = NULL;
246 if( same_symbolic_ ) data_.options.Fact = SLU::SamePattern_SameRowPerm;
251 #ifdef HAVE_AMESOS2_DEBUG 252 TEUCHOS_TEST_FOR_EXCEPTION( data_.A.ncol != as<int>(this->globalNumCols_),
254 "Error in converting to SuperLU SuperMatrix: wrong number of global columns." );
255 TEUCHOS_TEST_FOR_EXCEPTION( data_.A.nrow != as<int>(this->globalNumRows_),
257 "Error in converting to SuperLU SuperMatrix: wrong number of global rows." );
260 if( data_.options.Equil == SLU::YES ){
261 magnitude_type rowcnd, colcnd, amax;
265 function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
266 data_.C.getRawPtr(), &rowcnd, &colcnd,
268 TEUCHOS_TEST_FOR_EXCEPTION
269 (info2 < 0, std::runtime_error,
270 "SuperLU's gsequ function returned with status " << info2 <<
" < 0." 271 " This means that argument " << (-info2) <<
" given to the function" 272 " had an illegal value.");
274 if (info2 <= data_.A.nrow) {
275 TEUCHOS_TEST_FOR_EXCEPTION
276 (
true, std::runtime_error,
"SuperLU's gsequ function returned with " 277 "info = " << info2 <<
" > 0, and info <= A.nrow = " << data_.A.nrow
278 <<
". This means that row " << info2 <<
" of A is exactly zero.");
280 else if (info2 > data_.A.ncol) {
281 TEUCHOS_TEST_FOR_EXCEPTION
282 (
true, std::runtime_error,
"SuperLU's gsequ function returned with " 283 "info = " << info2 <<
" > 0, and info > A.ncol = " << data_.A.ncol
284 <<
". This means that column " << (info2 - data_.A.nrow) <<
" of " 285 "A is exactly zero.");
288 TEUCHOS_TEST_FOR_EXCEPTION
289 (
true, std::runtime_error,
"SuperLU's gsequ function returned " 290 "with info = " << info2 <<
" > 0, but its value is not in the " 291 "range permitted by the documentation. This should never happen " 292 "(it appears to be SuperLU's fault).");
297 function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
298 data_.C.getRawPtr(), rowcnd, colcnd,
299 amax, &(data_.equed));
308 SLU::sp_preorder(&(data_.options), &(data_.A), data_.perm_c.getRawPtr(),
309 data_.etree.getRawPtr(), &(data_.AC));
312 #ifdef HAVE_AMESOS2_TIMERS 313 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
316 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG 317 std::cout <<
"Superlu:: Before numeric factorization" << std::endl;
318 std::cout <<
"nzvals_ : " << nzvals_.toString() << std::endl;
319 std::cout <<
"rowind_ : " << rowind_.toString() << std::endl;
320 std::cout <<
"colptr_ : " << colptr_.toString() << std::endl;
323 if(ILU_Flag_==
false) {
324 function_map::gstrf(&(data_.options), &(data_.AC),
325 data_.relax, data_.panel_size, data_.etree.getRawPtr(),
326 NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
327 &(data_.L), &(data_.U), &(data_.stat), &info);
330 function_map::gsitrf(&(data_.options), &(data_.AC),
331 data_.relax, data_.panel_size, data_.etree.getRawPtr(),
332 NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
333 &(data_.L), &(data_.U), &(data_.stat), &info);
338 SLU::Destroy_CompCol_Permuted( &(data_.AC) );
341 this->setNnzLU( as<size_t>(((SLU::SCformat*)data_.L.Store)->nnz +
342 ((SLU::NCformat*)data_.U.Store)->nnz) );
346 Teuchos::broadcast(*(this->matrixA_->getComm()), 0, &info);
348 global_size_type info_st = as<global_size_type>(info);
349 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
351 "Factorization complete, but matrix is singular. Division by zero eminent");
352 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
354 "Memory allocation failure in Superlu factorization");
356 data_.options.Fact = SLU::FACTORED;
357 same_symbolic_ =
true;
363 template <
class Matrix,
class Vector>
370 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
371 const size_t nrhs = X->getGlobalNumVectors();
373 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
374 Teuchos::Array<slu_type> xValues(val_store_size);
375 Teuchos::Array<slu_type> bValues(val_store_size);
378 #ifdef HAVE_AMESOS2_TIMERS 379 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
380 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
383 slu_type>::do_get(B, bValues(),
385 ROOTED, this->rowIndexBase_);
390 magnitude_type rpg, rcond;
392 data_.ferr.resize(nrhs);
393 data_.berr.resize(nrhs);
396 #ifdef HAVE_AMESOS2_TIMERS 397 Teuchos::TimeMonitor mvConvTimer(this->timers_.vecConvTime_);
399 SLU::Dtype_t dtype = type_map::dtype;
400 int i_ld_rhs = as<int>(ld_rhs);
401 function_map::create_Dense_Matrix(&(data_.B), i_ld_rhs, as<int>(nrhs),
402 bValues.getRawPtr(), i_ld_rhs,
403 SLU::SLU_DN, dtype, SLU::SLU_GE);
404 function_map::create_Dense_Matrix(&(data_.X), i_ld_rhs, as<int>(nrhs),
405 xValues.getRawPtr(), i_ld_rhs,
406 SLU::SLU_DN, dtype, SLU::SLU_GE);
413 #ifdef HAVE_AMESOS2_TIMERS 414 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
417 if(ILU_Flag_==
false) {
418 function_map::gssvx(&(data_.options), &(data_.A),
419 data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
420 data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
421 data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
422 &(data_.X), &rpg, &rcond, data_.ferr.getRawPtr(),
423 data_.berr.getRawPtr(), &(data_.mem_usage), &(data_.stat), &ierr);
426 function_map::gsisx(&(data_.options), &(data_.A),
427 data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
428 data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
429 data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
430 &(data_.X), &rpg, &rcond, &(data_.mem_usage), &(data_.stat), &ierr);
436 SLU::Destroy_SuperMatrix_Store( &(data_.X) );
437 SLU::Destroy_SuperMatrix_Store( &(data_.B) );
438 data_.X.Store = NULL;
439 data_.B.Store = NULL;
443 Teuchos::broadcast(*(this->getComm()), 0, &ierr);
445 global_size_type ierr_st = as<global_size_type>(ierr);
446 TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
447 std::invalid_argument,
448 "Argument " << -ierr <<
" to SuperLU xgssvx had illegal value" );
449 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st <= this->globalNumCols_,
451 "Factorization complete, but U is exactly singular" );
452 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st > this->globalNumCols_ + 1,
454 "SuperLU allocated " << ierr - this->globalNumCols_ <<
" bytes of " 455 "memory before allocation failure occured." );
459 #ifdef HAVE_AMESOS2_TIMERS 460 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
466 ROOTED, this->rowIndexBase_);
474 template <
class Matrix,
class Vector>
481 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
485 template <
class Matrix,
class Vector>
490 using Teuchos::getIntegralValue;
491 using Teuchos::ParameterEntryValidator;
493 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
495 ILU_Flag_ = parameterList->get<
bool>(
"ILU_Flag",
false);
497 SLU::ilu_set_default_options(&(data_.options));
499 data_.options.PrintStat = SLU::NO;
502 data_.options.Trans = this->control_.useTranspose_ ? SLU::TRANS : SLU::NOTRANS;
504 if( parameterList->isParameter(
"Trans") ){
505 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"Trans").validator();
506 parameterList->getEntry(
"Trans").setValidator(trans_validator);
508 data_.options.Trans = getIntegralValue<SLU::trans_t>(*parameterList,
"Trans");
511 if( parameterList->isParameter(
"IterRefine") ){
512 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IterRefine").validator();
513 parameterList->getEntry(
"IterRefine").setValidator(refine_validator);
515 data_.options.IterRefine = getIntegralValue<SLU::IterRefine_t>(*parameterList,
"IterRefine");
518 if( parameterList->isParameter(
"ColPerm") ){
519 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry(
"ColPerm").validator();
520 parameterList->getEntry(
"ColPerm").setValidator(colperm_validator);
522 data_.options.ColPerm = getIntegralValue<SLU::colperm_t>(*parameterList,
"ColPerm");
525 data_.options.DiagPivotThresh = parameterList->get<
double>(
"DiagPivotThresh", 1.0);
527 bool equil = parameterList->get<
bool>(
"Equil",
true);
528 data_.options.Equil = equil ? SLU::YES : SLU::NO;
530 bool symmetric_mode = parameterList->get<
bool>(
"SymmetricMode",
false);
531 data_.options.SymmetricMode = symmetric_mode ? SLU::YES : SLU::NO;
534 if( parameterList->isParameter(
"RowPerm") ){
535 RCP<const ParameterEntryValidator> rowperm_validator = valid_params->getEntry(
"RowPerm").validator();
536 parameterList->getEntry(
"RowPerm").setValidator(rowperm_validator);
537 data_.options.RowPerm = getIntegralValue<SLU::rowperm_t>(*parameterList,
"RowPerm");
546 data_.options.ILU_DropTol = parameterList->get<
double>(
"ILU_DropTol", 0.0001);
548 data_.options.ILU_FillFactor = parameterList->get<
double>(
"ILU_FillFactor", 10.0);
550 if( parameterList->isParameter(
"ILU_Norm") ) {
551 RCP<const ParameterEntryValidator> norm_validator = valid_params->getEntry(
"ILU_Norm").validator();
552 parameterList->getEntry(
"ILU_Norm").setValidator(norm_validator);
553 data_.options.ILU_Norm = getIntegralValue<SLU::norm_t>(*parameterList,
"ILU_Norm");
556 if( parameterList->isParameter(
"ILU_MILU") ) {
557 RCP<const ParameterEntryValidator> milu_validator = valid_params->getEntry(
"ILU_MILU").validator();
558 parameterList->getEntry(
"ILU_MILU").setValidator(milu_validator);
559 data_.options.ILU_MILU = getIntegralValue<SLU::milu_t>(*parameterList,
"ILU_MILU");
562 data_.options.ILU_FillTol = parameterList->get<
double>(
"ILU_FillTol", 0.01);
567 template <
class Matrix,
class Vector>
568 Teuchos::RCP<const Teuchos::ParameterList>
572 using Teuchos::tuple;
573 using Teuchos::ParameterList;
574 using Teuchos::EnhancedNumberValidator;
575 using Teuchos::setStringToIntegralParameter;
576 using Teuchos::stringToIntegralParameterEntryValidator;
578 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
580 if( is_null(valid_params) ){
581 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
583 setStringToIntegralParameter<SLU::trans_t>(
"Trans",
"NOTRANS",
584 "Solve for the transpose system or not",
585 tuple<string>(
"TRANS",
"NOTRANS",
"CONJ"),
586 tuple<string>(
"Solve with transpose",
587 "Do not solve with transpose",
588 "Solve with the conjugate transpose"),
589 tuple<SLU::trans_t>(SLU::TRANS,
594 setStringToIntegralParameter<SLU::IterRefine_t>(
"IterRefine",
"NOREFINE",
595 "Type of iterative refinement to use",
596 tuple<string>(
"NOREFINE",
"SLU_SINGLE",
"SLU_DOUBLE"),
597 tuple<string>(
"Do not use iterative refinement",
598 "Do single iterative refinement",
599 "Do double iterative refinement"),
600 tuple<SLU::IterRefine_t>(SLU::NOREFINE,
606 setStringToIntegralParameter<SLU::colperm_t>(
"ColPerm",
"COLAMD",
607 "Specifies how to permute the columns of the " 608 "matrix for sparsity preservation",
609 tuple<string>(
"NATURAL",
"MMD_AT_PLUS_A",
610 "MMD_ATA",
"COLAMD"),
611 tuple<string>(
"Natural ordering",
612 "Minimum degree ordering on A^T + A",
613 "Minimum degree ordering on A^T A",
614 "Approximate minimum degree column ordering"),
615 tuple<SLU::colperm_t>(SLU::NATURAL,
621 Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
622 = Teuchos::rcp(
new EnhancedNumberValidator<double>(0.0, 1.0) );
623 pl->set(
"DiagPivotThresh", 1.0,
624 "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
625 diag_pivot_thresh_validator);
627 pl->set(
"Equil",
true,
"Whether to equilibrate the system before solve");
629 pl->set(
"SymmetricMode",
false,
630 "Specifies whether to use the symmetric mode. " 631 "Gives preference to diagonal pivots and uses " 632 "an (A^T + A)-based column permutation.");
636 setStringToIntegralParameter<SLU::rowperm_t>(
"RowPerm",
"LargeDiag",
637 "Type of row permutation strategy to use",
638 tuple<string>(
"NOROWPERM",
"LargeDiag",
"MY_PERMR"),
639 tuple<string>(
"Use natural ordering",
640 "Use weighted bipartite matching algorithm",
641 "Use the ordering given in perm_r input"),
642 tuple<SLU::rowperm_t>(SLU::NOROWPERM,
661 pl->set(
"ILU_DropTol", 0.0001,
"ILUT drop tolerance");
663 pl->set(
"ILU_FillFactor", 10.0,
"ILUT fill factor");
665 setStringToIntegralParameter<SLU::norm_t>(
"ILU_Norm",
"INF_NORM",
666 "Type of norm to use",
667 tuple<string>(
"ONE_NORM",
"TWO_NORM",
"INF_NORM"),
668 tuple<string>(
"1-norm",
"2-norm",
"inf-norm"),
669 tuple<SLU::norm_t>(SLU::ONE_NORM,SLU::TWO_NORM,SLU::INF_NORM),
672 setStringToIntegralParameter<SLU::milu_t>(
"ILU_MILU",
"SILU",
673 "Type of modified ILU to use",
674 tuple<string>(
"SILU",
"SMILU_1",
"SMILU_2",
"SMILU_3"),
675 tuple<string>(
"Regular ILU",
"MILU 1",
"MILU 2",
"MILU 3"),
676 tuple<SLU::milu_t>(SLU::SILU,SLU::SMILU_1,SLU::SMILU_2,
680 pl->set(
"ILU_FillTol", 0.01,
"ILUT fill tolerance");
682 pl->set(
"ILU_Flag",
false,
"ILU flag: if true, run ILU routines");
691 template <
class Matrix,
class Vector>
697 #ifdef HAVE_AMESOS2_TIMERS 698 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
702 if( current_phase == SYMBFACT )
return false;
705 if( data_.A.Store != NULL ){
706 SLU::Destroy_SuperMatrix_Store( &(data_.A) );
707 data_.A.Store = NULL;
712 nzvals_.resize(this->globalNumNonZeros_);
713 rowind_.resize(this->globalNumNonZeros_);
714 colptr_.resize(this->globalNumCols_ + 1);
719 #ifdef HAVE_AMESOS2_TIMERS 720 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
723 TEUCHOS_TEST_FOR_EXCEPTION( this->rowIndexBase_ != this->columnIndexBase_,
725 "Row and column maps have different indexbase ");
728 nzvals_(), rowind_(),
729 colptr_(), nnz_ret,
ROOTED,
731 this->rowIndexBase_);
735 SLU::Dtype_t dtype = type_map::dtype;
738 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
740 "Did not get the expected number of non-zero vals");
742 function_map::create_CompCol_Matrix( &(data_.A),
743 this->globalNumRows_, this->globalNumCols_,
748 SLU::SLU_NC, dtype, SLU::SLU_GE);
755 template<
class Matrix,
class Vector>
761 #endif // AMESOS2_SUPERLU_DEF_HPP Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
Superlu specific solve.
Definition: Amesos2_Superlu_def.hpp:365
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlu_def.hpp:693
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
int numericFactorization_impl()
Superlu specific numeric factorization.
Definition: Amesos2_Superlu_def.hpp:232
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
~Superlu()
Destructor.
Definition: Amesos2_Superlu_def.hpp:108
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlu_def.hpp:569
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Amesos2 Superlu declarations.
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Superlu(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Superlu_def.hpp:67
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
std::string description() const
Returns a short description of this Solver.
Definition: Amesos2_Superlu_def.hpp:131
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superlu_def.hpp:476
Definition: Amesos2_TypeDecl.hpp:127
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_Superlu_def.hpp:487
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superlu_def.hpp:185
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
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using Superlu.
Definition: Amesos2_Superlu_def.hpp:210
Amesos2 interface to the SuperLU package.
Definition: Amesos2_Superlu_decl.hpp:73