41 #ifndef IFPACK2_CRSRILUK_DEF_HPP
42 #define IFPACK2_CRSRILUK_DEF_HPP
44 #include "Ifpack2_LocalFilter.hpp"
45 #include "Tpetra_CrsMatrix.hpp"
46 #include "Ifpack2_LocalSparseTriangularSolver.hpp"
47 #include "Ifpack2_Details_getParamTryingTypes.hpp"
51 template<
class MatrixType>
57 isInitialized_ (false),
62 initializeTime_ (0.0),
73 template<
class MatrixType>
79 isInitialized_ (false),
84 initializeTime_ (0.0),
95 template<
class MatrixType>
98 template<
class MatrixType>
105 template<
class MatrixType>
114 isAllocated_ =
false;
115 isInitialized_ =
false;
117 A_local_ = Teuchos::null;
118 Graph_ = Teuchos::null;
125 if (! L_solver_.is_null ()) {
126 L_solver_->setMatrix (Teuchos::null);
128 if (! U_solver_.is_null ()) {
129 U_solver_->setMatrix (Teuchos::null);
141 template<
class MatrixType>
146 L_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor "
147 "is null. Please call initialize() (and preferably also compute()) "
148 "before calling this method. If the input matrix has not yet been set, "
149 "you must first call setMatrix() with a nonnull input matrix before you "
150 "may call initialize() or compute().");
155 template<
class MatrixType>
156 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
163 D_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor "
164 "(of diagonal entries) is null. Please call initialize() (and "
165 "preferably also compute()) before calling this method. If the input "
166 "matrix has not yet been set, you must first call setMatrix() with a "
167 "nonnull input matrix before you may call initialize() or compute().");
172 template<
class MatrixType>
177 U_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U factor "
178 "is null. Please call initialize() (and preferably also compute()) "
179 "before calling this method. If the input matrix has not yet been set, "
180 "you must first call setMatrix() with a nonnull input matrix before you "
181 "may call initialize() or compute().");
186 template<
class MatrixType>
189 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getNodeSmootherComplexity: "
190 "The input matrix A is null. Please call setMatrix() with a nonnull "
191 "input matrix, then call compute(), before calling this method.");
193 if(!L_.is_null() && !U_.is_null())
194 return A_->getNodeNumEntries() + L_->getNodeNumEntries() + U_->getNodeNumEntries();
201 template<
class MatrixType>
207 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
208 "The matrix is null. Please call setMatrix() with a nonnull input "
209 "before calling this method.");
213 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
214 "The computed graph is null. Please call initialize() before calling "
216 return Graph_->getL_Graph ()->getDomainMap ();
220 template<
class MatrixType>
226 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
227 "The matrix is null. Please call setMatrix() with a nonnull input "
228 "before calling this method.");
232 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
233 "The computed graph is null. Please call initialize() before calling "
235 return Graph_->getL_Graph ()->getRangeMap ();
239 template<
class MatrixType>
245 if (! isAllocated_) {
254 L_ =
rcp (
new crs_matrix_type (Graph_->getL_Graph ()));
255 U_ =
rcp (
new crs_matrix_type (Graph_->getU_Graph ()));
256 L_->setAllToScalar (STS::zero ());
257 U_->setAllToScalar (STS::zero ());
262 D_ =
rcp (
new vec_type (Graph_->getL_Graph ()->getRowMap ()));
268 template<
class MatrixType>
273 using Details::getParamTryingTypes;
274 const char prefix[] =
"Ifpack2::RILUK: ";
281 double overalloc = 2.;
292 const std::string paramName (
"fact: iluk level-of-fill");
293 getParamTryingTypes<int, int, global_ordinal_type, double, float>
294 (fillLevel, params, paramName, prefix);
299 const std::string paramName (
"fact: absolute threshold");
300 getParamTryingTypes<magnitude_type, magnitude_type, double>
301 (absThresh, params, paramName, prefix);
304 const std::string paramName (
"fact: relative threshold");
305 getParamTryingTypes<magnitude_type, magnitude_type, double>
306 (relThresh, params, paramName, prefix);
309 const std::string paramName (
"fact: relax value");
310 getParamTryingTypes<magnitude_type, magnitude_type, double>
311 (relaxValue, params, paramName, prefix);
314 const std::string paramName (
"fact: iluk overalloc");
315 getParamTryingTypes<double, double>
316 (overalloc, params, paramName, prefix);
320 L_solver_->setParameters(params);
321 U_solver_->setParameters(params);
327 LevelOfFill_ = fillLevel;
328 Overalloc_ = overalloc;
335 if (absThresh != -STM::one ()) {
336 Athresh_ = absThresh;
338 if (relThresh != -STM::one ()) {
339 Rthresh_ = relThresh;
341 if (relaxValue != -STM::one ()) {
342 RelaxValue_ = relaxValue;
347 template<
class MatrixType>
354 template<
class MatrixType>
361 template<
class MatrixType>
367 using Teuchos::rcp_dynamic_cast;
368 using Teuchos::rcp_implicit_cast;
373 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
374 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
381 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
382 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
383 if (! A_lf_r.is_null ()) {
384 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
390 return rcp (
new LocalFilter<row_matrix_type> (A));
395 template<
class MatrixType>
400 using Teuchos::rcp_const_cast;
401 using Teuchos::rcp_dynamic_cast;
402 using Teuchos::rcp_implicit_cast;
408 const char prefix[] =
"Ifpack2::RILUK::initialize: ";
411 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
412 "call setMatrix() with a nonnull input before calling this method.");
414 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
415 "fill complete. You may not invoke initialize() or compute() with this "
416 "matrix until the matrix is fill complete. If your matrix is a "
417 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
418 "range Maps, if appropriate) before calling this method.");
431 isInitialized_ =
false;
432 isAllocated_ =
false;
434 Graph_ = Teuchos::null;
436 A_local_ = makeLocalFilter (A_);
438 A_local_.is_null (), std::logic_error,
"Ifpack2::RILUK::initialize: "
439 "makeLocalFilter returned null; it failed to compute A_local. "
440 "Please report this bug to the Ifpack2 developers.");
448 RCP<const crs_matrix_type> A_local_crs =
450 if (A_local_crs.is_null ()) {
451 local_ordinal_type numRows = A_local_->getNodeNumRows();
452 Array<size_t> entriesPerRow(numRows);
453 for(local_ordinal_type i = 0; i < numRows; i++)
455 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
457 RCP<crs_matrix_type> A_local_crs_nc =
459 A_local_->getColMap (),
464 for(local_ordinal_type i = 0; i < numRows; i++)
466 size_t numEntries = 0;
467 A_local_->getLocalRowCopy(i, indices(), values(), numEntries);
468 ArrayView<const local_ordinal_type> indicesInsert(indices.data(), numEntries);
469 ArrayView<const scalar_type> valuesInsert(values.data(), numEntries);
470 A_local_crs_nc->insertLocalValues(i, indicesInsert, valuesInsert);
472 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
476 LevelOfFill_, 0, Overalloc_));
479 Graph_->initialize ();
481 checkOrderingConsistency (*A_local_);
482 L_solver_->setMatrix (L_);
483 L_solver_->initialize ();
484 U_solver_->setMatrix (U_);
485 U_solver_->initialize ();
492 isInitialized_ =
true;
494 initializeTime_ += timer.totalElapsedTime ();
497 template<
class MatrixType>
507 bool gidsAreConsistentlyOrdered=
true;
508 global_ordinal_type indexOfInconsistentGID=0;
509 for (global_ordinal_type i=0; i<rowGIDs.
size(); ++i) {
510 if (rowGIDs[i] != colGIDs[i]) {
511 gidsAreConsistentlyOrdered=
false;
512 indexOfInconsistentGID=i;
517 "The ordering of the local GIDs in the row and column maps is not the same"
518 << std::endl <<
"at index " << indexOfInconsistentGID
519 <<
". Consistency is required, as all calculations are done with"
520 << std::endl <<
"local indexing.");
523 template<
class MatrixType>
526 initAllValues (
const row_matrix_type& A)
533 using Teuchos::reduceAll;
534 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
536 size_t NumIn = 0, NumL = 0, NumU = 0;
537 bool DiagFound =
false;
538 size_t NumNonzeroDiags = 0;
539 size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
551 const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
556 L_->setAllToScalar (STS::zero ());
557 U_->setAllToScalar (STS::zero ());
560 D_->putScalar (STS::zero ());
561 ArrayRCP<scalar_type> DV = D_->get1dViewNonConst ();
563 RCP<const map_type> rowMap = L_->getRowMap ();
574 for (
size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
575 local_ordinal_type local_row = myRow;
579 A.getLocalRowCopy (local_row, InI(), InV(), NumIn);
587 for (
size_t j = 0; j < NumIn; ++j) {
588 const local_ordinal_type k = InI[j];
590 if (k == local_row) {
593 DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
597 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current "
598 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is "
599 "probably an artifact of the undocumented assumptions of the "
600 "original implementation (likely copied and pasted from Ifpack). "
601 "Nevertheless, the code I found here insisted on this being an error "
602 "state, so I will throw an exception here.");
604 else if (k < local_row) {
609 else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
621 DV[local_row] = Athresh_;
626 L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
630 L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
636 U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
640 U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
648 isInitialized_ =
true;
652 template<
class MatrixType>
656 const char prefix[] =
"Ifpack2::RILUK::compute: ";
662 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
663 "call setMatrix() with a nonnull input before calling this method.");
665 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
666 "fill complete. You may not invoke initialize() or compute() with this "
667 "matrix until the matrix is fill complete. If your matrix is a "
668 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
669 "range Maps, if appropriate) before calling this method.");
671 if (! isInitialized ()) {
684 initAllValues (*A_local_);
690 const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
692 size_t NumIn, NumL, NumU;
695 const size_t MaxNumEntries =
696 L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
700 size_t num_cols = U_->getColMap()->getNodeNumElements();
711 for (
size_t j = 0; j < num_cols; ++j) {
715 for (
size_t i = 0; i < L_->getNodeNumRows (); ++i) {
720 NumIn = MaxNumEntries;
721 L_->getLocalRowCopy (local_row, InI (), InV (), NumL);
724 InI[NumL] = local_row;
726 U_->getLocalRowCopy (local_row, InI (NumL+1, MaxNumEntries-NumL-1),
727 InV (NumL+1, MaxNumEntries-NumL-1), NumU);
731 for (
size_t j = 0; j < NumIn; ++j) {
737 for (
size_t jj = 0; jj < NumL; ++jj) {
743 U_->getLocalRowView(j, UUI, UUV);
746 if (RelaxValue_ == STM::zero ()) {
747 for (
size_t k = 0; k < NumUU; ++k) {
748 const int kk = colflag[UUI[k]];
753 InV[kk] -= multiplier * UUV[k];
758 for (
size_t k = 0; k < NumUU; ++k) {
762 const int kk = colflag[UUI[k]];
764 InV[kk] -= multiplier*UUV[k];
767 diagmod -= multiplier*UUV[k];
774 L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
779 if (RelaxValue_ != STM::zero ()) {
780 DV[i] += RelaxValue_*diagmod;
783 if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
784 if (STS::real (DV[i]) < STM::zero ()) {
785 DV[i] = -MinDiagonalValue;
788 DV[i] = MinDiagonalValue;
792 DV[i] = STS::one () / DV[i];
795 for (
size_t j = 0; j < NumU; ++j) {
796 InV[NumL+1+j] *= DV[i];
801 U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
805 for (
size_t j = 0; j < NumIn; ++j) {
806 colflag[InI[j]] = -1;
817 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
818 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
821 L_solver_->setMatrix (L_);
822 L_solver_->compute ();
823 U_solver_->setMatrix (U_);
824 U_solver_->compute ();
832 template<
class MatrixType>
835 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
836 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
842 using Teuchos::rcpFromRef;
845 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::apply: The matrix is "
846 "null. Please call setMatrix() with a nonnull input, then initialize() "
847 "and compute(), before calling this method.");
849 ! isComputed (), std::runtime_error,
850 "Ifpack2::RILUK::apply: If you have not yet called compute(), "
851 "you must call compute() before calling this method.");
852 TEUCHOS_TEST_FOR_EXCEPTION(
853 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
854 "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. "
855 "X.getNumVectors() = " << X.getNumVectors ()
856 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
857 TEUCHOS_TEST_FOR_EXCEPTION(
859 "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
860 "complex Scalar type. Please talk to the Ifpack2 developers to get this "
861 "fixed. There is a FIXME in this file about this very issue.");
862 #ifdef HAVE_IFPACK2_DEBUG
865 TEUCHOS_TEST_FOR_EXCEPTION( STM::isnaninf (D_nrm1), std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the stored diagonal is NaN or Inf.");
869 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
870 if (STM::isnaninf (norms[j])) {
875 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
877 #endif // HAVE_IFPACK2_DEBUG
885 if (alpha == one && beta == zero) {
888 L_solver_->apply (X, Y, mode);
892 Y.elementWiseMultiply (one, *D_, Y, zero);
894 U_solver_->apply (Y, Y, mode);
899 U_solver_->apply (X, Y, mode);
906 Y.elementWiseMultiply (one, *D_, Y, zero);
908 L_solver_->apply (Y, Y, mode);
922 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
923 apply (X, Y_tmp, mode);
924 Y.update (alpha, Y_tmp, beta);
929 #ifdef HAVE_IFPACK2_DEBUG
934 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
935 if (STM::isnaninf (norms[j])) {
940 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
942 #endif // HAVE_IFPACK2_DEBUG
949 template<
class MatrixType>
951 multiply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
952 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
955 const scalar_type zero = STS::zero ();
956 const scalar_type one = STS::one ();
959 U_->apply (X, Y, mode);
960 Y.update (one, X, one);
965 Y.elementWiseMultiply (one, *D_, Y, zero);
968 L_->apply (Y_tmp, Y, mode);
969 Y.update (one, Y_tmp, one);
972 L_->apply (X, Y, mode);
973 Y.update (one, X, one);
974 Y.elementWiseMultiply (one, *D_, Y, zero);
976 U_->apply (Y_tmp, Y, mode);
977 Y.update (one, Y_tmp, one);
982 template<
class MatrixType>
985 std::ostringstream os;
990 os <<
"\"Ifpack2::RILUK\": {";
991 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
992 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
994 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
997 os <<
"Matrix: null";
1000 os <<
"Global matrix dimensions: ["
1001 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
1002 <<
", Global nnz: " << A_->getGlobalNumEntries();
1005 if (! L_solver_.is_null ()) os <<
", " << L_solver_->description ();
1006 if (! U_solver_.is_null ()) os <<
", " << U_solver_->description ();
1014 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \
1015 template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_RILUK_def.hpp:224
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:261
const Tpetra::Vector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > & getD() const
Return the diagonal entries of the ILU factorization.
Definition: Ifpack2_RILUK_def.hpp:160
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_RILUK_def.hpp:396
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_RILUK_decl.hpp:273
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_RILUK_def.hpp:107
std::string description() const
A one-line description of this object.
Definition: Ifpack2_RILUK_def.hpp:983
void setParameters(const Teuchos::ParameterList ¶ms)
Definition: Ifpack2_RILUK_def.hpp:271
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:243
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_RILUK_def.hpp:187
Teuchos::RCP< const crs_matrix_type > getCrsMatrix() const
Return the input matrix A as a Tpetra::CrsMatrix, if possible; else throws.
Definition: Ifpack2_RILUK_def.hpp:356
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Tpetra::CrsMatrix specialization used by this class for representing L and U.
Definition: Ifpack2_RILUK_decl.hpp:276
const crs_matrix_type & getU() const
Return the U factor of the ILU factorization.
Definition: Ifpack2_RILUK_def.hpp:174
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:264
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_RILUK_def.hpp:205
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:258
virtual ~RILUK()
Destructor (declared virtual for memory safety).
Definition: Ifpack2_RILUK_def.hpp:96
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:77
void compute()
Compute the (numeric) incomplete factorization.
Definition: Ifpack2_RILUK_def.hpp:653
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the (inverse of the) incomplete factorization to X, resulting in Y.
Definition: Ifpack2_RILUK_def.hpp:835
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_RILUK_decl.hpp:267
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:97
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the input matrix.
Definition: Ifpack2_RILUK_def.hpp:349
double totalElapsedTime(bool readCurrentTime=false) const
const crs_matrix_type & getL() const
Return the L factor of the ILU factorization.
Definition: Ifpack2_RILUK_def.hpp:143
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:255