45#ifndef __IFPACK2_FASTILU_BASE_DEF_HPP__
46#define __IFPACK2_FASTILU_BASE_DEF_HPP__
49#include <KokkosKernels_Utils.hpp>
50#include <Kokkos_Timer.hpp>
51#include <Teuchos_TimeMonitor.hpp>
59template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 params_(Params::getDefaults()) {}
74template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
79 return mat_->getDomainMap();
82template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
87 return mat_->getRangeMap();
90template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92apply (
const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
93 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
94 Teuchos::ETransp mode,
98 const std::string timerName (
"Ifpack2::FastILU::apply");
99 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
100 if (timer.is_null ()) {
101 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
103 Teuchos::TimeMonitor timeMon (*timer);
105 if(!isInitialized() || !isComputed())
107 throw std::runtime_error(std::string(
"Called ") + getName() +
"::apply() without first calling initialize() and/or compute().");
109 if(X.getNumVectors() != Y.getNumVectors())
111 throw std::invalid_argument(getName() +
"::apply: X and Y have different numbers of vectors (pass X and Y with exactly matching dimensions)");
113 if(X.getLocalLength() != Y.getLocalLength())
115 throw std::invalid_argument(getName() +
"::apply: X and Y have different lengths (pass X and Y with exactly matching dimensions)");
119 int nvecs = X.getNumVectors();
120 auto nrowsX = X.getLocalLength();
121 auto nrowsY = Y.getLocalLength();
124 auto x2d = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
125 auto y2d = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
129 applyLocalPrec(x1d, y1d);
134 auto x2d = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
135 auto y2d = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
136 for(
int i = 0; i < nvecs; i++)
138 auto xColView1d = Kokkos::subview(x2d, Kokkos::ALL(), i);
139 auto yColView1d = Kokkos::subview(y2d, Kokkos::ALL(), i);
143 applyLocalPrec(x1d, y1d);
148template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
153 params_ = Params(List, getName());
156template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
160 const std::string timerName (
"Ifpack2::FastILU::initialize");
161 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
162 if (timer.is_null ()) {
163 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
165 Teuchos::TimeMonitor timeMon (*timer);
169 throw std::runtime_error(std::string(
"Called ") + getName() +
"::initialize() but matrix was null (call setMatrix() with a non-null matrix first)");
171 Kokkos::Timer copyTimer;
172 CrsArrayReader<Scalar, ImplScalar, LocalOrdinal, GlobalOrdinal, Node>::getStructure(mat_.get(), localRowPtrsHost_, localRowPtrs_, localColInds_);
173 crsCopyTime_ = copyTimer.seconds();
175 if (params_.use_metis)
177 const std::string timerNameMetis (
"Ifpack2::FastILU::Metis");
178 Teuchos::RCP<Teuchos::Time> timerMetis = Teuchos::TimeMonitor::lookupCounter (timerNameMetis);
179 if (timerMetis.is_null ()) {
180 timerMetis = Teuchos::TimeMonitor::getNewCounter (timerNameMetis);
182 Teuchos::TimeMonitor timeMonMetis (*timerMetis);
183 #ifdef HAVE_IFPACK2_METIS
184 idx_t nrows = localRowPtrsHost_.size() - 1;
187 metis_perm_ = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_perm"), nrows);
188 metis_iperm_ = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_iperm"), nrows);
191 auto localColIndsHost_ = Kokkos::create_mirror_view(localColInds_);
192 Kokkos::deep_copy(localColIndsHost_, localColInds_);
195 idx_t nnz = localColIndsHost_.size();
196 MetisArrayHost metis_rowptr;
197 MetisArrayHost metis_colidx;
199 bool metis_symmetrize =
true;
200 if (metis_symmetrize) {
202 using OrdinalArrayMirror =
typename OrdinalArray::host_mirror_type;
203 KokkosKernels::Impl::symmetrize_graph_symbolic_hashmap<
204 OrdinalArrayHost, OrdinalArrayMirror, MetisArrayHost, MetisArrayHost, Kokkos::HostSpace::execution_space>
205 (nrows, localRowPtrsHost_, localColIndsHost_, metis_rowptr, metis_colidx);
208 idx_t old_nnz = nnz = 0;
209 for (idx_t i = 0; i < nrows; i++) {
210 for (LocalOrdinal k = old_nnz; k < metis_rowptr(i+1); k++) {
211 if (metis_colidx(k) != i) {
212 metis_colidx(nnz) = metis_colidx(k);
216 old_nnz = metis_rowptr(i+1);
217 metis_rowptr(i+1) = nnz;
221 metis_rowptr = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_rowptr"), nrows+1);
222 metis_colidx = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_colidx"), nnz);
225 for (idx_t i = 0; i < nrows; i++) {
226 for (LocalOrdinal k = localRowPtrsHost_(i); k < localRowPtrsHost_(i+1); k++) {
227 if (localColIndsHost_(k) != i) {
228 metis_colidx(nnz) = localColIndsHost_(k);
232 metis_rowptr(i+1) = nnz;
237 int info = METIS_NodeND(&nrows, metis_rowptr.data(), metis_colidx.data(),
238 NULL, NULL, metis_perm_.data(), metis_iperm_.data());
239 if (METIS_OK != info) {
240 throw std::runtime_error(std::string(
"METIS_NodeND returned info = " + info));
244 throw std::runtime_error(std::string(
"TPL METIS is not enabled"));
253template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
260template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
266 throw std::runtime_error(getName() +
": initialize() must be called before compute()");
269 const std::string timerName (
"Ifpack2::FastILU::compute");
270 Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
271 if (timer.is_null ()) {
272 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
274 Teuchos::TimeMonitor timeMon (*timer);
277 Kokkos::Timer copyTimer;
278 CrsArrayReader<Scalar, ImplScalar, LocalOrdinal, GlobalOrdinal, Node>::getValues(mat_.get(), localValues_, localRowPtrsHost_);
279 crsCopyTime_ += copyTimer.seconds();
281 computedFlag_ =
true;
285template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 return computedFlag_;
292template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
293Teuchos::RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
300template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
307template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
314template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
321template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
328template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
335template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
349template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
354 throw std::runtime_error(std::string(
"Preconditioner type Ifpack2::Details::") + getName() +
" doesn't support checkLocalILU().");
357template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
362 throw std::runtime_error(std::string(
"Preconditioner type Ifpack2::Details::") + getName() +
" doesn't support checkLocalIC().");
365template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
368 std::ostringstream os;
370 os <<
"\"Ifpack2::Details::" << getName() <<
"\": {";
371 os <<
"Initialized: " << (isInitialized() ?
"true" :
"false") <<
", ";
372 os <<
"Computed: " << (isComputed() ?
"true" :
"false") <<
", ";
373 os <<
"Sweeps: " << getSweeps() <<
", ";
374 os <<
"Triangular solve type: " << getSpTrsvType() <<
", ";
375 if (getSpTrsvType() ==
"Fast") {
376 os <<
"# of triangular solve iterations: " << getNTrisol() <<
", ";
380 os <<
"Matrix: null";
384 os <<
"Global matrix dimensions: [" << mat_->getGlobalNumRows() <<
", " << mat_->getGlobalNumCols() <<
"]";
385 os <<
", Global nnz: " << mat_->getGlobalNumEntries();
390template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
392setMatrix(
const Teuchos::RCP<const TRowMatrix>& A)
396 throw std::invalid_argument(std::string(
"Ifpack2::Details::") + getName() +
"::setMatrix() called with a null matrix. Pass a non-null matrix.");
399 if(mat_.get() != A.get())
403 computedFlag_ =
false;
407template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
408typename FastILU_Base<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Params
414 p.sptrsv_algo = FastILU::SpTRSV::Fast;
426template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
427FastILU_Base<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
428Params::Params(
const Teuchos::ParameterList& pL, std::string precType)
430 *
this = getDefaults();
431 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude;
436 #define TYPE_ERROR(name, correctTypeName) {throw std::invalid_argument(precType + "::setParameters(): parameter \"" + name + "\" has the wrong type (must be " + correctTypeName + ")");}
437 #define CHECK_VALUE(param, member, cond, msg) {if(cond) {throw std::invalid_argument(precType + "::setParameters(): parameter \"" + param + "\" has value " + std::to_string(member) + " but " + msg);}}
440 if(pL.isParameter(
"metis"))
442 if(pL.isType<
bool>(
"metis"))
443 use_metis = pL.get<
bool>(
"metis");
445 TYPE_ERROR(
"metis",
"bool");
448 if(pL.isParameter(
"sweeps"))
450 if(pL.isType<
int>(
"sweeps"))
452 nFact = pL.get<
int>(
"sweeps");
453 CHECK_VALUE(
"sweeps", nFact, nFact < 1,
"must have a value of at least 1");
456 TYPE_ERROR(
"sweeps",
"int");
458 std::string sptrsv_type =
"Fast";
459 if(pL.isParameter(
"triangular solve type")) {
460 sptrsv_type = pL.get<std::string> (
"triangular solve type");
462 if (sptrsv_type ==
"Standard Host") {
463 sptrsv_algo = FastILU::SpTRSV::StandardHost;
464 }
else if (sptrsv_type ==
"Standard") {
465 sptrsv_algo = FastILU::SpTRSV::Standard;
469 if(pL.isParameter(
"triangular solve iterations"))
471 if(pL.isType<
int>(
"triangular solve iterations"))
473 nTrisol = pL.get<
int>(
"triangular solve iterations");
474 CHECK_VALUE(
"triangular solve iterations", nTrisol, nTrisol < 1,
"must have a value of at least 1");
477 TYPE_ERROR(
"triangular solve iterations",
"int");
480 if(pL.isParameter(
"level"))
482 if(pL.isType<
int>(
"level"))
484 level = pL.get<
int>(
"level");
486 else if(pL.isType<
double>(
"level"))
490 double dval = pL.get<
double>(
"level");
492 double fpart = modf(dval, &ipart);
494 CHECK_VALUE(
"level", level, fpart != 0,
"must be an integral value");
498 TYPE_ERROR(
"level",
"int");
500 CHECK_VALUE(
"level", level, level < 0,
"must be nonnegative");
503 if(pL.isParameter(
"damping factor"))
505 if(pL.isType<
double>(
"damping factor"))
506 omega = pL.get<
double>(
"damping factor");
507 else if(pL.isType<magnitude>(
"damping factor"))
508 omega = pL.get<magnitude>(
"damping factor");
510 TYPE_ERROR(
"damping factor",
"double or magnitude_type");
511 CHECK_VALUE(
"damping factor", omega, omega <= 0 || omega > 1,
"must be in the range (0, 1]");
514 if(pL.isParameter(
"shift"))
516 if(pL.isType<
double>(
"shift"))
517 shift = pL.get<
double>(
"shift");
518 else if(pL.isType<magnitude>(
"shift"))
519 shift = pL.get<magnitude>(
"shift");
521 TYPE_ERROR(
"shift",
"double or magnitude_type");
525 if(pL.isParameter(
"guess"))
527 if(pL.isType<
bool>(
"guess"))
528 guessFlag = pL.get<
bool>(
"guess");
530 TYPE_ERROR(
"guess",
"bool");
533 if(pL.isParameter(
"block size for ILU"))
535 if(pL.isType<
int>(
"block size for ILU"))
537 blockSizeILU = pL.get<
int>(
"block size for ILU");
538 CHECK_VALUE(
"block size for ILU", blockSizeILU, blockSizeILU < 1,
"must have a value of at least 1");
541 TYPE_ERROR(
"block size for ILU",
"int");
544 if(pL.isParameter(
"block size for SpTRSV"))
546 if(pL.isType<
int>(
"block size for SpTRSV"))
547 blockSize = pL.get<
int>(
"block size for SpTRSV");
549 TYPE_ERROR(
"block size for SpTRSV",
"int");
555#define IFPACK2_DETAILS_FASTILU_BASE_INSTANT(S, L, G, N) \
556template class Ifpack2::Details::FastILU_Base<S, L, G, N>;
Provides functions for retrieving local CRS arrays (row pointers, column indices, and values) from Tp...
The base class of the Ifpack2 FastILU wrappers (Filu, Fildl and Fic)
Definition Ifpack2_Details_FastILU_Base_decl.hpp:74
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >::impl_scalar_type ImplScalar
Kokkos scalar type.
Definition Ifpack2_Details_FastILU_Base_decl.hpp:81
double getInitializeTime() const
Get the time spent in the last initialize() call.
Definition Ifpack2_Details_FastILU_Base_def.hpp:323
virtual void checkLocalILU() const
Verify and print debug information about the underlying ILU preconditioner (only supported if this is...
Definition Ifpack2_Details_FastILU_Base_def.hpp:351
void setMatrix(const Teuchos::RCP< const TRowMatrix > &A)
Definition Ifpack2_Details_FastILU_Base_def.hpp:392
void compute()
Compute the preconditioner.
Definition Ifpack2_Details_FastILU_Base_def.hpp:262
double getComputeTime() const
Get the time spent in the last compute() call.
Definition Ifpack2_Details_FastILU_Base_def.hpp:330
double getCopyTime() const
Get the time spent deep copying local 3-array CRS out of the matrix.
Definition Ifpack2_Details_FastILU_Base_def.hpp:344
int getNumApply() const
Get the number of times apply() was called.
Definition Ifpack2_Details_FastILU_Base_def.hpp:316
void initialize()
Initialize the preconditioner.
Definition Ifpack2_Details_FastILU_Base_def.hpp:158
bool isInitialized() const
Whether initialize() has been called since the last time the matrix's structure was changed.
Definition Ifpack2_Details_FastILU_Base_def.hpp:255
FastILU_Base(Teuchos::RCP< const TRowMatrix > mat_)
Constructor.
Definition Ifpack2_Details_FastILU_Base_def.hpp:61
void setParameters(const Teuchos::ParameterList &List)
Validate parameters, and set defaults when parameters are not provided.
Definition Ifpack2_Details_FastILU_Base_def.hpp:150
void apply(const TMultiVec &X, TMultiVec &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Apply the preconditioner.
Definition Ifpack2_Details_FastILU_Base_def.hpp:92
double getApplyTime() const
Get the time spent in the last apply() call.
Definition Ifpack2_Details_FastILU_Base_def.hpp:337
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Get the domain map of the matrix.
Definition Ifpack2_Details_FastILU_Base_def.hpp:77
virtual void checkLocalIC() const
Verify and print debug information about the underlying IC preconditioner.
Definition Ifpack2_Details_FastILU_Base_def.hpp:359
int getNumInitialize() const
Get the number of times initialize() was called.
Definition Ifpack2_Details_FastILU_Base_def.hpp:302
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Get the range map of the matrix.
Definition Ifpack2_Details_FastILU_Base_def.hpp:85
Teuchos::RCP< const TRowMatrix > getMatrix() const
Get the current matrix.
Definition Ifpack2_Details_FastILU_Base_def.hpp:295
Kokkos::View< LocalOrdinal *, execution_space >::HostMirror OrdinalArrayHost
Array of LocalOrdinal on host.
Definition Ifpack2_Details_FastILU_Base_decl.hpp:93
std::string description() const
Return a brief description of the preconditioner, in YAML format.
Definition Ifpack2_Details_FastILU_Base_def.hpp:366
Kokkos::View< ImplScalar *, execution_space > ImplScalarArray
Array of Scalar on device.
Definition Ifpack2_Details_FastILU_Base_decl.hpp:95
int getNumCompute() const
Get the number of times compute() was called.
Definition Ifpack2_Details_FastILU_Base_def.hpp:309
bool isComputed() const
Whether compute() has been called since the last time the matrix's values or structure were changed.
Definition Ifpack2_Details_FastILU_Base_def.hpp:287
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74