52 #ifndef AMESOS2_SUPERLUMT_DEF_HPP 53 #define AMESOS2_SUPERLUMT_DEF_HPP 55 #include <Teuchos_Tuple.hpp> 56 #include <Teuchos_ParameterList.hpp> 57 #include <Teuchos_StandardParameterEntryValidators.hpp> 76 int sp_colorder(SuperMatrix*,
int*,superlumt_options_t*,SuperMatrix*);
96 template <
class Matrix,
class Vector>
98 Teuchos::RCP<Vector> X,
99 Teuchos::RCP<const Vector> B )
105 Teuchos::RCP<Teuchos::ParameterList> default_params
109 data_.options.lwork = 0;
113 data_.options.perm_r = data_.perm_r.getRawPtr();
114 data_.options.perm_c = data_.perm_c.getRawPtr();
119 data_.options.refact = SLUMT::NO;
120 data_.equed = SLUMT::NOEQUIL;
121 data_.rowequ =
false;
122 data_.colequ =
false;
124 data_.A.Store = NULL;
125 data_.AC.Store = NULL;
126 data_.BX.Store = NULL;
127 data_.L.Store = NULL;
128 data_.U.Store = NULL;
130 data_.stat.ops = NULL;
134 template <
class Matrix,
class Vector>
144 if ( data_.A.Store != NULL ){
146 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
150 if( data_.AC.Store != NULL ){
151 SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) );
154 if ( data_.L.Store != NULL ){
155 SLUMT::D::Destroy_SuperNode_SCP( &(data_.L) );
156 SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
159 free( data_.options.etree );
160 free( data_.options.colcnt_h );
161 free( data_.options.part_super_h );
166 if ( data_.BX.Store != NULL ){
172 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
175 if ( data_.stat.ops != NULL )
176 SLUMT::D::StatFree( &(data_.stat) );
179 template<
class Matrix,
class Vector>
184 int perm_spec = data_.options.ColPerm;
185 if( perm_spec != SLUMT::MY_PERMC && this->root_ ){
186 #ifdef HAVE_AMESOS2_TIMERS 187 Teuchos::TimeMonitor preOrderTimer(this->timers_.preOrderTime_);
190 SLUMT::S::get_perm_c(perm_spec, &(data_.A), data_.perm_c.getRawPtr());
199 template <
class Matrix,
class Vector>
207 data_.options.refact = SLUMT::NO;
209 if( this->status_.numericFactorizationDone() && this->root_ ){
214 SLUMT::D::Destroy_SuperNode_Matrix( &(data_.L) );
215 SLUMT::D::Destroy_CompCol_NCP( &(data_.U) );
216 data_.L.Store = NULL;
217 data_.U.Store = NULL;
224 template <
class Matrix,
class Vector>
230 #ifdef HAVE_AMESOS2_DEBUG 231 const int nprocs = data_.options.nprocs;
232 TEUCHOS_TEST_FOR_EXCEPTION( nprocs <= 0,
233 std::invalid_argument,
234 "The number of threads to spawn should be greater than 0." );
241 if( data_.options.fact == SLUMT::EQUILIBRATE ){
242 magnitude_type rowcnd, colcnd, amax;
245 function_map::gsequ(&(data_.A), data_.R.getRawPtr(),
246 data_.C.getRawPtr(), &rowcnd, &colcnd,
248 TEUCHOS_TEST_FOR_EXCEPTION( info != 0,
250 "SuperLU_MT gsequ returned with status " << info );
252 function_map::laqgs(&(data_.A), data_.R.getRawPtr(),
253 data_.C.getRawPtr(), rowcnd, colcnd,
254 amax, &(data_.equed));
256 data_.rowequ = (data_.equed == SLUMT::ROW) || (data_.equed == SLUMT::BOTH);
257 data_.colequ = (data_.equed == SLUMT::COL) || (data_.equed == SLUMT::BOTH);
259 data_.options.fact = SLUMT::DOFACT;
263 if( data_.AC.Store != NULL ){
264 SLUMT::D::Destroy_CompCol_Permuted( &(data_.AC) );
265 if( data_.options.refact == SLUMT::NO ){
266 free( data_.options.etree );
267 free( data_.options.colcnt_h );
268 free( data_.options.part_super_h );
270 data_.AC.Store = NULL;
274 SLUMT::sp_colorder(&(data_.A), data_.perm_c.getRawPtr(),
275 &(data_.options), &(data_.AC));
279 const int n = as<int>(this->globalNumCols_);
280 if( data_.stat.ops != NULL ){ SLUMT::D::StatFree( &(data_.stat) ); data_.stat.ops = NULL; }
281 SLUMT::D::StatAlloc(n, data_.options.nprocs, data_.options.panel_size,
282 data_.options.relax, &(data_.stat));
283 SLUMT::D::StatInit(n, data_.options.nprocs, &(data_.stat));
287 #ifdef HAVE_AMESOS2_TIMERS 288 Teuchos::TimeMonitor numFactTimer( this->timers_.numFactTime_ );
291 #ifdef HAVE_AMESOS2_VERBOSE_DEBUG 292 std::cout <<
"SuperLU_MT:: Before numeric factorization" << std::endl;
293 std::cout <<
"nzvals_ : " << nzvals_.toString() << std::endl;
294 std::cout <<
"rowind_ : " << rowind_.toString() << std::endl;
295 std::cout <<
"colptr_ : " << colptr_.toString() << std::endl;
298 function_map::gstrf(&(data_.options), &(data_.AC),
299 data_.perm_r.getRawPtr(), &(data_.L), &(data_.U),
300 &(data_.stat), &info);
304 this->setNnzLU( as<size_t>(((SLUMT::SCformat*)data_.L.Store)->nnz +
305 ((SLUMT::NCformat*)data_.U.Store)->nnz) );
309 const global_size_type info_st = as<global_size_type>(info);
310 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st <= this->globalNumCols_),
312 "Factorization complete, but matrix is singular. Division by zero eminent");
313 TEUCHOS_TEST_FOR_EXCEPTION( (info_st > 0) && (info_st > this->globalNumCols_),
315 "Memory allocation failure in SuperLU_MT factorization");
319 data_.options.fact = SLUMT::FACTORED;
320 data_.options.refact = SLUMT::YES;
323 Teuchos::broadcast(*(this->getComm()),0,&info);
328 template <
class Matrix,
class Vector>
335 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
336 const size_t nrhs = X->getGlobalNumVectors();
338 Teuchos::Array<slu_type> bxvals_(ld_rhs * nrhs);
339 size_t ldbx_ = as<size_t>(ld_rhs);
342 #ifdef HAVE_AMESOS2_TIMERS 343 Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
344 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
348 slu_type>::do_get(B, bxvals_(),
356 if( data_.BX.Store != NULL ){
357 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.BX) );
358 data_.BX.Store = NULL;
362 #ifdef HAVE_AMESOS2_TIMERS 363 Teuchos::TimeMonitor convTimer( this->timers_.vecConvTime_ );
365 SLUMT::Dtype_t dtype = type_map::dtype;
366 function_map::create_Dense_Matrix(&(data_.BX), as<int>(ld_rhs), as<int>(nrhs),
367 bxvals_.getRawPtr(), as<int>(ldbx_),
368 SLUMT::SLU_DN, dtype, SLUMT::SLU_GE);
371 if( data_.options.trans == SLUMT::NOTRANS ){
374 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
375 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
377 }
else if( data_.colequ ){
379 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
380 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
385 #ifdef HAVE_AMESOS2_TIMERS 386 Teuchos::TimeMonitor solveTimer( this->timers_.solveTime_ );
389 function_map::gstrs(data_.options.trans, &(data_.L),
390 &(data_.U), data_.perm_r.getRawPtr(),
391 data_.perm_c.getRawPtr(), &(data_.BX),
392 &(data_.stat), &info);
397 Teuchos::broadcast(*(this->getComm()),0,&info);
399 TEUCHOS_TEST_FOR_EXCEPTION( info < 0,
401 "Argument " << -info <<
" to gstrs had an illegal value" );
404 if( data_.options.trans == SLUMT::NOTRANS ){
407 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.C(),
408 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
410 }
else if( data_.rowequ ){
412 Util::scale(bxvals_(), as<size_t>(ld_rhs), ldbx_, data_.R(),
413 SLUMT::slu_mt_mult<slu_type,magnitude_type>());
418 #ifdef HAVE_AMESOS2_TIMERS 419 Teuchos::TimeMonitor redistTimer( this->timers_.vecRedistTime_ );
430 template <
class Matrix,
class Vector>
437 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
441 template <
class Matrix,
class Vector>
447 using Teuchos::getIntegralValue;
448 using Teuchos::ParameterEntryValidator;
450 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
453 data_.options.nprocs = parameterList->get<
int>(
"nprocs", 1);
455 data_.options.trans = this->control_.useTranspose_ ? SLUMT::TRANS : SLUMT::NOTRANS;
457 if( parameterList->isParameter(
"trans") ){
458 RCP<const ParameterEntryValidator> trans_validator = valid_params->getEntry(
"trans").validator();
459 parameterList->getEntry(
"trans").setValidator(trans_validator);
461 data_.options.trans = getIntegralValue<SLUMT::trans_t>(*parameterList,
"trans");
464 data_.options.panel_size = parameterList->get<
int>(
"panel_size", SLUMT::D::sp_ienv(1));
466 data_.options.relax = parameterList->get<
int>(
"relax", SLUMT::D::sp_ienv(2));
468 const bool equil = parameterList->get<
bool>(
"Equil",
true);
469 data_.options.fact = equil ? SLUMT::EQUILIBRATE : SLUMT::DOFACT;
471 const bool symmetric_mode = parameterList->get<
bool>(
"SymmetricMode",
false);
472 data_.options.SymmetricMode = symmetric_mode ? SLUMT::YES : SLUMT::NO;
474 const bool print_stat = parameterList->get<
bool>(
"PrintStat",
false);
475 data_.options.PrintStat = print_stat ? SLUMT::YES : SLUMT::NO;
477 data_.options.diag_pivot_thresh = parameterList->get<
double>(
"diag_pivot_thresh", 1.0);
479 if( parameterList->isParameter(
"ColPerm") ){
480 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry(
"ColPerm").validator();
481 parameterList->getEntry(
"ColPerm").setValidator(colperm_validator);
483 data_.options.ColPerm = getIntegralValue<SLUMT::colperm_t>(*parameterList,
"ColPerm");
487 data_.options.usepr = SLUMT::NO;
491 template <
class Matrix,
class Vector>
492 Teuchos::RCP<const Teuchos::ParameterList>
496 using Teuchos::tuple;
497 using Teuchos::ParameterList;
498 using Teuchos::EnhancedNumberValidator;
499 using Teuchos::setStringToIntegralParameter;
500 using Teuchos::stringToIntegralParameterEntryValidator;
502 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
504 if( is_null(valid_params) ){
505 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
507 Teuchos::RCP<EnhancedNumberValidator<int> > nprocs_validator
508 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
509 nprocs_validator->setMin(1);
510 pl->set(
"nprocs", 1,
"The number of processors to factorize with", nprocs_validator);
512 setStringToIntegralParameter<SLUMT::trans_t>(
"trans",
"NOTRANS",
513 "Solve for the transpose system or not",
514 tuple<string>(
"TRANS",
"NOTRANS",
"CONJ"),
515 tuple<string>(
"Solve with transpose",
516 "Do not solve with transpose",
517 "Solve with the conjugate transpose"),
518 tuple<SLUMT::trans_t>(SLUMT::TRANS,
523 Teuchos::RCP<EnhancedNumberValidator<int> > panel_size_validator
524 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
525 panel_size_validator->setMin(1);
526 pl->set(
"panel_size", SLUMT::D::sp_ienv(1),
527 "Specifies the number of consecutive columns to be treated as a unit of task.",
528 panel_size_validator);
530 Teuchos::RCP<EnhancedNumberValidator<int> > relax_validator
531 = Teuchos::rcp(
new EnhancedNumberValidator<int>() );
532 relax_validator->setMin(1);
533 pl->set(
"relax", SLUMT::D::sp_ienv(2),
534 "Specifies the number of columns to be grouped as a relaxed supernode",
537 pl->set(
"Equil",
true,
"Whether to equilibrate the system before solve");
539 Teuchos::RCP<EnhancedNumberValidator<double> > diag_pivot_thresh_validator
540 = Teuchos::rcp(
new EnhancedNumberValidator<double>(0.0, 1.0) );
541 pl->set(
"diag_pivot_thresh", 1.0,
542 "Specifies the threshold used for a diagonal entry to be an acceptable pivot",
543 diag_pivot_thresh_validator);
546 setStringToIntegralParameter<SLUMT::colperm_t>(
"ColPerm",
"COLAMD",
547 "Specifies how to permute the columns of the " 548 "matrix for sparsity preservation",
549 tuple<string>(
"NATURAL",
553 tuple<string>(
"Natural ordering",
554 "Minimum degree ordering on A^T + A",
555 "Minimum degree ordering on A^T A",
556 "Approximate minimum degree column ordering"),
557 tuple<SLUMT::colperm_t>(SLUMT::NATURAL,
558 SLUMT::MMD_AT_PLUS_A,
563 pl->set(
"SymmetricMode",
false,
"Specifies whether to use the symmetric mode");
568 pl->set(
"PrintStat",
false,
"Specifies whether to print the solver's statistics");
577 template <
class Matrix,
class Vector>
583 #ifdef HAVE_AMESOS2_TIMERS 584 Teuchos::TimeMonitor convTimer( this->timers_.mtxConvTime_ );
587 if( current_phase == SYMBFACT )
return false;
590 if( data_.A.Store != NULL ){
591 SLUMT::D::Destroy_SuperMatrix_Store( &(data_.A) );
592 data_.A.Store = NULL;
596 nzvals_.resize(this->globalNumNonZeros_);
597 rowind_.resize(this->globalNumNonZeros_);
598 colptr_.resize(this->globalNumCols_ + 1);
603 #ifdef HAVE_AMESOS2_TIMERS 604 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
610 nzvals_, rowind_, colptr_,
615 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != as<int>(this->globalNumNonZeros_),
617 "rank=0 failed to get all nonzero values in getCcs()");
619 SLUMT::Dtype_t dtype = type_map::dtype;
620 function_map::create_CompCol_Matrix(&(data_.A),
621 as<int>(this->globalNumRows_),
622 as<int>(this->globalNumCols_),
628 dtype, SLUMT::SLU_GE);
630 TEUCHOS_TEST_FOR_EXCEPTION( data_.A.Store == NULL,
632 "Failed to allocate matrix store" );
639 template<
class Matrix,
class Vector>
645 #endif // AMESOS2_SUPERLUMT_DEF_HPP bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superlumt_def.hpp:432
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
void scale(ArrayView< Scalar1 > vals, size_t l, size_t ld, ArrayView< Scalar2 > s)
Scales a 1-D representation of a multivector.
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
SuperLU_MT specific solve.
Definition: Amesos2_Superlumt_def.hpp:330
Amesos2 interface to the Multi-threaded version of SuperLU.
Definition: Amesos2_Superlumt_decl.hpp:77
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using SuperLU_MT.
Definition: Amesos2_Superlumt_def.hpp:201
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:243
Definition: Amesos2_TypeDecl.hpp:142
Utility functions for Amesos2.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this->setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:310
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superlumt_def.hpp:493
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Definition: Amesos2_Superlumt_def.hpp:63
A generic helper class for getting a CCS representation of a Matrix.
Definition: Amesos2_Util.hpp:580
Superlumt(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Superlumt_def.hpp:97
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superlumt_def.hpp:181
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
int numericFactorization_impl()
SuperLU_MT specific numeric factorization.
Definition: Amesos2_Superlumt_def.hpp:226
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_Superlumt_def.hpp:443
~Superlumt()
Destructor.
Definition: Amesos2_Superlumt_def.hpp:135
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:278
Definition: Amesos2_TypeDecl.hpp:127
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
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal structures.
Definition: Amesos2_Superlumt_def.hpp:579