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>
56 isInitialized_ (false),
61 initializeTime_ (0.0),
72 template<
class MatrixType>
77 isInitialized_ (false),
82 initializeTime_ (0.0),
93 template<
class MatrixType>
96 template<
class MatrixType>
103 template<
class MatrixType>
112 isAllocated_ =
false;
113 isInitialized_ =
false;
115 A_local_ = Teuchos::null;
116 Graph_ = Teuchos::null;
123 if (! L_solver_.is_null ()) {
124 L_solver_->setMatrix (Teuchos::null);
126 if (! U_solver_.is_null ()) {
127 U_solver_->setMatrix (Teuchos::null);
139 template<
class MatrixType>
144 L_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor "
145 "is null. Please call initialize() (and preferably also compute()) "
146 "before calling this method. If the input matrix has not yet been set, "
147 "you must first call setMatrix() with a nonnull input matrix before you "
148 "may call initialize() or compute().");
153 template<
class MatrixType>
154 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
161 D_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor "
162 "(of diagonal entries) is null. Please call initialize() (and "
163 "preferably also compute()) before calling this method. If the input "
164 "matrix has not yet been set, you must first call setMatrix() with a "
165 "nonnull input matrix before you may call initialize() or compute().");
170 template<
class MatrixType>
175 U_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U factor "
176 "is null. Please call initialize() (and preferably also compute()) "
177 "before calling this method. If the input matrix has not yet been set, "
178 "you must first call setMatrix() with a nonnull input matrix before you "
179 "may call initialize() or compute().");
184 template<
class MatrixType>
187 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getNodeSmootherComplexity: "
188 "The input matrix A is null. Please call setMatrix() with a nonnull "
189 "input matrix, then call compute(), before calling this method.");
191 if(!L_.is_null() && !U_.is_null())
192 return A_->getNodeNumEntries() + L_->getNodeNumEntries() + U_->getNodeNumEntries();
199 template<
class MatrixType>
205 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
206 "The matrix is null. Please call setMatrix() with a nonnull input "
207 "before calling this method.");
211 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
212 "The computed graph is null. Please call initialize() before calling "
214 return Graph_->getL_Graph ()->getDomainMap ();
218 template<
class MatrixType>
224 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
225 "The matrix is null. Please call setMatrix() with a nonnull input "
226 "before calling this method.");
230 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
231 "The computed graph is null. Please call initialize() before calling "
233 return Graph_->getL_Graph ()->getRangeMap ();
237 template<
class MatrixType>
243 if (! isAllocated_) {
252 L_ =
rcp (
new crs_matrix_type (Graph_->getL_Graph ()));
253 U_ =
rcp (
new crs_matrix_type (Graph_->getU_Graph ()));
254 L_->setAllToScalar (STS::zero ());
255 U_->setAllToScalar (STS::zero ());
260 D_ =
rcp (
new vec_type (Graph_->getL_Graph ()->getRowMap ()));
266 template<
class MatrixType>
271 using Details::getParamTryingTypes;
272 const char prefix[] =
"Ifpack2::RILUK: ";
289 const std::string paramName (
"fact: iluk level-of-fill");
290 getParamTryingTypes<int, int, global_ordinal_type, double, float>
291 (fillLevel, params, paramName, prefix);
296 const std::string paramName (
"fact: absolute threshold");
297 getParamTryingTypes<magnitude_type, magnitude_type, double>
298 (absThresh, params, paramName, prefix);
301 const std::string paramName (
"fact: relative threshold");
302 getParamTryingTypes<magnitude_type, magnitude_type, double>
303 (relThresh, params, paramName, prefix);
306 const std::string paramName (
"fact: relax value");
307 getParamTryingTypes<magnitude_type, magnitude_type, double>
308 (relaxValue, params, paramName, prefix);
312 L_solver_->setParameters(params);
313 U_solver_->setParameters(params);
319 LevelOfFill_ = fillLevel;
326 if (absThresh != -STM::one ()) {
327 Athresh_ = absThresh;
329 if (relThresh != -STM::one ()) {
330 Rthresh_ = relThresh;
332 if (relaxValue != -STM::one ()) {
333 RelaxValue_ = relaxValue;
338 template<
class MatrixType>
345 template<
class MatrixType>
352 template<
class MatrixType>
358 using Teuchos::rcp_dynamic_cast;
359 using Teuchos::rcp_implicit_cast;
364 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
365 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
372 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
373 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
374 if (! A_lf_r.is_null ()) {
375 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
381 return rcp (
new LocalFilter<row_matrix_type> (A));
386 template<
class MatrixType>
391 using Teuchos::rcp_const_cast;
392 using Teuchos::rcp_dynamic_cast;
393 using Teuchos::rcp_implicit_cast;
399 const char prefix[] =
"Ifpack2::RILUK::initialize: ";
402 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
403 "call setMatrix() with a nonnull input before calling this method.");
405 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
406 "fill complete. You may not invoke initialize() or compute() with this "
407 "matrix until the matrix is fill complete. If your matrix is a "
408 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
409 "range Maps, if appropriate) before calling this method.");
422 isInitialized_ =
false;
423 isAllocated_ =
false;
425 Graph_ = Teuchos::null;
427 A_local_ = makeLocalFilter (A_);
429 A_local_.is_null (), std::logic_error,
"Ifpack2::RILUK::initialize: "
430 "makeLocalFilter returned null; it failed to compute A_local. "
431 "Please report this bug to the Ifpack2 developers.");
439 RCP<const crs_matrix_type> A_local_crs =
441 if (A_local_crs.is_null ()) {
442 local_ordinal_type numRows = A_local_->getNodeNumRows();
443 Array<size_t> entriesPerRow(numRows);
444 for(local_ordinal_type i = 0; i < numRows; i++)
446 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
448 RCP<crs_matrix_type> A_local_crs_nc =
450 A_local_->getColMap (),
455 for(local_ordinal_type i = 0; i < numRows; i++)
457 size_t numEntries = 0;
458 A_local_->getLocalRowCopy(i, indices(), values(), numEntries);
459 ArrayView<const local_ordinal_type> indicesInsert(indices.data(), numEntries);
460 ArrayView<const scalar_type> valuesInsert(values.data(), numEntries);
461 A_local_crs_nc->insertLocalValues(i, indicesInsert, valuesInsert);
463 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
470 Graph_->initialize ();
472 checkOrderingConsistency (*A_local_);
473 L_solver_->setMatrix (L_);
474 L_solver_->initialize ();
475 U_solver_->setMatrix (U_);
476 U_solver_->initialize ();
483 isInitialized_ =
true;
485 initializeTime_ += timer.totalElapsedTime ();
488 template<
class MatrixType>
498 bool gidsAreConsistentlyOrdered=
true;
499 global_ordinal_type indexOfInconsistentGID=0;
500 for (global_ordinal_type i=0; i<rowGIDs.
size(); ++i) {
501 if (rowGIDs[i] != colGIDs[i]) {
502 gidsAreConsistentlyOrdered=
false;
503 indexOfInconsistentGID=i;
508 "The ordering of the local GIDs in the row and column maps is not the same"
509 << std::endl <<
"at index " << indexOfInconsistentGID
510 <<
". Consistency is required, as all calculations are done with"
511 << std::endl <<
"local indexing.");
514 template<
class MatrixType>
517 initAllValues (
const row_matrix_type& A)
525 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
527 size_t NumIn = 0, NumL = 0, NumU = 0;
528 bool DiagFound =
false;
529 size_t NumNonzeroDiags = 0;
530 size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
542 const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
547 L_->setAllToScalar (STS::zero ());
548 U_->setAllToScalar (STS::zero ());
551 D_->putScalar (STS::zero ());
552 ArrayRCP<scalar_type> DV = D_->get1dViewNonConst ();
554 RCP<const map_type> rowMap = L_->getRowMap ();
565 for (
size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
566 local_ordinal_type local_row = myRow;
570 A.getLocalRowCopy (local_row, InI(), InV(), NumIn);
578 for (
size_t j = 0; j < NumIn; ++j) {
579 const local_ordinal_type k = InI[j];
581 if (k == local_row) {
584 DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
588 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current "
589 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is "
590 "probably an artifact of the undocumented assumptions of the "
591 "original implementation (likely copied and pasted from Ifpack). "
592 "Nevertheless, the code I found here insisted on this being an error "
593 "state, so I will throw an exception here.");
595 else if (k < local_row) {
600 else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
612 DV[local_row] = Athresh_;
617 L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
621 L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
627 U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
631 U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
639 isInitialized_ =
true;
643 template<
class MatrixType>
647 const char prefix[] =
"Ifpack2::RILUK::compute: ";
653 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
654 "call setMatrix() with a nonnull input before calling this method.");
656 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
657 "fill complete. You may not invoke initialize() or compute() with this "
658 "matrix until the matrix is fill complete. If your matrix is a "
659 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
660 "range Maps, if appropriate) before calling this method.");
662 if (! isInitialized ()) {
675 initAllValues (*A_local_);
681 const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
683 size_t NumIn, NumL, NumU;
686 const size_t MaxNumEntries =
687 L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
691 size_t num_cols = U_->getColMap()->getNodeNumElements();
702 for (
size_t j = 0; j < num_cols; ++j) {
706 for (
size_t i = 0; i < L_->getNodeNumRows (); ++i) {
711 NumIn = MaxNumEntries;
712 L_->getLocalRowCopy (local_row, InI (), InV (), NumL);
715 InI[NumL] = local_row;
717 U_->getLocalRowCopy (local_row, InI (NumL+1, MaxNumEntries-NumL-1),
718 InV (NumL+1, MaxNumEntries-NumL-1), NumU);
722 for (
size_t j = 0; j < NumIn; ++j) {
728 for (
size_t jj = 0; jj < NumL; ++jj) {
734 U_->getLocalRowView(j, UUI, UUV);
737 if (RelaxValue_ == STM::zero ()) {
738 for (
size_t k = 0; k < NumUU; ++k) {
739 const int kk = colflag[UUI[k]];
744 InV[kk] -= multiplier * UUV[k];
749 for (
size_t k = 0; k < NumUU; ++k) {
753 const int kk = colflag[UUI[k]];
755 InV[kk] -= multiplier*UUV[k];
758 diagmod -= multiplier*UUV[k];
765 L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
770 if (RelaxValue_ != STM::zero ()) {
771 DV[i] += RelaxValue_*diagmod;
774 if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
775 if (STS::real (DV[i]) < STM::zero ()) {
776 DV[i] = -MinDiagonalValue;
779 DV[i] = MinDiagonalValue;
783 DV[i] = STS::one () / DV[i];
786 for (
size_t j = 0; j < NumU; ++j) {
787 InV[NumL+1+j] *= DV[i];
792 U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
796 for (
size_t j = 0; j < NumIn; ++j) {
797 colflag[InI[j]] = -1;
808 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
809 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
812 L_solver_->setMatrix (L_);
813 L_solver_->compute ();
814 U_solver_->setMatrix (U_);
815 U_solver_->compute ();
823 template<
class MatrixType>
826 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
827 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
833 using Teuchos::rcpFromRef;
836 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::apply: The matrix is "
837 "null. Please call setMatrix() with a nonnull input, then initialize() "
838 "and compute(), before calling this method.");
840 ! isComputed (), std::runtime_error,
841 "Ifpack2::RILUK::apply: If you have not yet called compute(), "
842 "you must call compute() before calling this method.");
843 TEUCHOS_TEST_FOR_EXCEPTION(
844 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
845 "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. "
846 "X.getNumVectors() = " << X.getNumVectors ()
847 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
848 TEUCHOS_TEST_FOR_EXCEPTION(
850 "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
851 "complex Scalar type. Please talk to the Ifpack2 developers to get this "
852 "fixed. There is a FIXME in this file about this very issue.");
853 #ifdef HAVE_IFPACK2_DEBUG
856 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.");
860 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
861 if (STM::isnaninf (norms[j])) {
866 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
868 #endif // HAVE_IFPACK2_DEBUG
876 if (alpha == one && beta == zero) {
879 L_solver_->apply (X, Y, mode);
883 Y.elementWiseMultiply (one, *D_, Y, zero);
885 U_solver_->apply (Y, Y, mode);
890 U_solver_->apply (X, Y, mode);
897 Y.elementWiseMultiply (one, *D_, Y, zero);
899 L_solver_->apply (Y, Y, mode);
913 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
914 apply (X, Y_tmp, mode);
915 Y.update (alpha, Y_tmp, beta);
920 #ifdef HAVE_IFPACK2_DEBUG
925 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
926 if (STM::isnaninf (norms[j])) {
931 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
933 #endif // HAVE_IFPACK2_DEBUG
940 template<
class MatrixType>
942 multiply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
943 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
946 const scalar_type zero = STS::zero ();
947 const scalar_type one = STS::one ();
950 U_->apply (X, Y, mode);
951 Y.update (one, X, one);
956 Y.elementWiseMultiply (one, *D_, Y, zero);
959 L_->apply (Y_tmp, Y, mode);
960 Y.update (one, Y_tmp, one);
963 L_->apply (X, Y, mode);
964 Y.update (one, X, one);
965 Y.elementWiseMultiply (one, *D_, Y, zero);
967 U_->apply (Y_tmp, Y, mode);
968 Y.update (one, Y_tmp, one);
973 template<
class MatrixType>
976 std::ostringstream os;
981 os <<
"\"Ifpack2::RILUK\": {";
982 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
983 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
985 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
988 os <<
"Matrix: null";
991 os <<
"Global matrix dimensions: ["
992 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
993 <<
", Global nnz: " << A_->getGlobalNumEntries();
996 if (! L_solver_.is_null ()) os <<
", " << L_solver_->description ();
997 if (! U_solver_.is_null ()) os <<
", " << U_solver_->description ();
1005 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \
1006 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:222
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:158
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_RILUK_def.hpp:387
#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:105
std::string description() const
A one-line description of this object.
Definition: Ifpack2_RILUK_def.hpp:974
void setParameters(const Teuchos::ParameterList ¶ms)
Definition: Ifpack2_RILUK_def.hpp:269
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:185
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:347
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:172
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:203
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:94
"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:644
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
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:826
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:340
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:141
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:255