52 #ifndef AMESOS2_SUPERLUDIST_DEF_HPP
53 #define AMESOS2_SUPERLUDIST_DEF_HPP
55 #include <Teuchos_Tuple.hpp>
56 #include <Teuchos_StandardParameterEntryValidators.hpp>
57 #include <Teuchos_DefaultMpiComm.hpp>
58 #include <Teuchos_Details_MpiTypeTraits.hpp>
68 template <
class Matrix,
class Vector>
70 Teuchos::RCP<Vector> X,
71 Teuchos::RCP<const Vector> B)
76 , force_symbfact_(false)
77 , is_contiguous_(true)
82 using Teuchos::MpiComm;
83 using Teuchos::outArg;
84 using Teuchos::ParameterList;
85 using Teuchos::parameterList;
88 using Teuchos::rcp_dynamic_cast;
89 using Teuchos::REDUCE_SUM;
90 using Teuchos::reduceAll;
91 typedef global_ordinal_type GO;
92 typedef Tpetra::Map<local_ordinal_type, GO, node_type> map_type;
98 RCP<const Comm<int> > comm = this->
getComm ();
99 const int myRank = comm->getRank ();
100 const int numProcs = comm->getSize ();
102 SLUD::int_t nprow, npcol;
109 RCP<const Comm<int> > matComm = this->
matrixA_->getComm ();
110 TEUCHOS_TEST_FOR_EXCEPTION(
111 matComm.is_null (), std::logic_error,
"Amesos2::Superlustdist "
112 "constructor: The matrix's communicator is null!");
113 RCP<const MpiComm<int> > matMpiComm =
114 rcp_dynamic_cast<
const MpiComm<int> > (matComm);
119 TEUCHOS_TEST_FOR_EXCEPTION(
120 matMpiComm.is_null (), std::logic_error,
"Amesos2::Superlustdist "
121 "constructor: The matrix's communicator is not an MpiComm!");
122 TEUCHOS_TEST_FOR_EXCEPTION(
123 matMpiComm->getRawMpiComm ().is_null (), std::logic_error,
"Amesos2::"
124 "Superlustdist constructor: The matrix's communicator claims to be a "
125 "Teuchos::MpiComm<int>, but its getRawPtrComm() method returns "
126 "Teuchos::null! This means that the underlying MPI_Comm doesn't even "
127 "exist, which likely implies that the Teuchos::MpiComm was constructed "
128 "incorrectly. It means something different than if the MPI_Comm were "
130 MPI_Comm rawMpiComm = (* (matMpiComm->getRawMpiComm ())) ();
131 data_.mat_comm = rawMpiComm;
137 SLUD::superlu_gridinit(data_.mat_comm, nprow, npcol, &(data_.grid));
146 SLUD::set_default_options_dist(&data_.options);
148 RCP<ParameterList> default_params =
153 data_.options.Fact = SLUD::DOFACT;
154 data_.equed = SLUD::NOEQUIL;
155 data_.options.SolveInitialized = SLUD::NO;
156 data_.options.RefineInitialized = SLUD::NO;
157 data_.rowequ =
false;
158 data_.colequ =
false;
161 data_.largediag_mc64_job = 4;
171 data_.symb_comm = MPI_COMM_NULL;
176 data_.domains = (
int) ( pow(2.0, floor(log10((
double)nprow*npcol)/log10(2.0))) );
178 const int color = (myRank < data_.domains) ? 0 : MPI_UNDEFINED;
179 MPI_Comm_split (data_.mat_comm, color, myRank, &(data_.symb_comm));
191 int myProcParticipates = 0;
192 if (myRank < nprow * npcol) {
194 myProcParticipates = 1;
199 int numParticipatingProcs = 0;
200 reduceAll<int, int> (*comm, REDUCE_SUM, myProcParticipates,
201 outArg (numParticipatingProcs));
202 TEUCHOS_TEST_FOR_EXCEPTION(
203 this->globalNumRows_ != 0 && numParticipatingProcs == 0,
204 std::logic_error,
"Amesos2::Superludist constructor: The matrix has "
205 << this->globalNumRows_ <<
" > 0 global row(s), but no processes in the "
206 "communicator want to participate in its factorization! nprow = "
207 << nprow <<
" and npcol = " << npcol <<
".");
210 size_t myNumRows = 0;
213 const GO quotient = (numParticipatingProcs == 0) ? static_cast<GO> (0) :
214 GNR /
static_cast<GO
> (numParticipatingProcs);
216 GNR - quotient *
static_cast<GO
> (numParticipatingProcs);
217 const GO lclNumRows = (
static_cast<GO
> (myRank) < remainder) ?
218 (quotient +
static_cast<GO
> (1)) : quotient;
219 myNumRows =
static_cast<size_t> (lclNumRows);
223 const GO indexBase = this->
matrixA_->getRowMap ()->getIndexBase ();
225 rcp (
new map_type (this->globalNumRows_, myNumRows, indexBase, comm));
231 data_.A.Store = NULL;
232 function_map::LUstructInit(this->globalNumRows_, this->globalNumCols_, &(data_.lu));
233 SLUD::PStatInit(&(data_.stat));
236 data_.scale_perm.perm_r = data_.perm_r.getRawPtr();
237 data_.scale_perm.perm_c = data_.perm_c.getRawPtr();
241 template <
class Matrix,
class Vector>
253 if ( this->status_.getNumPreOrder() > 0 ){
255 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
256 SUPERLU_FREE( data_.sizes );
257 SUPERLU_FREE( data_.fstVtxSep );
260 free( data_.fstVtxSep );
266 if( data_.A.Store != NULL ){
267 SLUD::Destroy_SuperMatrix_Store_dist( &(data_.A) );
271 if ( this->status_.getNumNumericFact() > 0 ){
272 function_map::Destroy_LU(this->globalNumRows_, &(data_.grid), &(data_.lu));
274 function_map::LUstructFree(&(data_.lu));
280 if ( this->status_.symbolicFactorizationDone() &&
281 !this->status_.numericFactorizationDone() ){
282 if ( data_.pslu_freeable.xlsub != NULL ){
283 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
284 SUPERLU_FREE( data_.pslu_freeable.xlsub );
285 SUPERLU_FREE( data_.pslu_freeable.lsub );
287 free( data_.pslu_freeable.xlsub );
288 free( data_.pslu_freeable.lsub );
291 if ( data_.pslu_freeable.xusub != NULL ){
292 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
293 SUPERLU_FREE( data_.pslu_freeable.xusub );
294 SUPERLU_FREE( data_.pslu_freeable.usub );
296 free( data_.pslu_freeable.xusub );
297 free( data_.pslu_freeable.usub );
300 if ( data_.pslu_freeable.supno_loc != NULL ){
301 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
302 SUPERLU_FREE( data_.pslu_freeable.supno_loc );
303 SUPERLU_FREE( data_.pslu_freeable.xsup_beg_loc );
304 SUPERLU_FREE( data_.pslu_freeable.xsup_end_loc );
306 free( data_.pslu_freeable.supno_loc );
307 free( data_.pslu_freeable.xsup_beg_loc );
308 free( data_.pslu_freeable.xsup_end_loc );
311 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
312 SUPERLU_FREE( data_.pslu_freeable.globToLoc );
314 free( data_.pslu_freeable.globToLoc );
318 SLUD::PStatFree( &(data_.stat) ) ;
323 if ( data_.options.SolveInitialized == SLUD::YES )
324 function_map::SolveFinalize(&(data_.options), &(data_.solve_struct));
328 SLUD::superlu_gridexit(&(data_.grid));
332 if ( data_.symb_comm != MPI_COMM_NULL ) MPI_Comm_free(&(data_.symb_comm));
335 template<
class Matrix,
class Vector>
339 int job = data_.largediag_mc64_job;
342 data_.R1.resize(data_.A.nrow);
343 data_.C1.resize(data_.A.ncol);
346 SLUD::NCformat *GAstore = (SLUD::NCformat*) GA.Store;
347 SLUD::int_t* colptr = GAstore->colptr;
348 SLUD::int_t* rowind = GAstore->rowind;
349 SLUD::int_t nnz = GAstore->nnz;
350 slu_type *a_GA = (slu_type *) GAstore->nzval;
351 MPI_Datatype mpi_dtype = Teuchos::Details::MpiTypeTraits<magnitude_type>::getType();
352 MPI_Datatype mpi_itype = Teuchos::Details::MpiTypeTraits<SLUD::int_t>::getType();
355 if ( !data_.grid.iam ) {
356 iinfo = function_map::ldperm_dist(job, data_.A.nrow, nnz, colptr, rowind, a_GA,
357 data_.perm_r.getRawPtr(), data_.R1.getRawPtr(), data_.C1.getRawPtr());
359 MPI_Bcast( &iinfo, 1, MPI_INT, 0, data_.grid.comm );
361 MPI_Bcast( data_.perm_r.getRawPtr(), data_.A.nrow, mpi_itype, 0, data_.grid.comm );
362 if ( job == 5 && data_.options.Equil ) {
363 MPI_Bcast( data_.R1.getRawPtr(), data_.A.nrow, mpi_dtype, 0, data_.grid.comm );
364 MPI_Bcast( data_.C1.getRawPtr(), data_.A.ncol, mpi_dtype, 0, data_.grid.comm );
368 MPI_Bcast( &iinfo, 1, mpi_int_t, 0, data_.grid.comm );
370 MPI_Bcast( data_.perm_r.getRawPtr(), data_.A.nrow, mpi_itype, 0, data_.grid.comm );
371 if ( job == 5 && data_.options.Equil ) {
372 MPI_Bcast( data_.R1.getRawPtr(), data_.A.nrow, mpi_dtype, 0, data_.grid.comm );
373 MPI_Bcast( data_.C1.getRawPtr(), data_.A.ncol, mpi_dtype, 0, data_.grid.comm );
377 TEUCHOS_TEST_FOR_EXCEPTION( iinfo != 0,
379 "SuperLU_DIST pre-ordering failed to compute row perm with "
380 << iinfo << std::endl);
384 for (SLUD::int_t i = 0; i < data_.A.nrow; ++i) data_.R1[i] = exp(data_.R1[i]);
385 for (SLUD::int_t i = 0; i < data_.A.ncol; ++i) data_.C1[i] = exp(data_.C1[i]);
390 template<
class Matrix,
class Vector>
394 if (data_.options.RowPerm == SLUD::NOROWPERM) {
395 SLUD::int_t slu_rows_ub = Teuchos::as<SLUD::int_t>(this->globalNumRows_);
396 for( SLUD::int_t i = 0; i < slu_rows_ub; ++i ) data_.perm_r[i] = i;
398 else if (data_.options.RowPerm == SLUD::LargeDiag_MC64) {
399 if (!force_symbfact_)
401 return (EXIT_SUCCESS + 1);
411 if( this->status_.getNumPreOrder() > 0 ){
412 #if defined(AMESOS2_ENABLES_SUPERLUDIST_VERSION5_AND_HIGHER)
413 SUPERLU_FREE( data_.sizes );
414 SUPERLU_FREE( data_.fstVtxSep );
417 free( data_.fstVtxSep );
422 #ifdef HAVE_AMESOS2_TIMERS
423 Teuchos::TimeMonitor preOrderTime( this->timers_.preOrderTime_ );
425 info = SLUD::get_perm_c_parmetis( &(data_.A),
426 data_.perm_r.getRawPtr(), data_.perm_c.getRawPtr(),
427 data_.grid.nprow * data_.grid.npcol, data_.domains,
428 &(data_.sizes), &(data_.fstVtxSep),
429 &(data_.grid), &(data_.symb_comm) );
431 TEUCHOS_TEST_FOR_EXCEPTION( info > 0.0,
433 "SuperLU_DIST pre-ordering ran out of memory after allocating "
434 << info <<
" bytes of memory" );
446 template <
class Matrix,
class Vector>
451 if (!force_symbfact_) {
452 if (data_.options.RowPerm == SLUD::LargeDiag_MC64) {
454 return (EXIT_SUCCESS + 1);
462 #ifdef HAVE_AMESOS2_TIMERS
463 Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
466 #if (SUPERLU_DIST_MAJOR_VERSION > 7)
467 info = SLUD::symbfact_dist(&(data_.options), (data_.grid.nprow) * (data_.grid.npcol),
468 data_.domains, &(data_.A), data_.perm_c.getRawPtr(),
469 data_.perm_r.getRawPtr(), data_.sizes,
470 data_.fstVtxSep, &(data_.pslu_freeable),
471 &(data_.grid.comm), &(data_.symb_comm),
475 info = SLUD::symbfact_dist((data_.grid.nprow) * (data_.grid.npcol),
476 data_.domains, &(data_.A), data_.perm_c.getRawPtr(),
477 data_.perm_r.getRawPtr(), data_.sizes,
478 data_.fstVtxSep, &(data_.pslu_freeable),
479 &(data_.grid.comm), &(data_.symb_comm),
483 TEUCHOS_TEST_FOR_EXCEPTION( info > 0.0,
485 "SuperLU_DIST symbolic factorization ran out of memory after"
486 " allocating " << info <<
" bytes of memory" );
488 same_symbolic_ =
false;
489 same_solve_struct_ =
false;
495 template <
class Matrix,
class Vector>
501 SLUD::SuperMatrix GA;
502 bool need_value =
false;
505 if( data_.options.Equil == SLUD::YES ) {
506 SLUD::int_t info = 0;
509 data_.R.resize(this->globalNumRows_);
510 data_.C.resize(this->globalNumCols_);
511 function_map::gsequ_loc(&(data_.A), data_.R.getRawPtr(), data_.C.getRawPtr(),
512 &(data_.rowcnd), &(data_.colcnd), &(data_.amax), &info, &(data_.grid));
515 function_map::laqgs_loc(&(data_.A), data_.R.getRawPtr(), data_.C.getRawPtr(),
516 data_.rowcnd, data_.colcnd, data_.amax,
519 data_.rowequ = (data_.equed == SLUD::ROW) || (data_.equed == SLUD::BOTH);
520 data_.colequ = (data_.equed == SLUD::COL) || (data_.equed == SLUD::BOTH);
523 if (data_.options.RowPerm == SLUD::LargeDiag_MC64) {
526 SLUD::D::pdCompRow_loc_to_CompCol_global(
true, &data_.A, &data_.grid, &GA);
529 computeRowPermutationLargeDiagMC64(GA);
532 force_symbfact_ =
true;
534 symbolicFactorization_impl();
535 force_symbfact_ =
false;
539 if (data_.largediag_mc64_job == 5)
541 SLUD::NRformat_loc *Astore = (SLUD::NRformat_loc*) data_.A.Store;
542 slu_type *a = (slu_type*) Astore->nzval;
543 SLUD::int_t m_loc = Astore->m_loc;
544 SLUD::int_t fst_row = Astore->fst_row;
545 SLUD::int_t i, j, irow = fst_row, icol;
549 SLUD::slu_dist_mult<slu_type, magnitude_type> mult_op;
550 for (j = 0; j < m_loc; ++j) {
551 for (i = rowptr_view_.data()[j]; i < rowptr_view_.data()[j+1]; ++i) {
552 icol = colind_view_.data()[i];
553 a[i] = mult_op(a[i], data_.R1[irow] * data_.C1[icol]);
559 if ( data_.rowequ )
for (i = 0; i < data_.A.nrow; ++i) data_.R[i] *= data_.R1[i];
560 else for (i = 0; i < data_.A.nrow; ++i) data_.R[i] = data_.R1[i];
561 if ( data_.colequ )
for (i = 0; i < data_.A.ncol; ++i) data_.C[i] *= data_.C1[i];
562 else for (i = 0; i < data_.A.ncol; ++i) data_.C[i] = data_.C1[i];
564 data_.rowequ = data_.colequ = 1;
570 size_t nnz_loc = ((SLUD::NRformat_loc*)data_.A.Store)->nnz_loc;
571 for(
size_t i = 0; i < nnz_loc; ++i ) colind_view_(i) = data_.perm_c[colind_view_(i)];
574 if( same_symbolic_ ){
579 #if (SUPERLU_DIST_MAJOR_VERSION > 7)
580 data_.options.Fact = SLUD::SamePattern_SameRowPerm;
581 function_map::pdistribute(&(data_.options),
582 as<SLUD::int_t>(this->globalNumRows_),
583 &(data_.A), &(data_.scale_perm),
584 &(data_.glu_freeable), &(data_.lu),
587 function_map::pdistribute(SLUD::SamePattern_SameRowPerm,
588 as<SLUD::int_t>(this->globalNumRows_),
589 &(data_.A), &(data_.scale_perm),
590 &(data_.glu_freeable), &(data_.lu),
594 #if (SUPERLU_DIST_MAJOR_VERSION > 7)
595 data_.options.Fact = SLUD::DOFACT;
596 function_map::dist_psymbtonum(&(data_.options),
597 as<SLUD::int_t>(this->globalNumRows_),
598 &(data_.A), &(data_.scale_perm),
599 &(data_.pslu_freeable), &(data_.lu),
602 function_map::dist_psymbtonum(SLUD::DOFACT,
603 as<SLUD::int_t>(this->globalNumRows_),
604 &(data_.A), &(data_.scale_perm),
605 &(data_.pslu_freeable), &(data_.lu),
611 bool notran = (data_.options.Trans == SLUD::NOTRANS);
612 magnitude_type anorm = function_map::plangs((notran ? (
char *)
"1" : (
char *)
"I"), &(data_.A), &(data_.grid));
616 #ifdef HAVE_AMESOS2_TIMERS
617 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
619 function_map::gstrf(&(data_.options), this->globalNumRows_,
620 this->globalNumCols_, anorm, &(data_.lu),
621 &(data_.grid), &(data_.stat), &info);
625 TEUCHOS_TEST_FOR_EXCEPTION( info > 0,
627 "L and U factors have been computed but U("
628 << info <<
"," << info <<
") is exactly zero "
629 "(i.e. U is singular)");
633 SLUD::Destroy_CompCol_Matrix_dist(&GA);
639 data_.options.Fact = SLUD::FACTORED;
640 same_symbolic_ =
true;
646 template <
class Matrix,
class Vector>
655 const size_t local_len_rhs = superlu_rowmap_->getLocalNumElements();
656 const global_size_type nrhs = X->getGlobalNumVectors();
657 const global_ordinal_type first_global_row_b = superlu_rowmap_->getMinGlobalIndex();
660 bvals_.resize(nrhs * local_len_rhs);
661 xvals_.resize(nrhs * local_len_rhs);
667 #ifdef HAVE_AMESOS2_TIMERS
668 Teuchos::TimeMonitor convTimer(this->timers_.vecConvTime_);
677 #ifdef HAVE_AMESOS2_TIMERS
678 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
683 copy_helper::do_get(B,
686 Teuchos::ptrInArg(*superlu_rowmap_));
707 if( !same_solve_struct_ ){
708 if( data_.options.SolveInitialized == SLUD::YES ){
709 function_map::SolveFinalize(&(data_.options), &(data_.solve_struct));
711 function_map::SolveInit(&(data_.options), &(data_.A), data_.perm_r.getRawPtr(),
712 data_.perm_c.getRawPtr(), as<SLUD::int_t>(nrhs), &(data_.lu),
713 &(data_.grid), &(data_.solve_struct));
717 same_solve_struct_ =
true;
721 if (data_.options.Equil == SLUD::YES && data_.rowequ) {
722 SLUD::int_t ld = as<SLUD::int_t>(local_len_rhs);
723 SLUD::slu_dist_mult<slu_type, magnitude_type> mult_op;
724 for(global_size_type j = 0; j < nrhs; ++j) {
725 for(
size_t i = 0; i < local_len_rhs; ++i) {
726 bvals_[i + j*ld] = mult_op(bvals_[i + j*ld], data_.R[first_global_row_b + i]);
734 #ifdef HAVE_AMESOS2_TIMERS
735 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
738 #if (SUPERLU_DIST_MAJOR_VERSION > 7)
739 function_map::gstrs(&(data_.options), as<SLUD::int_t>(this->globalNumRows_), &(data_.lu),
740 &(data_.scale_perm), &(data_.grid), bvals_.getRawPtr(),
741 as<SLUD::int_t>(local_len_rhs), as<SLUD::int_t>(first_global_row_b),
742 as<SLUD::int_t>(local_len_rhs), as<int>(nrhs),
743 &(data_.solve_struct), &(data_.stat), &ierr);
745 function_map::gstrs(as<SLUD::int_t>(this->globalNumRows_), &(data_.lu),
746 &(data_.scale_perm), &(data_.grid), bvals_.getRawPtr(),
747 as<SLUD::int_t>(local_len_rhs), as<SLUD::int_t>(first_global_row_b),
748 as<SLUD::int_t>(local_len_rhs), as<int>(nrhs),
749 &(data_.solve_struct), &(data_.stat), &ierr);
753 TEUCHOS_TEST_FOR_EXCEPTION( ierr < 0,
755 "Argument " << -ierr <<
" to gstrs had an illegal value" );
770 #ifdef HAVE_AMESOS2_TIMERS
771 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
773 SLUD::int_t ld = as<SLUD::int_t>(local_len_rhs);
774 function_map::permute_Dense_Matrix(as<SLUD::int_t>(first_global_row_b),
775 as<SLUD::int_t>(local_len_rhs),
776 data_.solve_struct.row_to_proc,
777 data_.solve_struct.inv_perm_c,
778 bvals_.getRawPtr(), ld,
779 xvals_.getRawPtr(), ld,
785 if (data_.options.Equil == SLUD::YES && data_.colequ) {
786 SLUD::int_t ld = as<SLUD::int_t>(local_len_rhs);
787 SLUD::slu_dist_mult<slu_type, magnitude_type> mult_op;
788 for(global_size_type j = 0; j < nrhs; ++j) {
789 for(
size_t i = 0; i < local_len_rhs; ++i) {
790 xvals_[i + j*ld] = mult_op(xvals_[i + j*ld], data_.C[first_global_row_b + i]);
798 #ifdef HAVE_AMESOS2_TIMERS
799 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
802 put_helper::do_put(X,
805 Teuchos::ptrInArg(*superlu_rowmap_));
812 template <
class Matrix,
class Vector>
817 return( this->globalNumRows_ == this->globalNumCols_ );
821 template <
class Matrix,
class Vector>
827 using Teuchos::getIntegralValue;
828 using Teuchos::ParameterEntryValidator;
830 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
832 if( parameterList->isParameter(
"npcol") || parameterList->isParameter(
"nprow") ){
833 TEUCHOS_TEST_FOR_EXCEPTION( !(parameterList->isParameter(
"nprow") &&
834 parameterList->isParameter(
"npcol")),
835 std::invalid_argument,
836 "nprow and npcol must be set together" );
838 SLUD::int_t nprow = parameterList->template get<SLUD::int_t>(
"nprow");
839 SLUD::int_t npcol = parameterList->template get<SLUD::int_t>(
"npcol");
841 TEUCHOS_TEST_FOR_EXCEPTION( nprow * npcol > this->getComm()->getSize(),
842 std::invalid_argument,
843 "nprow and npcol combination invalid" );
845 if( (npcol != data_.grid.npcol) || (nprow != data_.grid.nprow) ){
847 SLUD::superlu_gridexit(&(data_.grid));
849 SLUD::superlu_gridinit(data_.mat_comm, nprow, npcol, &(data_.grid));
853 TEUCHOS_TEST_FOR_EXCEPTION( this->control_.useTranspose_,
854 std::invalid_argument,
855 "SuperLU_DIST does not support solving the tranpose system" );
857 data_.options.Trans = SLUD::NOTRANS;
860 bool equil = parameterList->get<
bool>(
"Equil",
false);
861 data_.options.Equil = equil ? SLUD::YES : SLUD::NO;
863 if( parameterList->isParameter(
"RowPerm") ){
864 RCP<const ParameterEntryValidator> rowperm_validator = valid_params->getEntry(
"RowPerm").validator();
865 parameterList->getEntry(
"RowPerm").setValidator(rowperm_validator);
867 data_.options.RowPerm = getIntegralValue<SLUD::rowperm_t>(*parameterList,
"RowPerm");
870 if( parameterList->isParameter(
"LargeDiag_MC64-Options") ){
871 data_.largediag_mc64_job = parameterList->template get<int>(
"LargeDiag_MC64-Options");
874 if( parameterList->isParameter(
"ColPerm") ){
875 RCP<const ParameterEntryValidator> colperm_validator = valid_params->getEntry(
"ColPerm").validator();
876 parameterList->getEntry(
"ColPerm").setValidator(colperm_validator);
878 data_.options.ColPerm = getIntegralValue<SLUD::colperm_t>(*parameterList,
"ColPerm");
887 data_.options.IterRefine = SLUD::NOREFINE;
889 bool replace_tiny = parameterList->get<
bool>(
"ReplaceTinyPivot",
true);
890 data_.options.ReplaceTinyPivot = replace_tiny ? SLUD::YES : SLUD::NO;
892 if( parameterList->isParameter(
"IsContiguous") ){
893 is_contiguous_ = parameterList->get<
bool>(
"IsContiguous");
898 template <
class Matrix,
class Vector>
899 Teuchos::RCP<const Teuchos::ParameterList>
903 using Teuchos::tuple;
904 using Teuchos::ParameterList;
905 using Teuchos::EnhancedNumberValidator;
906 using Teuchos::setStringToIntegralParameter;
907 using Teuchos::setIntParameter;
908 using Teuchos::stringToIntegralParameterEntryValidator;
910 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
912 if( is_null(valid_params) ){
913 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
915 Teuchos::RCP<EnhancedNumberValidator<SLUD::int_t> > col_row_validator
916 = Teuchos::rcp(
new EnhancedNumberValidator<SLUD::int_t>() );
917 col_row_validator->setMin(1);
919 pl->set(
"npcol", data_.grid.npcol,
920 "Number of columns in the processor grid. "
921 "Must be set with nprow", col_row_validator);
922 pl->set(
"nprow", data_.grid.nprow,
923 "Number of rows in the SuperLU_DIST processor grid. "
924 "Must be set together with npcol", col_row_validator);
927 setStringToIntegralParameter<SLUD::trans_t>(
"Trans",
"NOTRANS",
928 "Solve for the transpose system or not",
929 tuple<string>(
"NOTRANS"),
930 tuple<string>(
"Do not solve with transpose"),
931 tuple<SLUD::trans_t>(SLUD::NOTRANS),
935 pl->set(
"Equil",
false,
"Whether to equilibrate the system before solve");
948 pl->set(
"ReplaceTinyPivot",
true,
949 "Specifies whether to replace tiny diagonals during LU factorization");
952 setStringToIntegralParameter<SLUD::rowperm_t>(
"RowPerm",
"NOROWPERM",
953 "Specifies how to permute the rows of the "
954 "matrix for sparsity preservation",
955 tuple<string>(
"NOROWPERM",
"LargeDiag_MC64"),
956 tuple<string>(
"Natural ordering",
957 "Duff/Koster algorithm"),
958 tuple<SLUD::rowperm_t>(SLUD::NOROWPERM,
959 SLUD::LargeDiag_MC64),
962 setIntParameter(
"LargeDiag_MC64-Options", 4,
"Options for RowPerm-LargeDiag_MC64", pl.getRawPtr());
965 setStringToIntegralParameter<SLUD::colperm_t>(
"ColPerm",
"PARMETIS",
966 "Specifies how to permute the columns of the "
967 "matrix for sparsity preservation",
968 tuple<string>(
"NATURAL",
"PARMETIS"),
969 tuple<string>(
"Natural ordering",
970 "ParMETIS ordering on A^T + A"),
971 tuple<SLUD::colperm_t>(SLUD::NATURAL,
975 pl->set(
"IsContiguous",
true,
"Whether GIDs contiguous");
984 template <
class Matrix,
class Vector>
988 SLUD::int_t& npcol)
const {
989 TEUCHOS_TEST_FOR_EXCEPTION( nprocs < 1,
990 std::invalid_argument,
991 "Number of MPI processes must be at least 1" );
992 SLUD::int_t c, r = 1;
993 while( r*r <= nprocs ) r++;
996 while( (r--)*c != nprocs ){
1002 if( r > 1 || nprocs < 9){
1009 template <
class Matrix,
class Vector>
1013 using Teuchos::Array;
1014 using Teuchos::ArrayView;
1015 using Teuchos::ptrInArg;
1020 #ifdef HAVE_AMESOS2_TIMERS
1021 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
1025 if( data_.A.Store != NULL ){
1026 SLUD::Destroy_SuperMatrix_Store_dist( &(data_.A) );
1027 data_.A.Store = NULL;
1030 Teuchos::RCP<const MatrixAdapter<Matrix> > redist_mat
1031 = this->matrixA_->get(ptrInArg(*superlu_rowmap_));
1033 int_t l_nnz, l_rows, g_rows, g_cols, fst_global_row;
1034 l_nnz = as<int_t>(redist_mat->getLocalNNZ());
1035 l_rows = as<int_t>(redist_mat->getLocalNumRows());
1036 g_rows = as<int_t>(redist_mat->getGlobalNumRows());
1038 fst_global_row = as<int_t>(superlu_rowmap_->getMinGlobalIndex());
1040 Kokkos::resize(nzvals_view_, l_nnz);
1041 Kokkos::resize(colind_view_, l_nnz);
1042 Kokkos::resize(rowptr_view_, l_rows + 1);
1045 #ifdef HAVE_AMESOS2_TIMERS
1046 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
1050 host_value_type_array,host_ordinal_type_array, host_size_type_array >::do_get(
1052 nzvals_view_, colind_view_, rowptr_view_,
1054 ptrInArg(*superlu_rowmap_),
1059 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != l_nnz,
1061 "Did not get the expected number of non-zero vals");
1064 SLUD::Dtype_t dtype = type_map::dtype;
1067 function_map::create_CompRowLoc_Matrix(&(data_.A),
1069 l_nnz, l_rows, fst_global_row,
1070 nzvals_view_.data(),
1071 colind_view_.data(),
1072 rowptr_view_.data(),
1074 dtype, SLUD::SLU_GE);
1081 template<
class Matrix,
class Vector>
1087 #endif // AMESOS2_SUPERLUDIST_DEF_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Return a const parameter list of all of the valid parameters that this->setParameterList(...) will accept.
Definition: Amesos2_SolverCore_def.hpp:558
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers...
Definition: Amesos2_SolverCore_decl.hpp:105
void computeRowPermutationLargeDiagMC64(SLUD::SuperMatrix &GA)
Compute the row permutation for option LargeDiag-MC64.
Definition: Amesos2_Superludist_def.hpp:337
Amesos2 interface to the distributed memory version of SuperLU.
Definition: Amesos2_Superludist_decl.hpp:90
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
global_size_type globalNumCols_
Number of global columns in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:479
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList) override
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:526
Superludist(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_Superludist_def.hpp:69
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_Superludist_def.hpp:823
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:267
Utility functions for Amesos2.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const override
Returns a pointer to the Teuchos::Comm communicator with this operator.
Definition: Amesos2_SolverCore_decl.hpp:363
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_Superludist_def.hpp:900
Provides definition of SuperLU_DIST types as well as conversions and type traits. ...
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_Superludist_def.hpp:392
Similar to get_ccs_helper , but used to get a CRS representation of the given matrix.
Definition: Amesos2_Util.hpp:659
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_Superludist_def.hpp:814
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > superlu_rowmap_
Maps rows of the matrix to processors in the SuperLU_DIST processor grid.
Definition: Amesos2_Superludist_decl.hpp:343
~Superludist()
Destructor.
Definition: Amesos2_Superludist_def.hpp:242
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using SuperLU_DIST.
Definition: Amesos2_Superludist_def.hpp:448
void get_default_grid_size(int nprocs, SLUD::int_t &nprow, SLUD::int_t &npcol) const
Definition: Amesos2_Superludist_def.hpp:986
bool in_grid_
true if this processor is in SuperLU_DISTS's 2D process grid
Definition: Amesos2_Superludist_decl.hpp:335
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:373
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
SuperLU_DIST specific solve.
Definition: Amesos2_Superludist_def.hpp:648
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal solver structures.
Definition: Amesos2_Superludist_def.hpp:1011
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
int numericFactorization_impl()
SuperLU_DIST specific numeric factorization.
Definition: Amesos2_Superludist_def.hpp:497