45 #ifndef __IFPACK2_FASTILU_BASE_DEF_HPP__
46 #define __IFPACK2_FASTILU_BASE_DEF_HPP__
49 #include "Tpetra_BlockCrsMatrix.hpp"
50 #include "Tpetra_BlockCrsMatrix_Helpers.hpp"
51 #include "Ifpack2_Details_getCrsMatrix.hpp"
52 #include <KokkosKernels_Utils.hpp>
53 #include <Kokkos_Timer.hpp>
62 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 params_(Params::getDefaults()) {}
77 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 return mat_->getDomainMap();
85 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 return mat_->getRangeMap();
93 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 apply (
const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
96 Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
101 const std::string timerName (
"Ifpack2::FastILU::apply");
108 if(!isInitialized() || !isComputed())
110 throw std::runtime_error(std::string(
"Called ") + getName() +
"::apply() without first calling initialize() and/or compute().");
112 if(X.getNumVectors() != Y.getNumVectors())
114 throw std::invalid_argument(getName() +
"::apply: X and Y have different numbers of vectors (pass X and Y with exactly matching dimensions)");
116 if(X.getLocalLength() != Y.getLocalLength())
118 throw std::invalid_argument(getName() +
"::apply: X and Y have different lengths (pass X and Y with exactly matching dimensions)");
122 int nvecs = X.getNumVectors();
123 auto nrowsX = X.getLocalLength();
124 auto nrowsY = Y.getLocalLength();
127 auto x2d = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
128 auto y2d = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
132 applyLocalPrec(x1d, y1d);
137 auto x2d = X.getLocalViewDevice(Tpetra::Access::ReadOnly);
138 auto y2d = Y.getLocalViewDevice(Tpetra::Access::ReadWrite);
139 for(
int i = 0; i < nvecs; i++)
141 auto xColView1d = Kokkos::subview(x2d, Kokkos::ALL(), i);
142 auto yColView1d = Kokkos::subview(y2d, Kokkos::ALL(), i);
143 ImplScalarArray x1d (const_cast<ImplScalar*>(xColView1d.data()), nrowsX);
144 ImplScalarArray y1d (const_cast<ImplScalar*>(yColView1d.data()), nrowsY);
146 applyLocalPrec(x1d, y1d);
151 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
156 params_ = Params(List, getName());
159 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
163 return params_.blockCrs && params_.blockCrsSize > 1;
166 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
170 const std::string timerName (
"Ifpack2::FastILU::initialize");
179 throw std::runtime_error(std::string(
"Called ") + getName() +
"::initialize() but matrix was null (call setMatrix() with a non-null matrix first)");
183 auto crs_matrix = Ifpack2::Details::getCrsMatrix(this->mat_);
185 if (params_.fillBlocks) {
187 auto crs_matrix_block_filled = Tpetra::fillLogicalBlocks(*crs_matrix, params_.blockCrsSize);
188 auto bcrs_matrix = Tpetra::convertToBlockCrsMatrix(*crs_matrix_block_filled, params_.blockCrsSize);
193 auto bcrs_matrix = Tpetra::convertToBlockCrsMatrix(*crs_matrix, params_.blockCrsSize);
198 Kokkos::Timer copyTimer;
199 CrsArrayReader<Scalar, ImplScalar, LocalOrdinal, GlobalOrdinal, Node>::getStructure(mat_.get(), localRowPtrsHost_, localRowPtrs_, localColInds_);
200 CrsArrayReader<Scalar, ImplScalar, LocalOrdinal, GlobalOrdinal, Node>::getValues(mat_.get(), localValues_, localRowPtrsHost_);
201 crsCopyTime_ = copyTimer.seconds();
203 if (params_.use_metis)
205 assert(!params_.blockCrs);
206 const std::string timerNameMetis (
"Ifpack2::FastILU::Metis");
212 #ifdef HAVE_IFPACK2_METIS
213 idx_t nrows = localRowPtrsHost_.size() - 1;
216 metis_perm_ = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_perm"), nrows);
217 metis_iperm_ = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_iperm"), nrows);
220 auto localColIndsHost_ = Kokkos::create_mirror_view(localColInds_);
221 Kokkos::deep_copy(localColIndsHost_, localColInds_);
224 idx_t nnz = localColIndsHost_.size();
225 MetisArrayHost metis_rowptr;
226 MetisArrayHost metis_colidx;
228 bool metis_symmetrize =
true;
229 if (metis_symmetrize) {
231 using OrdinalArrayMirror =
typename OrdinalArray::host_mirror_type;
232 KokkosKernels::Impl::symmetrize_graph_symbolic_hashmap<
233 OrdinalArrayHost, OrdinalArrayMirror, MetisArrayHost, MetisArrayHost, Kokkos::HostSpace::execution_space>
234 (nrows, localRowPtrsHost_, localColIndsHost_, metis_rowptr, metis_colidx);
237 idx_t old_nnz = nnz = 0;
238 for (idx_t i = 0; i < nrows; i++) {
239 for (LocalOrdinal k = old_nnz; k < metis_rowptr(i+1); k++) {
240 if (metis_colidx(k) != i) {
241 metis_colidx(nnz) = metis_colidx(k);
245 old_nnz = metis_rowptr(i+1);
246 metis_rowptr(i+1) = nnz;
250 metis_rowptr = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_rowptr"), nrows+1);
251 metis_colidx = MetisArrayHost(Kokkos::ViewAllocateWithoutInitializing(
"metis_colidx"), nnz);
254 for (idx_t i = 0; i < nrows; i++) {
255 for (LocalOrdinal k = localRowPtrsHost_(i); k < localRowPtrsHost_(i+1); k++) {
256 if (localColIndsHost_(k) != i) {
257 metis_colidx(nnz) = localColIndsHost_(k);
261 metis_rowptr(i+1) = nnz;
266 int info = METIS_NodeND(&nrows, metis_rowptr.data(), metis_colidx.data(),
267 NULL, NULL, metis_perm_.data(), metis_iperm_.data());
268 if (METIS_OK != info) {
269 throw std::runtime_error(std::string(
"METIS_NodeND returned info = " + info));
273 throw std::runtime_error(std::string(
"TPL METIS is not enabled"));
282 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
295 throw std::runtime_error(getName() +
": initialize() must be called before compute()");
298 const std::string timerName (
"Ifpack2::FastILU::compute");
306 Kokkos::Timer copyTimer;
307 CrsArrayReader<Scalar, ImplScalar, LocalOrdinal, GlobalOrdinal, Node>::getValues(mat_.get(), localValues_, localRowPtrsHost_);
308 crsCopyTime_ += copyTimer.seconds();
310 computedFlag_ =
true;
314 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
318 return computedFlag_;
321 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
329 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
336 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
343 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
350 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
357 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
364 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
371 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
378 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
383 throw std::runtime_error(std::string(
"Preconditioner type Ifpack2::Details::") + getName() +
" doesn't support checkLocalILU().");
386 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
391 throw std::runtime_error(std::string(
"Preconditioner type Ifpack2::Details::") + getName() +
" doesn't support checkLocalIC().");
394 template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
397 std::ostringstream os;
399 os <<
"\"Ifpack2::Details::" << getName() <<
"\": {";
400 os <<
"Initialized: " << (isInitialized() ?
"true" :
"false") <<
", ";
401 os <<
"Computed: " << (isComputed() ?
"true" :
"false") <<
", ";
402 os <<
"Sweeps: " << getSweeps() <<
", ";
403 os <<
"Triangular solve type: " << getSpTrsvType() <<
", ";
404 if (getSpTrsvType() ==
"Fast") {
405 os <<
"# of triangular solve iterations: " << getNTrisol() <<
", ";
409 os <<
"Matrix: null";
413 os <<
"Global matrix dimensions: [" << mat_->getGlobalNumRows() <<
", " << mat_->getGlobalNumCols() <<
"]";
414 os <<
", Global nnz: " << mat_->getGlobalNumEntries();
419 template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
425 throw std::invalid_argument(std::string(
"Ifpack2::Details::") + getName() +
"::setMatrix() called with a null matrix. Pass a non-null matrix.");
428 if(mat_.get() != A.
get())
432 computedFlag_ =
false;
436 template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
443 p.sptrsv_algo = FastILU::SpTRSV::Fast;
454 p.fillBlocks =
false;
458 template<
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
typename Node>
459 FastILU_Base<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
462 *
this = getDefaults();
467 #define TYPE_ERROR(name, correctTypeName) {throw std::invalid_argument(precType + "::setParameters(): parameter \"" + name + "\" has the wrong type (must be " + correctTypeName + ")");}
468 #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);}}
473 if(pL.
isType<
bool>(
"metis"))
474 use_metis = pL.
get<
bool>(
"metis");
476 TYPE_ERROR(
"metis",
"bool");
481 if(pL.
isType<
int>(
"sweeps"))
483 nFact = pL.
get<
int>(
"sweeps");
484 CHECK_VALUE(
"sweeps", nFact, nFact < 1,
"must have a value of at least 1");
487 TYPE_ERROR(
"sweeps",
"int");
489 std::string sptrsv_type =
"Fast";
491 sptrsv_type = pL.
get<std::string> (
"triangular solve type");
493 if (sptrsv_type ==
"Standard Host") {
494 sptrsv_algo = FastILU::SpTRSV::StandardHost;
495 }
else if (sptrsv_type ==
"Standard") {
496 sptrsv_algo = FastILU::SpTRSV::Standard;
502 if(pL.
isType<
int>(
"triangular solve iterations"))
504 nTrisol = pL.
get<
int>(
"triangular solve iterations");
505 CHECK_VALUE(
"triangular solve iterations", nTrisol, nTrisol < 1,
"must have a value of at least 1");
508 TYPE_ERROR(
"triangular solve iterations",
"int");
513 if(pL.
isType<
int>(
"level"))
515 level = pL.
get<
int>(
"level");
517 else if(pL.
isType<
double>(
"level"))
521 double dval = pL.
get<
double>(
"level");
523 double fpart = modf(dval, &ipart);
525 CHECK_VALUE(
"level", level, fpart != 0,
"must be an integral value");
529 TYPE_ERROR(
"level",
"int");
531 CHECK_VALUE(
"level", level, level < 0,
"must be nonnegative");
535 if(pL.
isType<
double>(
"damping factor"))
536 omega = pL.
get<
double>(
"damping factor");
538 TYPE_ERROR(
"damping factor",
"double");
542 if(pL.
isType<
double>(
"shift"))
543 shift = pL.
get<
double>(
"shift");
545 TYPE_ERROR(
"shift",
"double");
550 if(pL.
isType<
bool>(
"guess"))
551 guessFlag = pL.
get<
bool>(
"guess");
553 TYPE_ERROR(
"guess",
"bool");
558 if(pL.
isType<
int>(
"block size for ILU"))
560 blockSizeILU = pL.
get<
int>(
"block size for ILU");
561 CHECK_VALUE(
"block size for ILU", blockSizeILU, blockSizeILU < 1,
"must have a value of at least 1");
564 TYPE_ERROR(
"block size for ILU",
"int");
569 if(pL.
isType<
int>(
"block size for SpTRSV"))
570 blockSize = pL.
get<
int>(
"block size for SpTRSV");
572 TYPE_ERROR(
"block size for SpTRSV",
"int");
577 if(pL.
isType<
bool>(
"block crs"))
578 blockCrs = pL.
get<
bool>(
"block crs");
580 TYPE_ERROR(
"block crs",
"bool");
585 if(pL.
isType<
int>(
"block crs block size"))
586 blockCrsSize = pL.
get<
int>(
"block crs block size");
588 TYPE_ERROR(
"block crs block size",
"int");
593 if(pL.
isType<
bool>(
"fill blocks for input"))
594 blockCrsSize = pL.
get<
bool>(
"fill blocks for input");
596 TYPE_ERROR(
"fill blocks for input",
"bool");
603 #define IFPACK2_DETAILS_FASTILU_BASE_INSTANT(S, L, G, N) \
604 template class Ifpack2::Details::FastILU_Base<S, L, G, N>;
int getNumCompute() const
Get the number of times compute() was called.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:338
virtual void checkLocalIC() const
Verify and print debug information about the underlying IC preconditioner.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:388
double getComputeTime() const
Get the time spent in the last compute() call.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:359
Kokkos::View< LocalOrdinal *, execution_space >::HostMirror OrdinalArrayHost
Array of LocalOrdinal on host.
Definition: Ifpack2_Details_FastILU_Base_decl.hpp:93
T & get(const std::string &name, T def_value)
double getCopyTime() const
Get the time spent deep copying local 3-array CRS out of the matrix.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:373
double getInitializeTime() const
Get the time spent in the last initialize() call.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:352
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const
Get the range map of the matrix.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:88
double getApplyTime() const
Get the time spent in the last apply() call.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:366
Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const
Get the domain map of the matrix.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:80
bool isParameter(const std::string &name) const
Teuchos::RCP< const TRowMatrix > getMatrix() const
Get the current matrix.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:324
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:380
int getNumApply() const
Get the number of times apply() was called.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:345
The base class of the Ifpack2 FastILU wrappers (Filu, Fildl and Fic)
Definition: Ifpack2_Details_FastILU_Base_decl.hpp:71
void setParameters(const Teuchos::ParameterList &List)
Validate parameters, and set defaults when parameters are not provided.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:153
std::string description() const
Return a brief description of the preconditioner, in YAML format.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:395
void compute()
Compute the preconditioner.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:291
Kokkos::View< ImplScalar *, execution_space > ImplScalarArray
Array of Scalar on device.
Definition: Ifpack2_Details_FastILU_Base_decl.hpp:95
FastILU_Base(Teuchos::RCP< const TRowMatrix > mat_)
Constructor.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:64
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:316
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:284
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:168
bool isType(const std::string &name) const
void setMatrix(const Teuchos::RCP< const TRowMatrix > &A)
Definition: Ifpack2_Details_FastILU_Base_def.hpp:421
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:95
int getNumInitialize() const
Get the number of times initialize() was called.
Definition: Ifpack2_Details_FastILU_Base_def.hpp:331
Provides functions for retrieving local CRS arrays (row pointers, column indices, and values) from Tp...