45#ifndef AMESOS2_STRUMPACK_DEF_HPP
46#define AMESOS2_STRUMPACK_DEF_HPP
48#include <Teuchos_Tuple.hpp>
49#include <Teuchos_StandardParameterEntryValidators.hpp>
56#include <Teuchos_DefaultMpiComm.hpp>
57#include "StrumpackSparseSolverMPIDist.hpp"
59#include "StrumpackSparseSolver.hpp"
65 template <
class Matrix,
class Vector>
67 Teuchos::RCP<Vector> X,
68 Teuchos::RCP<const Vector> B)
74 using Teuchos::MpiComm;
76 using Teuchos::ParameterList;
77 using Teuchos::parameterList;
80 using Teuchos::rcp_dynamic_cast;
81 typedef global_ordinal_type GO;
83 typedef Tpetra::Map<local_ordinal_type, GO, node_type> map_type;
85 RCP<const Comm<int> > comm = this->
getComm ();
88 RCP<const MpiComm<int> > mpiComm =
89 rcp_dynamic_cast<const MpiComm<int> > (comm);
90 TEUCHOS_TEST_FOR_EXCEPTION
91 (mpiComm.is_null (), std::logic_error,
"Amesos2::STRUMPACK "
92 "constructor: The matrix's communicator is not an MpiComm!");
93 MPI_Comm rawMpiComm = (* (mpiComm->getRawMpiComm ())) ();
95 sp_ = Teuchos::RCP<strumpack::StrumpackSparseSolverMPIDist<scalar_type,GO>>
96 (
new strumpack::StrumpackSparseSolverMPIDist<scalar_type,GO>(rawMpiComm, this->
control_.verbose_));
98 sp_ = Teuchos::RCP<strumpack::StrumpackSparseSolver<scalar_type,GO>>
99 (
new strumpack::StrumpackSparseSolver<scalar_type,GO>(this->
control_.verbose_, this->root_));
107 RCP<ParameterList> default_params =
112 const size_t myNumRows = this->
matrixA_->getLocalNumRows();
113 const GO indexBase = this->
matrixA_->getRowMap ()->getIndexBase ();
115 rcp (
new map_type (this->
globalNumRows_, myNumRows, indexBase, comm));
124 template <
class Matrix,
class Vector>
133 template<
class Matrix,
class Vector>
137#ifdef HAVE_AMESOS2_TIMERS
138 Teuchos::TimeMonitor preOrderTime( this->timers_.preOrderTime_ );
153 template <
class Matrix,
class Vector>
157#ifdef HAVE_AMESOS2_TIMERS
158 Teuchos::TimeMonitor symFactTime( this->timers_.symFactTime_ );
161 strumpack::ReturnCode ret = sp_->reorder();
162 TEUCHOS_TEST_FOR_EXCEPTION( ret == strumpack::ReturnCode::MATRIX_NOT_SET,
164 "STRUMPACK matrix reordering failed: "
165 "matrix was not set." );
166 TEUCHOS_TEST_FOR_EXCEPTION( ret == strumpack::ReturnCode::REORDERING_ERROR,
168 "STRUMPACK matrix reordering failed." );
178 template <
class Matrix,
class Vector>
182#ifdef HAVE_AMESOS2_TIMERS
183 Teuchos::TimeMonitor numFactTimer(this->timers_.numFactTime_);
186 strumpack::ReturnCode ret = sp_->factor();
188 TEUCHOS_TEST_FOR_EXCEPTION( ret != strumpack::ReturnCode::SUCCESS,
190 "Error in STRUMPACK factorization." );
200 template <
class Matrix,
class Vector>
209 const size_t local_len_rhs = strumpack_rowmap_->getLocalNumElements();
210 const global_size_type nrhs = X->getGlobalNumVectors();
213 bvals_.resize(nrhs * local_len_rhs);
214 xvals_.resize(nrhs * local_len_rhs);
218#ifdef HAVE_AMESOS2_TIMERS
219 Teuchos::TimeMonitor convTimer(this->timers_.vecConvTime_);
223#ifdef HAVE_AMESOS2_TIMERS
224 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
228 copy_helper::do_get(B,
231 Teuchos::ptrInArg(*strumpack_rowmap_));
236#ifdef HAVE_AMESOS2_TIMERS
237 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
239 strumpack::DenseMatrixWrapper<scalar_type>
240 Bsp(local_len_rhs, nrhs, bvals_().getRawPtr(), local_len_rhs),
241 Xsp(local_len_rhs, nrhs, xvals_().getRawPtr(), local_len_rhs);
242 strumpack::ReturnCode ret =sp_->solve(Bsp, Xsp);
244 TEUCHOS_TEST_FOR_EXCEPTION( ret != strumpack::ReturnCode::SUCCESS,
246 "Error in STRUMPACK solve" );
251#ifdef HAVE_AMESOS2_TIMERS
252 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
257 put_helper::do_put(X,
260 Teuchos::ptrInArg(*strumpack_rowmap_));
264 const global_size_type ld_rhs = this->root_ ? X->getGlobalLength() : 0;
265 const size_t nrhs = X->getGlobalNumVectors();
266 bvals_.resize(nrhs * ld_rhs);
267 xvals_.resize(nrhs * ld_rhs);
270#ifdef HAVE_AMESOS2_TIMERS
271 Teuchos::TimeMonitor solveTimer(this->timers_.solveTime_);
274 strumpack::DenseMatrixWrapper<scalar_type>
275 Bsp(ld_rhs, nrhs, bvals_().getRawPtr(), ld_rhs),
276 Xsp(ld_rhs, nrhs, xvals_().getRawPtr(), ld_rhs);
277 strumpack::ReturnCode ret =sp_->solve(Bsp, Xsp);
279 TEUCHOS_TEST_FOR_EXCEPTION( ret != strumpack::ReturnCode::SUCCESS,
281 "Error in STRUMPACK solve" );
285#ifdef HAVE_AMESOS2_TIMERS
286 Teuchos::TimeMonitor redistTimer(this->timers_.vecRedistTime_);
292 ROOTED, this->rowIndexBase_);
299 template <
class Matrix,
class Vector>
305 return( this->globalNumRows_ == this->globalNumCols_ );
307 return( this->matrixA_->getGlobalNumRows() == this->matrixA_->getGlobalNumCols() );
317 template <
class Matrix,
class Vector>
323 using Teuchos::getIntegralValue;
324 using Teuchos::ParameterEntryValidator;
326 RCP<const Teuchos::ParameterList> valid_params = getValidParameters_impl();
328 if( parameterList->isParameter(
"Matching") ){
329 RCP<const ParameterEntryValidator> matching_validator = valid_params->getEntry(
"Matching").validator();
330 parameterList->getEntry(
"Matching").setValidator(matching_validator);
332 sp_->options().set_matching(getIntegralValue<strumpack::MatchingJob>(*parameterList,
"Matching"));
335 if( parameterList->isParameter(
"Ordering") ){
336 RCP<const ParameterEntryValidator> reordering_validator = valid_params->getEntry(
"Ordering").validator();
337 parameterList->getEntry(
"Ordering").setValidator(reordering_validator);
339 sp_->options().set_reordering_method(getIntegralValue<strumpack::ReorderingStrategy>(*parameterList,
"Ordering"));
342 if( parameterList->isParameter(
"ReplaceTinyPivot") ){
343 RCP<const ParameterEntryValidator> replacepivot_validator = valid_params->getEntry(
"ReplaceTinyPivot").validator();
344 parameterList->getEntry(
"ReplaceTinyPivot").setValidator(replacepivot_validator);
346 if( replacepivot_validator) {
347 sp_->options().enable_replace_tiny_pivots();
350 sp_->options().disable_replace_tiny_pivots();
354 if( parameterList->isParameter(
"IterRefine") ){
355 RCP<const ParameterEntryValidator> iter_refine_validator = valid_params->getEntry(
"IterRefine").validator();
356 parameterList->getEntry(
"IterRefine").setValidator(iter_refine_validator);
358 sp_->options().set_Krylov_solver(getIntegralValue<strumpack::KrylovSolver>(*parameterList,
"IterRefine"));
361 if( parameterList->isParameter(
"Compression") ){
362 RCP<const ParameterEntryValidator> compression_validator = valid_params->getEntry(
"Compression").validator();
363 parameterList->getEntry(
"Compression").setValidator(compression_validator);
365 sp_->options().set_compression(getIntegralValue<strumpack::CompressionType>(*parameterList,
"Compression"));
368 TEUCHOS_TEST_FOR_EXCEPTION( this->control_.useTranspose_,
369 std::invalid_argument,
370 "STRUMPACK does not support solving the tranpose system" );
374 template <
class Matrix,
class Vector>
375 Teuchos::RCP<const Teuchos::ParameterList>
379 using Teuchos::tuple;
380 using Teuchos::ParameterList;
381 using Teuchos::EnhancedNumberValidator;
382 using Teuchos::setStringToIntegralParameter;
383 using Teuchos::stringToIntegralParameterEntryValidator;
385 static Teuchos::RCP<const Teuchos::ParameterList> valid_params;
387 if( is_null(valid_params) ){
388 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
390 setStringToIntegralParameter<strumpack::MatchingJob>(
"Matching",
"NONE",
391 "Specifies how to permute the "
392 "matrix for numerical stability",
393 tuple<string>(
"NONE",
"MAX_CARDINALITY",
"MAX_SMALLEST_DIAGONAL",
"MAX_SMALLEST_DIAGONAL_2",
"MAX_DIAGONAL_SUM",
"MAX_DIAGONAL_PRODUCT_SCALING",
"COMBBLAS"),
394 tuple<string>(
"NONE",
"MAX_CARDINALITY",
"MAX_SMALLEST_DIAGONAL",
"MAX_SMALLEST_DIAGONAL_2",
"MAX_DIAGONAL_SUM",
"MAX_DIAGONAL_PRODUCT_SCALING",
"COMBBLAS"),
395 tuple<strumpack::MatchingJob>(strumpack::MatchingJob::NONE,
396 strumpack::MatchingJob::MAX_CARDINALITY,
397 strumpack::MatchingJob::MAX_SMALLEST_DIAGONAL,
398 strumpack::MatchingJob::MAX_SMALLEST_DIAGONAL_2,
399 strumpack::MatchingJob::MAX_DIAGONAL_SUM,
400 strumpack::MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING,
401 strumpack::MatchingJob::COMBBLAS),
404#if defined(STRUMPACK_USE_PARMETIS)
405 std::string default_ordering(
"PARMETIS");
407 std::string default_ordering(
"METIS");
409 setStringToIntegralParameter<strumpack::ReorderingStrategy>(
"Ordering", default_ordering,
410 "Specifies how to permute the "
411 "matrix for sparsity preservation",
412 tuple<string>(
"NATURAL",
"PARMETIS",
"METIS",
"SCOTCH",
"GEOMETRIC",
"PTSCOTCH",
"RCM"),
413 tuple<string>(
"Natural ordering",
417 "Geometric ordering",
420 tuple<strumpack::ReorderingStrategy>(strumpack::ReorderingStrategy::NATURAL,
421 strumpack::ReorderingStrategy::PARMETIS,
422 strumpack::ReorderingStrategy::METIS,
423 strumpack::ReorderingStrategy::SCOTCH,
424 strumpack::ReorderingStrategy::GEOMETRIC,
425 strumpack::ReorderingStrategy::PTSCOTCH,
426 strumpack::ReorderingStrategy::RCM),
429 pl->set(
"ReplaceTinyPivot",
true,
"Specifies whether to replace tiny diagonals during LU factorization");
434 setStringToIntegralParameter<strumpack::KrylovSolver>(
"IterRefine",
"DIRECT",
435 "Type of iterative refinement to use",
436 tuple<string>(
"AUTO",
"DIRECT",
"REFINE",
"PREC_GMRES",
"GMRES",
"PREC_BICGSTAB",
"BICGSTAB"),
437 tuple<string>(
"Use iterative refinement if no compression is used, otherwise use GMRes.",
438 "Single application of the multifrontal solver.",
439 "Iterative refinement.",
440 "Preconditioned GMRes.",
441 "UN-preconditioned GMRes.",
442 "Preconditioned BiCGStab.",
443 "UN-preconditioned BiCGStab."),
444 tuple<strumpack::KrylovSolver>(strumpack::KrylovSolver::AUTO,
445 strumpack::KrylovSolver::DIRECT,
446 strumpack::KrylovSolver::REFINE,
447 strumpack::KrylovSolver::PREC_GMRES,
448 strumpack::KrylovSolver::GMRES,
449 strumpack::KrylovSolver::PREC_BICGSTAB,
450 strumpack::KrylovSolver::BICGSTAB),
455 setStringToIntegralParameter<strumpack::CompressionType>(
"Compression",
"NONE",
456 "Type of compression to use",
457 tuple<string>(
"NONE",
"HSS",
"BLR",
"HODLR",
"LOSSLESS",
"LOSSY"),
458 tuple<string>(
"No compression, purely direct solver.",
459 "HSS compression of frontal matrices.",
460 "Block low-rank compression of fronts.",
461 "Hierarchically Off-diagonal Low-Rank compression of frontal matrices.",
462 "Lossless compresssion.",
463 "Lossy compresssion."),
464 tuple<strumpack::CompressionType>(strumpack::CompressionType::NONE,
465 strumpack::CompressionType::HSS,
466 strumpack::CompressionType::BLR,
467 strumpack::CompressionType::HODLR,
468 strumpack::CompressionType::LOSSLESS,
469 strumpack::CompressionType::LOSSY),
487 template <
class Matrix,
class Vector>
490 using Teuchos::Array;
491 using Teuchos::ArrayView;
492 using Teuchos::ptrInArg;
494 using Teuchos::rcp_dynamic_cast;
498 using Teuchos::MpiComm;
500 #ifdef HAVE_AMESOS2_TIMERS
501 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
504 Teuchos::RCP<const MatrixAdapter<Matrix> > redist_mat
505 = this->matrixA_->get(ptrInArg(*strumpack_rowmap_));
507 typedef global_ordinal_type GO;
509 l_nnz = as<GO>(redist_mat->getLocalNNZ());
510 l_rows = as<GO>(redist_mat->getLocalNumRows());
512 RCP<const Comm<int> > comm = this->getComm ();
513 RCP<const MpiComm<int> > mpiComm =
514 rcp_dynamic_cast<const MpiComm<int> > (comm);
516 const int numProcs = comm->getSize ();
517 const int myRank = comm->getRank ();
518 Array<GO> dist(numProcs+1);
520 dist[myRank+1] = as<GO>(strumpack_rowmap_->getMaxGlobalIndex()) + 1;
522 MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, dist.data()+1,
sizeof(GO), MPI_BYTE,
523 mpiComm->getRawMpiComm()->operator()());
524 nzvals_.resize(l_nnz);
525 colind_.resize(l_nnz);
526 rowptr_.resize(l_rows + 1);
531#ifdef HAVE_AMESOS2_TIMERS
532 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
535 Util::get_crs_helper<
537 scalar_type, GO, GO >::do_get(redist_mat.ptr(),
538 nzvals_(), colind_(), rowptr_(),
540 ptrInArg(*strumpack_rowmap_),
546 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != l_nnz,
548 "Did not get the expected number of non-zero vals");
551 sp_->set_distributed_csr_matrix
552 (l_rows, rowptr_.getRawPtr(), colind_.getRawPtr(),
553 nzvals_.getRawPtr(), dist.getRawPtr(),
false);
555#ifdef HAVE_AMESOS2_TIMERS
556 Teuchos::TimeMonitor convTimer(this->timers_.mtxConvTime_);
559 typedef global_ordinal_type GO;
563 nzvals_.resize(this->globalNumNonZeros_);
564 colind_.resize(this->globalNumNonZeros_);
565 rowptr_.resize(this->globalNumRows_ + 1);
568#ifdef HAVE_AMESOS2_TIMERS
569 Teuchos::TimeMonitor mtxRedistTimer( this->timers_.mtxRedistTime_ );
572 Util::get_crs_helper<
574 scalar_type, GO, GO >::do_get(this->matrixA_.ptr(),
575 nzvals_(), colind_(), rowptr_(),
581 TEUCHOS_TEST_FOR_EXCEPTION( nnz_ret != this->globalNumNonZeros_,
583 "Did not get the expected number of non-zero vals");
586 sp_->set_csr_matrix(this->globalNumRows_, rowptr_.getRawPtr(), colind_.getRawPtr(),
587 nzvals_.getRawPtr(),
false);
594 template<
class Matrix,
class Vector>
@ ROOTED
Definition: Amesos2_TypeDecl.hpp:127
@ ARBITRARY
Definition: Amesos2_TypeDecl.hpp:143
Utility functions for Amesos2.
A Matrix adapter interface for Amesos2.
Definition: Amesos2_MatrixAdapter_decl.hpp:76
Amesos2 interface to STRUMPACK direct solver and preconditioner.
Definition: Amesos2_STRUMPACK_decl.hpp:72
void setParameters_impl(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Definition: Amesos2_STRUMPACK_def.hpp:319
~STRUMPACK()
Destructor.
Definition: Amesos2_STRUMPACK_def.hpp:125
bool loadA_impl(EPhase current_phase)
Reads matrix data into internal solver structures.
Definition: Amesos2_STRUMPACK_def.hpp:489
bool matrixShapeOK_impl() const
Determines whether the shape of the matrix is OK for this solver.
Definition: Amesos2_STRUMPACK_def.hpp:301
int solve_impl(const Teuchos::Ptr< MultiVecAdapter< Vector > > X, const Teuchos::Ptr< const MultiVecAdapter< Vector > > B) const
STRUMPACK specific solve.
Definition: Amesos2_STRUMPACK_def.hpp:202
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters_impl() const
Definition: Amesos2_STRUMPACK_def.hpp:376
int preOrdering_impl()
Performs pre-ordering on the matrix to increase efficiency.
Definition: Amesos2_STRUMPACK_def.hpp:135
int numericFactorization_impl()
STRUMPACK specific numeric factorization.
Definition: Amesos2_STRUMPACK_def.hpp:180
int symbolicFactorization_impl()
Perform symbolic factorization of the matrix using STRUMPACK.
Definition: Amesos2_STRUMPACK_def.hpp:155
STRUMPACK(Teuchos::RCP< const Matrix > A, Teuchos::RCP< Vector > X, Teuchos::RCP< const Vector > B)
Initialize from Teuchos::RCP.
Definition: Amesos2_STRUMPACK_def.hpp:66
Amesos2::SolverCore: A templated interface for interaction with third-party direct sparse solvers.
Definition: Amesos2_SolverCore_decl.hpp:106
Teuchos::RCP< const MatrixAdapter< Matrix > > matrixA_
The LHS operator.
Definition: Amesos2_SolverCore_decl.hpp:455
super_type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶meterList)
Set/update internal variables and solver options.
Definition: Amesos2_SolverCore_def.hpp:518
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a const parameter list of all of the valid parameters that this->setParameterList(....
Definition: Amesos2_SolverCore_def.hpp:550
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
global_size_type globalNumRows_
Number of global rows in matrixA_.
Definition: Amesos2_SolverCore_decl.hpp:476
Control control_
Parameters for solving.
Definition: Amesos2_SolverCore_decl.hpp:494
EPhase
Used to indicate a phase in the direct solution.
Definition: Amesos2_TypeDecl.hpp:65
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176
Helper class for getting 1-D copies of multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:267
Helper class for putting 1-D data arrays into multivectors.
Definition: Amesos2_MultiVecAdapter_decl.hpp:373