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( info2 != 0,
270 "SuperLU gsequ returned with status " << info2 );
273 function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
274 data_.C.getRawPtr(), rowcnd, colcnd,
275 amax, &(data_.equed));
284 SLU::sp_preorder(&(data_.options), &(data_.A), data_.perm_c.getRawPtr(),
285 data_.etree.getRawPtr(), &(data_.AC));
288 #ifdef HAVE_AMESOS2_TIMERS 289 Teuchos::TimeMonitor numFactTimer(this->
timers_.numFactTime_);
292 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG 293 std::cout <<
"Superlu:: Before numeric factorization" << std::endl;
294 std::cout <<
"nzvals_ : " <<
nzvals_.toString() << std::endl;
295 std::cout <<
"rowind_ : " <<
rowind_.toString() << std::endl;
296 std::cout <<
"colptr_ : " <<
colptr_.toString() << std::endl;
299 if(ILU_Flag_==
false) {
300 function_map::gstrf(&(data_.options), &(data_.AC),
301 data_.relax, data_.panel_size, data_.etree.getRawPtr(),
302 NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
303 &(data_.L), &(data_.U), &(data_.stat), &info);
306 function_map::gsitrf(&(data_.options), &(data_.AC),
307 data_.relax, data_.panel_size, data_.etree.getRawPtr(),
308 NULL, 0, data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
309 &(data_.L), &(data_.U), &(data_.stat), &info);
314 SLU::Destroy_CompCol_Permuted( &(data_.AC) );
317 this->
setNnzLU( as<size_t>(((SLU::SCformat*)data_.L.Store)->nnz +
318 ((SLU::NCformat*)data_.U.Store)->nnz) );
322 Teuchos::broadcast(*(this->
matrixA_->getComm()), 0, &info);
324 global_size_type info_st = as<global_size_type>(info);
325 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->
globalNumCols_),
327 "Factorization complete, but matrix is singular. Division by zero eminent");
328 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
330 "Memory allocation failure in Superlu factorization");
332 data_.options.Fact = SLU::FACTORED;
333 same_symbolic_ =
true;
339 template <
class Matrix,
class Vector>
346 const global_size_type ld_rhs = this->
root_ ? X->getGlobalLength() : 0;
347 const size_t nrhs = X->getGlobalNumVectors();
349 const size_t val_store_size = as<size_t>(ld_rhs * nrhs);
350 Teuchos::Array<slu_type> xValues(val_store_size);
351 Teuchos::Array<slu_type> bValues(val_store_size);
354 #ifdef HAVE_AMESOS2_TIMERS 355 Teuchos::TimeMonitor mvConvTimer(this->
timers_.vecConvTime_);
356 Teuchos::TimeMonitor redistTimer( this->
timers_.vecRedistTime_ );
359 slu_type>::do_get(B, bValues(),
366 magnitude_type rpg, rcond;
368 data_.ferr.resize(nrhs);
369 data_.berr.resize(nrhs);
372 #ifdef HAVE_AMESOS2_TIMERS 373 Teuchos::TimeMonitor mvConvTimer(this->
timers_.vecConvTime_);
375 SLU::Dtype_t dtype = type_map::dtype;
376 int i_ld_rhs = as<int>(ld_rhs);
377 function_map::create_Dense_Matrix(&(data_.B), i_ld_rhs, as<int>(nrhs),
378 bValues.getRawPtr(), i_ld_rhs,
379 SLU::SLU_DN, dtype, SLU::SLU_GE);
380 function_map::create_Dense_Matrix(&(data_.X), i_ld_rhs, as<int>(nrhs),
381 xValues.getRawPtr(), i_ld_rhs,
382 SLU::SLU_DN, dtype, SLU::SLU_GE);
389 #ifdef HAVE_AMESOS2_TIMERS 390 Teuchos::TimeMonitor solveTimer(this->
timers_.solveTime_);
393 if(ILU_Flag_==
false) {
394 function_map::gssvx(&(data_.options), &(data_.A),
395 data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
396 data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
397 data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
398 &(data_.X), &rpg, &rcond, data_.ferr.getRawPtr(),
399 data_.berr.getRawPtr(), &(data_.mem_usage), &(data_.stat), &ierr);
402 function_map::gsisx(&(data_.options), &(data_.A),
403 data_.perm_c.getRawPtr(), data_.perm_r.getRawPtr(),
404 data_.etree.getRawPtr(), &(data_.equed), data_.R.getRawPtr(),
405 data_.C.getRawPtr(), &(data_.L), &(data_.U), NULL, 0, &(data_.B),
406 &(data_.X), &rpg, &rcond, &(data_.mem_usage), &(data_.stat), &ierr);
412 SLU::Destroy_SuperMatrix_Store( &(data_.X) );
413 SLU::Destroy_SuperMatrix_Store( &(data_.B) );
414 data_.X.Store = NULL;
415 data_.B.Store = NULL;
419 Teuchos::broadcast(*(this->
getComm()), 0, &ierr);
421 global_size_type ierr_st = as<global_size_type>(ierr);
422 TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
423 std::invalid_argument,
424 "Argument " << -ierr <<
" to SuperLU xgssvx had illegal value" );
425 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st <= this->
globalNumCols_,
427 "Factorization complete, but U is exactly singular" );
428 TEUCHOS_TEST_FOR_EXCEPTION( ierr > 0 && ierr_st > this->globalNumCols_ + 1,
430 "SuperLU allocated " << ierr - this->globalNumCols_ <<
" bytes of " 431 "memory before allocation failure occured." );
435 #ifdef HAVE_AMESOS2_TIMERS 436 Teuchos::TimeMonitor redistTimer(this->
timers_.vecRedistTime_);
450 template <
class Matrix,
class Vector>
457 return( this->
matrixA_->getGlobalNumRows() == this->
matrixA_->getGlobalNumCols() );
461 template <
class Matrix,
class Vector>
466 using Teuchos::getIntegralValue;
467 using Teuchos::ParameterEntryValidator;
471 ILU_Flag_ = parameterList->get<
bool>(
"ILU_Flag",
false);
473 SLU::ilu_set_default_options(&(data_.options));
475 data_.options.PrintStat = SLU::NO;
478 data_.options.Trans = this->
control_.useTranspose_ ? SLU::TRANS : SLU::NOTRANS;
480 if( parameterList->isParameter(
"Trans") ){
481 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"Trans").validator();
482 parameterList->getEntry(
"Trans").setValidator(trans_validator);
484 data_.options.Trans = getIntegralValue<SLU::trans_t>(*parameterList,
"Trans");
487 if( parameterList->isParameter(
"IterRefine") ){
488 RCP<const ParameterEntryValidator> refine_validator = valid_params->getEntry(
"IterRefine").validator();
489 parameterList->getEntry(
"IterRefine").setValidator(refine_validator);
491 data_.options.IterRefine = getIntegralValue<SLU::IterRefine_t>(*parameterList,
"IterRefine");
494 if( parameterList->isParameter(
"ColPerm") ){
495 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry(
"ColPerm").validator();
496 parameterList->getEntry(
"ColPerm").setValidator(colperm_validator);
498 data_.options.ColPerm = getIntegralValue<SLU::colperm_t>(*parameterList,
"ColPerm");
501 data_.options.DiagPivotThresh = parameterList->get<
double>(
"DiagPivotThresh", 1.0);
503 bool equil = parameterList->get<
bool>(
"Equil",
true);
504 data_.options.Equil = equil ? SLU::YES : SLU::NO;
506 bool symmetric_mode = parameterList->get<
bool>(
"SymmetricMode",
false);
507 data_.options.SymmetricMode = symmetric_mode ? SLU::YES : SLU::NO;
510 if( parameterList->isParameter(
"RowPerm") ){
511 RCP<const ParameterEntryValidator> rowperm_validator = valid_params->getEntry(
"RowPerm").validator();
512 parameterList->getEntry(
"RowPerm").setValidator(rowperm_validator);
513 data_.options.RowPerm = getIntegralValue<SLU::rowperm_t>(*parameterList,
"RowPerm");
522 data_.options.ILU_DropTol = parameterList->get<
double>(
"ILU_DropTol", 0.0001);
524 data_.options.ILU_FillFactor = parameterList->get<
double>(
"ILU_FillFactor", 10.0);
526 if( parameterList->isParameter(
"ILU_Norm") ) {
527 RCP<const ParameterEntryValidator> norm_validator = valid_params->getEntry(
"ILU_Norm").validator();
528 parameterList->getEntry(
"ILU_Norm").setValidator(norm_validator);
529 data_.options.ILU_Norm = getIntegralValue<SLU::norm_t>(*parameterList,
"ILU_Norm");
532 if( parameterList->isParameter(
"ILU_MILU") ) {
533 RCP<const ParameterEntryValidator> milu_validator = valid_params->getEntry(
"ILU_MILU").validator();
534 parameterList->getEntry(
"ILU_MILU").setValidator(milu_validator);
535 data_.options.ILU_MILU = getIntegralValue<SLU::milu_t>(*parameterList,
"ILU_MILU");
538 data_.options.ILU_FillTol = parameterList->get<
double>(
"ILU_FillTol", 0.01);
543 template <
class Matrix,
class Vector>
544 Teuchos::RCP<const Teuchos::ParameterList>
548 using Teuchos::tuple;
549 using Teuchos::ParameterList;
550 using Teuchos::EnhancedNumberValidator;
551 using Teuchos::setStringToIntegralParameter;
552 using Teuchos::stringToIntegralParameterEntryValidator;
554 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
556 if( is_null(valid_params) ){
557 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
559 setStringToIntegralParameter<SLU::trans_t>(
"Trans",
"NOTRANS",
560 "Solve for the transpose system or not",
561 tuple<string>(
"TRANS",
"NOTRANS",
"CONJ"),
562 tuple<string>(
"Solve with transpose",
563 "Do not solve with transpose",
564 "Solve with the conjugate transpose"),
565 tuple<SLU::trans_t>(SLU::TRANS,
570 setStringToIntegralParameter<SLU::IterRefine_t>(
"IterRefine",
"NOREFINE",
571 "Type of iterative refinement to use",
572 tuple<string>(
"NOREFINE",
"SLU_SINGLE",
"SLU_DOUBLE"),
573 tuple<string>(
"Do not use iterative refinement",
574 "Do single iterative refinement",
575 "Do double iterative refinement"),
576 tuple<SLU::IterRefine_t>(SLU::NOREFINE,
582 setStringToIntegralParameter<SLU::colperm_t>(
"ColPerm",
"COLAMD",
583 "Specifies how to permute the columns of the " 584 "matrix for sparsity preservation",
585 tuple<string>(
"NATURAL",
"MMD_AT_PLUS_A",
586 "MMD_ATA",
"COLAMD"),
587 tuple<string>(
"Natural ordering",
588 "Minimum degree ordering on A^T + A",
589 "Minimum degree ordering on A^T A",
590 "Approximate minimum degree column ordering"),
591 tuple<SLU::colperm_t>(SLU::NATURAL,
597 Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
598 = Teuchos::rcp(
new EnhancedNumberValidator<double>(0.0, 1.0) );
599 pl->set(
"DiagPivotThresh", 1.0,
600 "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
601 diag_pivot_thresh_validator);
603 pl->set(
"Equil",
true,
"Whether to equilibrate the system before solve");
605 pl->set(
"SymmetricMode",
false,
606 "Specifies whether to use the symmetric mode. " 607 "Gives preference to diagonal pivots and uses " 608 "an (A^T + A)-based column permutation.");
612 setStringToIntegralParameter<SLU::rowperm_t>(
"RowPerm",
"LargeDiag",
613 "Type of row permutation strategy to use",
614 tuple<string>(
"NOROWPERM",
"LargeDiag",
"MY_PERMR"),
615 tuple<string>(
"Use natural ordering",
616 "Use weighted bipartite matching algorithm",
617 "Use the ordering given in perm_r input"),
618 tuple<SLU::rowperm_t>(SLU::NOROWPERM,
637 pl->set(
"ILU_DropTol", 0.0001,
"ILUT drop tolerance");
639 pl->set(
"ILU_FillFactor", 10.0,
"ILUT fill factor");
641 setStringToIntegralParameter<SLU::norm_t>(
"ILU_Norm",
"INF_NORM",
642 "Type of norm to use",
643 tuple<string>(
"ONE_NORM",
"TWO_NORM",
"INF_NORM"),
644 tuple<string>(
"1-norm",
"2-norm",
"inf-norm"),
645 tuple<SLU::norm_t>(SLU::ONE_NORM,SLU::TWO_NORM,SLU::INF_NORM),
648 setStringToIntegralParameter<SLU::milu_t>(
"ILU_MILU",
"SILU",
649 "Type of modified ILU to use",
650 tuple<string>(
"SILU",
"SMILU_1",
"SMILU_2",
"SMILU_3"),
651 tuple<string>(
"Regular ILU",
"MILU 1",
"MILU 2",
"MILU 3"),
652 tuple<SLU::milu_t>(SLU::SILU,SLU::SMILU_1,SLU::SMILU_2,
656 pl->set(
"ILU_FillTol", 0.01,
"ILUT fill tolerance");
658 pl->set(
"ILU_Flag",
false,
"ILU flag: if true, run ILU routines");
667 template <
class Matrix,
class Vector>
673 #ifdef HAVE_AMESOS2_TIMERS 674 Teuchos::TimeMonitor convTimer(this->
timers_.mtxConvTime_);
678 if( current_phase == SYMBFACT )
return false;
681 if( data_.A.Store != NULL ){
682 SLU::Destroy_SuperMatrix_Store( &(data_.A) );
683 data_.A.Store = NULL;
695 #ifdef HAVE_AMESOS2_TIMERS 696 Teuchos::TimeMonitor mtxRedistTimer( this->
timers_.mtxRedistTime_ );
701 "Row and column maps have different indexbase ");
711 SLU::Dtype_t dtype = type_map::dtype;
716 "Did not get the expected number of non-zero vals");
718 function_map::create_CompCol_Matrix( &(data_.A),
719 this->globalNumRows_, this->globalNumCols_,
724 SLU::SLU_NC, dtype, SLU::SLU_GE);
731 template<
class Matrix,
class Vector>
737 #endif // AMESOS2_SUPERLU_DEF_HPP Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
Teuchos::Array< slu_type > nzvals_
Stores the values of the nonzero entries for SuperLU.
Definition: Amesos2_Superlu_decl.hpp:259
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlu_def.hpp:669
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
Teuchos::Array< int > colptr_
Stores the row indices of the nonzero entries.
Definition: Amesos2_Superlu_decl.hpp:263
int numericFactorization_impl()
Superlu specific numeric factorization.
Definition: Amesos2_Superlu_def.hpp:232
global_size_type columnIndexBase_
Index base of column map of matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:488
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::Array< int > rowind_
Stores the location in Ai_ and Aval_ that starts row j.
Definition: Amesos2_Superlu_decl.hpp:261
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
Amesos2 Superlu declarations.
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superlu_def.hpp:452
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
Control control_
Parameters for solving.
Definition: Amesos2_SolverCore_decl.hpp:495
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlu_def.hpp:545
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
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:341
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
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
global_size_type globalNumNonZeros_
Number of global non-zero values in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:482
void setNnzLU(size_t nnz)
Set the number of non-zero values in the and factors.
Definition: Amesos2_SolverCore_decl.hpp:452
Definition: Amesos2_TypeDecl.hpp:127
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_Superlu_def.hpp:463
Timers timers_
Various timing statistics.
Definition: Amesos2_SolverCore_decl.hpp:498
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
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
std::string description() const
Returns a short description of this Solver.
Definition: Amesos2_Superlu_def.hpp:131
global_size_type rowIndexBase_
Index base of rowmap of matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:485
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