41#ifndef IFPACK2_CRSRILUK_DEF_HPP
42#define IFPACK2_CRSRILUK_DEF_HPP
44#include "Ifpack2_LocalFilter.hpp"
45#include "Tpetra_CrsMatrix.hpp"
46#include "Teuchos_StandardParameterEntryValidators.hpp"
47#include "Ifpack2_LocalSparseTriangularSolver.hpp"
48#include "Ifpack2_Details_getParamTryingTypes.hpp"
49#include "Ifpack2_Details_getCrsMatrix.hpp"
50#include "Kokkos_Sort.hpp"
51#include "KokkosSparse_Utils.hpp"
52#include "KokkosKernels_Sorting.hpp"
63 static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
65 type_strs[0] =
"Serial";
66 type_strs[1] =
"KSPILUK";
68 type_enums[0] = Serial;
69 type_enums[1] = KSPILUK;
74template<
class MatrixType>
80 isInitialized_ (false),
85 initializeTime_ (0.0),
91 isKokkosKernelsSpiluk_(false)
97template<
class MatrixType>
102 isAllocated_ (false),
103 isInitialized_ (false),
108 initializeTime_ (0.0),
114 isKokkosKernelsSpiluk_(false)
120template<
class MatrixType>
123 if (Teuchos::nonnull (KernelHandle_))
125 KernelHandle_->destroy_spiluk_handle();
129template<
class MatrixType>
133 L_solver_->setObjectLabel(
"lower");
135 U_solver_->setObjectLabel(
"upper");
138template<
class MatrixType>
146 if (A.getRawPtr () != A_.getRawPtr ()) {
147 isAllocated_ =
false;
148 isInitialized_ =
false;
150 A_local_ = Teuchos::null;
151 Graph_ = Teuchos::null;
158 if (! L_solver_.is_null ()) {
159 L_solver_->setMatrix (Teuchos::null);
161 if (! U_solver_.is_null ()) {
162 U_solver_->setMatrix (Teuchos::null);
174template<
class MatrixType>
178 TEUCHOS_TEST_FOR_EXCEPTION(
179 L_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor "
180 "is null. Please call initialize() (and preferably also compute()) "
181 "before calling this method. If the input matrix has not yet been set, "
182 "you must first call setMatrix() with a nonnull input matrix before you "
183 "may call initialize() or compute().");
188template<
class MatrixType>
189const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
195 TEUCHOS_TEST_FOR_EXCEPTION(
196 D_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor "
197 "(of diagonal entries) is null. Please call initialize() (and "
198 "preferably also compute()) before calling this method. If the input "
199 "matrix has not yet been set, you must first call setMatrix() with a "
200 "nonnull input matrix before you may call initialize() or compute().");
205template<
class MatrixType>
209 TEUCHOS_TEST_FOR_EXCEPTION(
210 U_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U factor "
211 "is null. Please call initialize() (and preferably also compute()) "
212 "before calling this method. If the input matrix has not yet been set, "
213 "you must first call setMatrix() with a nonnull input matrix before you "
214 "may call initialize() or compute().");
219template<
class MatrixType>
221 TEUCHOS_TEST_FOR_EXCEPTION(
222 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getNodeSmootherComplexity: "
223 "The input matrix A is null. Please call setMatrix() with a nonnull "
224 "input matrix, then call compute(), before calling this method.");
226 if(!L_.is_null() && !U_.is_null())
227 return A_->getLocalNumEntries() + L_->getLocalNumEntries() + U_->getLocalNumEntries();
233template<
class MatrixType>
234Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
238 TEUCHOS_TEST_FOR_EXCEPTION(
239 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
240 "The matrix is null. Please call setMatrix() with a nonnull input "
241 "before calling this method.");
244 TEUCHOS_TEST_FOR_EXCEPTION(
245 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: "
246 "The computed graph is null. Please call initialize() before calling "
248 return Graph_->getL_Graph ()->getDomainMap ();
252template<
class MatrixType>
253Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
257 TEUCHOS_TEST_FOR_EXCEPTION(
258 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
259 "The matrix is null. Please call setMatrix() with a nonnull input "
260 "before calling this method.");
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: "
265 "The computed graph is null. Please call initialize() before calling "
267 return Graph_->getL_Graph ()->getRangeMap ();
271template<
class MatrixType>
277 if (! isAllocated_) {
286 L_ = rcp (
new crs_matrix_type (Graph_->getL_Graph ()));
287 U_ = rcp (
new crs_matrix_type (Graph_->getU_Graph ()));
288 L_->setAllToScalar (STS::zero ());
289 U_->setAllToScalar (STS::zero ());
294 D_ = rcp (
new vec_type (Graph_->getL_Graph ()->getRowMap ()));
300template<
class MatrixType>
306 using Teuchos::ParameterList;
307 using Teuchos::Array;
308 using Details::getParamTryingTypes;
309 const char prefix[] =
"Ifpack2::RILUK: ";
316 double overalloc = 2.;
327 const std::string paramName (
"fact: iluk level-of-fill");
328 getParamTryingTypes<int, int, global_ordinal_type, double, float>
329 (fillLevel, params, paramName, prefix);
334 const std::string paramName (
"fact: absolute threshold");
335 getParamTryingTypes<magnitude_type, magnitude_type, double>
336 (absThresh, params, paramName, prefix);
339 const std::string paramName (
"fact: relative threshold");
340 getParamTryingTypes<magnitude_type, magnitude_type, double>
341 (relThresh, params, paramName, prefix);
344 const std::string paramName (
"fact: relax value");
345 getParamTryingTypes<magnitude_type, magnitude_type, double>
346 (relaxValue, params, paramName, prefix);
349 const std::string paramName (
"fact: iluk overalloc");
350 getParamTryingTypes<double, double>
351 (overalloc, params, paramName, prefix);
355 Details::IlukImplType::Enum ilukimplType = Details::IlukImplType::Serial;
357 static const char typeName[] =
"fact: type";
359 if ( ! params.isType<std::string>(typeName))
break;
362 Array<std::string> ilukimplTypeStrs;
363 Array<Details::IlukImplType::Enum> ilukimplTypeEnums;
364 Details::IlukImplType::loadPLTypeOption (ilukimplTypeStrs, ilukimplTypeEnums);
365 Teuchos::StringToIntegralParameterEntryValidator<Details::IlukImplType::Enum>
366 s2i(ilukimplTypeStrs (), ilukimplTypeEnums (), typeName,
false);
368 ilukimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
371 if (ilukimplType == Details::IlukImplType::KSPILUK) {
372 this->isKokkosKernelsSpiluk_ =
true;
375 this->isKokkosKernelsSpiluk_ =
false;
379 L_solver_->setParameters(params);
380 U_solver_->setParameters(params);
386 LevelOfFill_ = fillLevel;
387 Overalloc_ = overalloc;
394 if (absThresh != -STM::one ()) {
395 Athresh_ = absThresh;
397 if (relThresh != -STM::one ()) {
398 Rthresh_ = relThresh;
400 if (relaxValue != -STM::one ()) {
401 RelaxValue_ = relaxValue;
406template<
class MatrixType>
407Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
409 return Teuchos::rcp_implicit_cast<const row_matrix_type> (A_);
413template<
class MatrixType>
414Teuchos::RCP<const typename RILUK<MatrixType>::crs_matrix_type>
416 return Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_,
true);
420template<
class MatrixType>
421Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
426 using Teuchos::rcp_dynamic_cast;
427 using Teuchos::rcp_implicit_cast;
432 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
433 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
440 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
441 rcp_dynamic_cast<const LocalFilter<row_matrix_type> > (A);
442 if (! A_lf_r.is_null ()) {
443 return rcp_implicit_cast<const row_matrix_type> (A_lf_r);
454template<
class MatrixType>
459 using Teuchos::rcp_const_cast;
460 using Teuchos::rcp_dynamic_cast;
461 using Teuchos::rcp_implicit_cast;
462 using Teuchos::Array;
463 using Teuchos::ArrayView;
467 const char prefix[] =
"Ifpack2::RILUK::initialize: ";
469 TEUCHOS_TEST_FOR_EXCEPTION
470 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
471 "call setMatrix() with a nonnull input before calling this method.");
472 TEUCHOS_TEST_FOR_EXCEPTION
473 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
474 "fill complete. You may not invoke initialize() or compute() with this "
475 "matrix until the matrix is fill complete. If your matrix is a "
476 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
477 "range Maps, if appropriate) before calling this method.");
479 Teuchos::Time timer (
"RILUK::initialize");
480 double startTime = timer.wallTime();
482 Teuchos::TimeMonitor timeMon (timer);
491 isInitialized_ =
false;
492 isAllocated_ =
false;
494 Graph_ = Teuchos::null;
496 A_local_ = makeLocalFilter (A_);
497 TEUCHOS_TEST_FOR_EXCEPTION(
498 A_local_.is_null (), std::logic_error,
"Ifpack2::RILUK::initialize: "
499 "makeLocalFilter returned null; it failed to compute A_local. "
500 "Please report this bug to the Ifpack2 developers.");
508 if (this->isKokkosKernelsSpiluk_) {
509 this->KernelHandle_ = Teuchos::rcp (
new kk_handle_type ());
510 KernelHandle_->create_spiluk_handle( KokkosSparse::Experimental::SPILUKAlgorithm::SEQLVLSCHD_TP1,
511 A_local_->getLocalNumRows(),
512 2*A_local_->getLocalNumEntries()*(LevelOfFill_+1),
513 2*A_local_->getLocalNumEntries()*(LevelOfFill_+1) );
517 RCP<const crs_matrix_type> A_local_crs = Details::getCrsMatrix(A_local_);
518 if(A_local_crs.is_null()) {
519 local_ordinal_type numRows = A_local_->getLocalNumRows();
520 Array<size_t> entriesPerRow(numRows);
522 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
524 RCP<crs_matrix_type> A_local_crs_nc =
525 rcp (
new crs_matrix_type (A_local_->getRowMap (),
526 A_local_->getColMap (),
529 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getLocalMaxNumRowEntries());
530 nonconst_values_host_view_type values(
"values",A_local_->getLocalMaxNumRowEntries());
531 for(local_ordinal_type i = 0; i < numRows; i++) {
532 size_t numEntries = 0;
533 A_local_->getLocalRowCopy(i, indices, values, numEntries);
534 A_local_crs_nc->insertLocalValues(i, numEntries,
reinterpret_cast<scalar_type*
>(values.data()), indices.data());
536 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
537 A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
540 LevelOfFill_, 0, Overalloc_));
543 if (this->isKokkosKernelsSpiluk_) Graph_->initialize (KernelHandle_);
544 else Graph_->initialize ();
547 checkOrderingConsistency (*A_local_);
548 L_solver_->setMatrix (L_);
549 L_solver_->initialize ();
553#if !defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) || !defined(KOKKOS_ENABLE_CUDA) || (CUDA_VERSION < 11030)
554 L_solver_->compute ();
556 U_solver_->setMatrix (U_);
557 U_solver_->initialize ();
558#if !defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) || !defined(KOKKOS_ENABLE_CUDA) || (CUDA_VERSION < 11030)
559 U_solver_->compute ();
567 isInitialized_ =
true;
569 initializeTime_ += (timer.wallTime() - startTime);
572template<
class MatrixType>
580 Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getLocalElementList();
581 Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getLocalElementList();
582 bool gidsAreConsistentlyOrdered=
true;
583 global_ordinal_type indexOfInconsistentGID=0;
584 for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
585 if (rowGIDs[i] != colGIDs[i]) {
586 gidsAreConsistentlyOrdered=
false;
587 indexOfInconsistentGID=i;
591 TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==
false, std::runtime_error,
592 "The ordering of the local GIDs in the row and column maps is not the same"
593 << std::endl <<
"at index " << indexOfInconsistentGID
594 <<
". Consistency is required, as all calculations are done with"
595 << std::endl <<
"local indexing.");
598template<
class MatrixType>
601initAllValues (
const row_matrix_type& A)
603 using Teuchos::ArrayRCP;
607 using Teuchos::REDUCE_SUM;
608 using Teuchos::reduceAll;
609 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
611 size_t NumIn = 0, NumL = 0, NumU = 0;
612 bool DiagFound =
false;
613 size_t NumNonzeroDiags = 0;
614 size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
618 nonconst_local_inds_host_view_type InI(
"InI",MaxNumEntries);
619 Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
620 Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
621 nonconst_values_host_view_type InV(
"InV",MaxNumEntries);
622 Teuchos::Array<scalar_type> LV(MaxNumEntries);
623 Teuchos::Array<scalar_type> UV(MaxNumEntries);
626 const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
631 L_->setAllToScalar (STS::zero ());
632 U_->setAllToScalar (STS::zero ());
635 D_->putScalar (STS::zero ());
636 auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
638 RCP<const map_type> rowMap = L_->getRowMap ();
648 Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getLocalElementList();
649 for (
size_t myRow=0; myRow<A.getLocalNumRows(); ++myRow) {
650 local_ordinal_type local_row = myRow;
654 A.getLocalRowCopy (local_row, InI, InV, NumIn);
662 for (
size_t j = 0; j < NumIn; ++j) {
663 const local_ordinal_type k = InI[j];
665 if (k == local_row) {
668 DV(local_row) += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
671 TEUCHOS_TEST_FOR_EXCEPTION(
672 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current "
673 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is "
674 "probably an artifact of the undocumented assumptions of the "
675 "original implementation (likely copied and pasted from Ifpack). "
676 "Nevertheless, the code I found here insisted on this being an error "
677 "state, so I will throw an exception here.");
679 else if (k < local_row) {
684 else if (Teuchos::as<size_t>(k) <= rowMap->getLocalNumElements()) {
696 DV(local_row) = Athresh_;
701 L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
705 L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
711 U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
715 U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
723 isInitialized_ =
true;
727template<
class MatrixType>
732 using Teuchos::rcp_const_cast;
733 using Teuchos::rcp_dynamic_cast;
734 using Teuchos::Array;
735 using Teuchos::ArrayView;
736 const char prefix[] =
"Ifpack2::RILUK::compute: ";
741 TEUCHOS_TEST_FOR_EXCEPTION
742 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please "
743 "call setMatrix() with a nonnull input before calling this method.");
744 TEUCHOS_TEST_FOR_EXCEPTION
745 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not "
746 "fill complete. You may not invoke initialize() or compute() with this "
747 "matrix until the matrix is fill complete. If your matrix is a "
748 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
749 "range Maps, if appropriate) before calling this method.");
751 if (! isInitialized ()) {
755 Teuchos::Time timer (
"RILUK::compute");
758 Teuchos::TimeMonitor timeMon (timer);
759 double startTime = timer.wallTime();
763 if (!this->isKokkosKernelsSpiluk_) {
766 initAllValues (*A_local_);
772 const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
774 size_t NumIn, NumL, NumU;
777 const size_t MaxNumEntries =
778 L_->getLocalMaxNumRowEntries () + U_->getLocalMaxNumRowEntries () + 1;
780 Teuchos::Array<local_ordinal_type> InI(MaxNumEntries);
781 Teuchos::Array<scalar_type> InV(MaxNumEntries);
782 size_t num_cols = U_->getColMap()->getLocalNumElements();
783 Teuchos::Array<int> colflag(num_cols);
785 auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
789 for (
size_t j = 0; j < num_cols; ++j) {
792 using IST =
typename row_matrix_type::impl_scalar_type;
793 for (
size_t i = 0; i < L_->getLocalNumRows (); ++i) {
797 local_inds_host_view_type UUI;
798 values_host_view_type UUV;
802 NumIn = MaxNumEntries;
803 nonconst_local_inds_host_view_type InI_v(InI.data(),MaxNumEntries);
804 nonconst_values_host_view_type InV_v(
reinterpret_cast<IST*
>(InV.data()),MaxNumEntries);
806 L_->getLocalRowCopy (local_row, InI_v , InV_v, NumL);
809 InI[NumL] = local_row;
811 nonconst_local_inds_host_view_type InI_sub(InI.data()+NumL+1,MaxNumEntries-NumL-1);
812 nonconst_values_host_view_type InV_sub(
reinterpret_cast<IST*
>(InV.data())+NumL+1,MaxNumEntries-NumL-1);
814 U_->getLocalRowCopy (local_row, InI_sub,InV_sub, NumU);
818 for (
size_t j = 0; j < NumIn; ++j) {
824 for (
size_t jj = 0; jj < NumL; ++jj) {
826 IST multiplier = InV[jj];
830 U_->getLocalRowView(j, UUI, UUV);
833 if (RelaxValue_ == STM::zero ()) {
834 for (
size_t k = 0; k < NumUU; ++k) {
835 const int kk = colflag[UUI[k]];
840 InV[kk] -=
static_cast<scalar_type>(multiplier * UUV[k]);
846 for (
size_t k = 0; k < NumUU; ++k) {
850 const int kk = colflag[UUI[k]];
852 InV[kk] -=
static_cast<scalar_type>(multiplier*UUV[k]);
855 diagmod -=
static_cast<scalar_type>(multiplier*UUV[k]);
863 L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
868 if (RelaxValue_ != STM::zero ()) {
869 DV(i) += RelaxValue_*diagmod;
872 if (STS::magnitude (DV(i)) > STS::magnitude (MaxDiagonalValue)) {
873 if (STS::real (DV(i)) < STM::zero ()) {
874 DV(i) = -MinDiagonalValue;
877 DV(i) = MinDiagonalValue;
884 for (
size_t j = 0; j < NumU; ++j) {
890 U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
894 for (
size_t j = 0; j < NumIn; ++j) {
895 colflag[InI[j]] = -1;
906 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
907 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
910 L_solver_->setMatrix (L_);
911 L_solver_->compute ();
912 U_solver_->setMatrix (U_);
913 U_solver_->compute ();
917 RCP<const crs_matrix_type> A_local_crs = Details::getCrsMatrix(A_local_);
918 if(A_local_crs.is_null()) {
920 Array<size_t> entriesPerRow(numRows);
922 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
924 RCP<crs_matrix_type> A_local_crs_nc =
926 A_local_->getColMap (),
929 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getLocalMaxNumRowEntries());
930 nonconst_values_host_view_type values(
"values",A_local_->getLocalMaxNumRowEntries());
932 size_t numEntries = 0;
933 A_local_->getLocalRowCopy(i, indices, values, numEntries);
934 A_local_crs_nc->insertLocalValues(i, numEntries,
reinterpret_cast<scalar_type*
>(values.data()),indices.data());
936 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
937 A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
939 auto lclMtx = A_local_crs->getLocalMatrixDevice();
940 A_local_rowmap_ = lclMtx.graph.row_map;
941 A_local_entries_ = lclMtx.graph.entries;
942 A_local_values_ = lclMtx.values;
948 if (L_->isStaticGraph () || L_->isLocallyIndexed ()) {
949 L_->setAllToScalar (STS::zero ());
950 U_->setAllToScalar (STS::zero ());
953 using row_map_type =
typename crs_matrix_type::local_matrix_device_type::row_map_type;
956 auto lclL = L_->getLocalMatrixDevice();
957 row_map_type L_rowmap = lclL.graph.row_map;
958 auto L_entries = lclL.graph.entries;
959 auto L_values = lclL.values;
961 auto lclU = U_->getLocalMatrixDevice();
962 row_map_type U_rowmap = lclU.graph.row_map;
963 auto U_entries = lclU.graph.entries;
964 auto U_values = lclU.values;
966 KokkosSparse::Experimental::spiluk_numeric( KernelHandle_.getRawPtr(), LevelOfFill_,
967 A_local_rowmap_, A_local_entries_, A_local_values_,
968 L_rowmap, L_entries, L_values, U_rowmap, U_entries, U_values );
971 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
972 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
974 L_solver_->setMatrix (L_);
975 L_solver_->compute ();
976 U_solver_->setMatrix (U_);
977 U_solver_->compute ();
982 computeTime_ += (timer.wallTime() - startTime);
986template<
class MatrixType>
989apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
990 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
991 Teuchos::ETransp mode,
996 using Teuchos::rcpFromRef;
998 TEUCHOS_TEST_FOR_EXCEPTION(
999 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::apply: The matrix is "
1000 "null. Please call setMatrix() with a nonnull input, then initialize() "
1001 "and compute(), before calling this method.");
1002 TEUCHOS_TEST_FOR_EXCEPTION(
1003 ! isComputed (), std::runtime_error,
1004 "Ifpack2::RILUK::apply: If you have not yet called compute(), "
1005 "you must call compute() before calling this method.");
1006 TEUCHOS_TEST_FOR_EXCEPTION(
1007 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
1008 "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. "
1009 "X.getNumVectors() = " << X.getNumVectors ()
1010 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
1011 TEUCHOS_TEST_FOR_EXCEPTION(
1012 STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
1013 "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
1014 "complex Scalar type. Please talk to the Ifpack2 developers to get this "
1015 "fixed. There is a FIXME in this file about this very issue.");
1016#ifdef HAVE_IFPACK2_DEBUG
1019 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.");
1020 Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
1023 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
1024 if (STM::isnaninf (norms[j])) {
1029 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
1036 Teuchos::Time timer (
"RILUK::apply");
1037 double startTime = timer.wallTime();
1039 Teuchos::TimeMonitor timeMon (timer);
1040 if (alpha == one && beta == zero) {
1041 if (mode == Teuchos::NO_TRANS) {
1042#if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA) && (CUDA_VERSION >= 11030)
1046 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1049 L_solver_->apply (X, Y_tmp, mode);
1051 if (!this->isKokkosKernelsSpiluk_) {
1054 Y_tmp.elementWiseMultiply (one, *D_, Y_tmp, zero);
1057 U_solver_->apply (Y_tmp, Y, mode);
1060 L_solver_->apply (X, Y, mode);
1062 if (!this->isKokkosKernelsSpiluk_) {
1065 Y.elementWiseMultiply (one, *D_, Y, zero);
1068 U_solver_->apply (Y, Y, mode);
1072#if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) && defined(KOKKOS_ENABLE_CUDA) && (CUDA_VERSION >= 11030)
1076 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1079 U_solver_->apply (X, Y_tmp, mode);
1081 if (!this->isKokkosKernelsSpiluk_) {
1087 Y_tmp.elementWiseMultiply (one, *D_, Y_tmp, zero);
1090 L_solver_->apply (Y_tmp, Y, mode);
1093 U_solver_->apply (X, Y, mode);
1095 if (!this->isKokkosKernelsSpiluk_) {
1101 Y.elementWiseMultiply (one, *D_, Y, zero);
1104 L_solver_->apply (Y, Y, mode);
1109 if (alpha == zero) {
1119 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1120 apply (X, Y_tmp, mode);
1121 Y.update (alpha, Y_tmp, beta);
1126#ifdef HAVE_IFPACK2_DEBUG
1128 Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
1131 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
1132 if (STM::isnaninf (norms[j])) {
1137 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
1142 applyTime_ += (timer.wallTime() - startTime);
1179template<
class MatrixType>
1182 std::ostringstream os;
1187 os <<
"\"Ifpack2::RILUK\": {";
1188 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
1189 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
1191 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
1193 if (A_.is_null ()) {
1194 os <<
"Matrix: null";
1197 os <<
"Global matrix dimensions: ["
1198 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
1199 <<
", Global nnz: " << A_->getGlobalNumEntries();
1202 if (! L_solver_.is_null ()) os <<
", " << L_solver_->description ();
1203 if (! U_solver_.is_null ()) os <<
", " << U_solver_->description ();
1211#define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \
1212 template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:100
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:163
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:88
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:254
crs_matrix_type::impl_scalar_type impl_scalar_type
Scalar type stored in Kokkos::Views (CrsMatrix and MultiVector)
Definition: Ifpack2_RILUK_decl.hpp:293
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_RILUK_def.hpp:455
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:290
virtual ~RILUK()
Destructor (declared virtual for memory safety).
Definition: Ifpack2_RILUK_def.hpp:121
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:989
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_RILUK_decl.hpp:269
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:256
const crs_matrix_type & getU() const
Return the U factor of the ILU factorization.
Definition: Ifpack2_RILUK_def.hpp:207
std::string description() const
A one-line description of this object.
Definition: Ifpack2_RILUK_def.hpp:1180
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_RILUK_def.hpp:140
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:257
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:193
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:266
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the input matrix.
Definition: Ifpack2_RILUK_def.hpp:408
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:415
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:237
void setParameters(const Teuchos::ParameterList ¶ms)
Definition: Ifpack2_RILUK_def.hpp:303
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:260
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:263
const crs_matrix_type & getL() const
Return the L factor of the ILU factorization.
Definition: Ifpack2_RILUK_def.hpp:176
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_RILUK_def.hpp:220
void compute()
Compute the (numeric) incomplete factorization.
Definition: Ifpack2_RILUK_def.hpp:728
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74