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 "Kokkos_Sort.hpp" 50 #include "KokkosKernels_SparseUtils.hpp" 51 #include "KokkosKernels_Sorting.hpp" 62 static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
64 type_strs[0] =
"Serial";
65 type_strs[1] =
"KSPILUK";
67 type_enums[0] = Serial;
68 type_enums[1] = KSPILUK;
73 template<
class MatrixType>
79 isInitialized_ (false),
84 initializeTime_ (0.0),
90 isKokkosKernelsSpiluk_(false)
96 template<
class MatrixType>
101 isAllocated_ (false),
102 isInitialized_ (false),
107 initializeTime_ (0.0),
113 isKokkosKernelsSpiluk_(false)
119 template<
class MatrixType>
122 if (Teuchos::nonnull (KernelHandle_))
124 KernelHandle_->destroy_spiluk_handle();
128 template<
class MatrixType>
135 template<
class MatrixType>
143 if (A.getRawPtr () != A_.getRawPtr ()) {
144 isAllocated_ =
false;
145 isInitialized_ =
false;
147 A_local_ = Teuchos::null;
148 Graph_ = Teuchos::null;
155 if (! L_solver_.is_null ()) {
158 if (! U_solver_.is_null ()) {
159 U_solver_->setMatrix (Teuchos::null);
171 template<
class MatrixType>
175 TEUCHOS_TEST_FOR_EXCEPTION(
176 L_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor " 177 "is null. Please call initialize() (and preferably also compute()) " 178 "before calling this method. If the input matrix has not yet been set, " 179 "you must first call setMatrix() with a nonnull input matrix before you " 180 "may call initialize() or compute().");
185 template<
class MatrixType>
186 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
192 TEUCHOS_TEST_FOR_EXCEPTION(
193 D_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor " 194 "(of diagonal entries) is null. Please call initialize() (and " 195 "preferably also compute()) before calling this method. If the input " 196 "matrix has not yet been set, you must first call setMatrix() with a " 197 "nonnull input matrix before you may call initialize() or compute().");
202 template<
class MatrixType>
206 TEUCHOS_TEST_FOR_EXCEPTION(
207 U_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U factor " 208 "is null. Please call initialize() (and preferably also compute()) " 209 "before calling this method. If the input matrix has not yet been set, " 210 "you must first call setMatrix() with a nonnull input matrix before you " 211 "may call initialize() or compute().");
216 template<
class MatrixType>
218 TEUCHOS_TEST_FOR_EXCEPTION(
219 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getNodeSmootherComplexity: " 220 "The input matrix A is null. Please call setMatrix() with a nonnull " 221 "input matrix, then call compute(), before calling this method.");
223 if(!L_.is_null() && !U_.is_null())
224 return A_->getNodeNumEntries() + L_->getNodeNumEntries() + U_->getNodeNumEntries();
230 template<
class MatrixType>
231 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: " 237 "The matrix is null. Please call setMatrix() with a nonnull input " 238 "before calling this method.");
241 TEUCHOS_TEST_FOR_EXCEPTION(
242 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: " 243 "The computed graph is null. Please call initialize() before calling " 249 template<
class MatrixType>
250 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
254 TEUCHOS_TEST_FOR_EXCEPTION(
255 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: " 256 "The matrix is null. Please call setMatrix() with a nonnull input " 257 "before calling this method.");
260 TEUCHOS_TEST_FOR_EXCEPTION(
261 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: " 262 "The computed graph is null. Please call initialize() before calling " 268 template<
class MatrixType>
274 if (! isAllocated_) {
283 L_ = rcp (
new crs_matrix_type (Graph_->getL_Graph ()));
284 U_ = rcp (
new crs_matrix_type (Graph_->getU_Graph ()));
285 L_->setAllToScalar (STS::zero ());
286 U_->setAllToScalar (STS::zero ());
291 D_ = rcp (
new vec_type (Graph_->getL_Graph ()->getRowMap ()));
297 template<
class MatrixType>
300 setParameters (
const Teuchos::ParameterList& params)
303 using Teuchos::ParameterList;
304 using Teuchos::Array;
305 using Details::getParamTryingTypes;
306 const char prefix[] =
"Ifpack2::RILUK: ";
313 double overalloc = 2.;
324 const std::string paramName (
"fact: iluk level-of-fill");
325 getParamTryingTypes<int, int, global_ordinal_type, double, float>
326 (fillLevel, params, paramName, prefix);
331 const std::string paramName (
"fact: absolute threshold");
332 getParamTryingTypes<magnitude_type, magnitude_type, double>
333 (absThresh, params, paramName, prefix);
336 const std::string paramName (
"fact: relative threshold");
337 getParamTryingTypes<magnitude_type, magnitude_type, double>
338 (relThresh, params, paramName, prefix);
341 const std::string paramName (
"fact: relax value");
342 getParamTryingTypes<magnitude_type, magnitude_type, double>
343 (relaxValue, params, paramName, prefix);
346 const std::string paramName (
"fact: iluk overalloc");
347 getParamTryingTypes<double, double>
348 (overalloc, params, paramName, prefix);
352 Details::IlukImplType::Enum ilukimplType = Details::IlukImplType::Serial;
354 static const char typeName[] =
"fact: type";
356 if ( ! params.isType<std::string>(typeName))
break;
359 Array<std::string> ilukimplTypeStrs;
360 Array<Details::IlukImplType::Enum> ilukimplTypeEnums;
361 Details::IlukImplType::loadPLTypeOption (ilukimplTypeStrs, ilukimplTypeEnums);
362 Teuchos::StringToIntegralParameterEntryValidator<Details::IlukImplType::Enum>
363 s2i(ilukimplTypeStrs (), ilukimplTypeEnums (), typeName,
false);
365 ilukimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
368 if (ilukimplType == Details::IlukImplType::KSPILUK) {
369 this->isKokkosKernelsSpiluk_ =
true;
372 this->isKokkosKernelsSpiluk_ =
false;
376 L_solver_->setParameters(params);
377 U_solver_->setParameters(params);
383 LevelOfFill_ = fillLevel;
384 Overalloc_ = overalloc;
391 if (absThresh != -STM::one ()) {
392 Athresh_ = absThresh;
394 if (relThresh != -STM::one ()) {
395 Rthresh_ = relThresh;
397 if (relaxValue != -STM::one ()) {
398 RelaxValue_ = relaxValue;
403 template<
class MatrixType>
404 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
410 template<
class MatrixType>
411 Teuchos::RCP<const typename RILUK<MatrixType>::crs_matrix_type>
417 template<
class MatrixType>
418 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
423 using Teuchos::rcp_dynamic_cast;
424 using Teuchos::rcp_implicit_cast;
429 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
430 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
437 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
438 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
439 if (! A_lf_r.is_null ()) {
440 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
446 return rcp (
new LocalFilter<row_matrix_type> (A));
451 template<
class MatrixType>
456 using Teuchos::rcp_const_cast;
457 using Teuchos::rcp_dynamic_cast;
458 using Teuchos::rcp_implicit_cast;
459 using Teuchos::Array;
460 using Teuchos::ArrayView;
464 const char prefix[] =
"Ifpack2::RILUK::initialize: ";
466 TEUCHOS_TEST_FOR_EXCEPTION
467 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please " 468 "call setMatrix() with a nonnull input before calling this method.");
469 TEUCHOS_TEST_FOR_EXCEPTION
470 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not " 471 "fill complete. You may not invoke initialize() or compute() with this " 472 "matrix until the matrix is fill complete. If your matrix is a " 473 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and " 474 "range Maps, if appropriate) before calling this method.");
476 Teuchos::Time timer (
"RILUK::initialize");
477 double startTime = timer.wallTime();
479 Teuchos::TimeMonitor timeMon (timer);
488 isInitialized_ =
false;
489 isAllocated_ =
false;
491 Graph_ = Teuchos::null;
493 A_local_ = makeLocalFilter (A_);
494 TEUCHOS_TEST_FOR_EXCEPTION(
495 A_local_.is_null (), std::logic_error,
"Ifpack2::RILUK::initialize: " 496 "makeLocalFilter returned null; it failed to compute A_local. " 497 "Please report this bug to the Ifpack2 developers.");
505 if (this->isKokkosKernelsSpiluk_) {
506 this->KernelHandle_ = Teuchos::rcp (
new kk_handle_type ());
507 KernelHandle_->create_spiluk_handle( KokkosSparse::Experimental::SPILUKAlgorithm::SEQLVLSCHD_TP1,
508 A_local_->getNodeNumRows(),
509 2*A_local_->getNodeNumEntries()*(LevelOfFill_+1),
510 2*A_local_->getNodeNumEntries()*(LevelOfFill_+1) );
514 RCP<const crs_matrix_type> A_local_crs =
516 if (A_local_crs.is_null ()) {
518 Array<size_t> entriesPerRow(numRows);
520 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
522 RCP<crs_matrix_type> A_local_crs_nc =
524 A_local_->getColMap (),
527 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getNodeMaxNumRowEntries());
528 nonconst_values_host_view_type values(
"values",A_local_->getNodeMaxNumRowEntries());
530 size_t numEntries = 0;
531 A_local_->getLocalRowCopy(i, indices, values, numEntries);
532 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
534 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
538 LevelOfFill_, 0, Overalloc_));
541 if (this->isKokkosKernelsSpiluk_) Graph_->initialize (KernelHandle_);
542 else Graph_->initialize ();
545 checkOrderingConsistency (*A_local_);
546 L_solver_->setMatrix (L_);
547 L_solver_->initialize ();
548 U_solver_->setMatrix (U_);
549 U_solver_->initialize ();
556 isInitialized_ =
true;
558 initializeTime_ += (timer.wallTime() - startTime);
561 template<
class MatrixType>
569 Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getNodeElementList();
570 Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getNodeElementList();
571 bool gidsAreConsistentlyOrdered=
true;
572 global_ordinal_type indexOfInconsistentGID=0;
573 for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
574 if (rowGIDs[i] != colGIDs[i]) {
575 gidsAreConsistentlyOrdered=
false;
576 indexOfInconsistentGID=i;
580 TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==
false, std::runtime_error,
581 "The ordering of the local GIDs in the row and column maps is not the same" 582 << std::endl <<
"at index " << indexOfInconsistentGID
583 <<
". Consistency is required, as all calculations are done with" 584 << std::endl <<
"local indexing.");
587 template<
class MatrixType>
590 initAllValues (
const row_matrix_type& A)
592 using Teuchos::ArrayRCP;
596 using Teuchos::REDUCE_SUM;
597 using Teuchos::reduceAll;
598 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
600 size_t NumIn = 0, NumL = 0, NumU = 0;
601 bool DiagFound =
false;
602 size_t NumNonzeroDiags = 0;
603 size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
607 nonconst_local_inds_host_view_type InI(
"InI",MaxNumEntries);
608 Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
609 Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
610 nonconst_values_host_view_type InV(
"InV",MaxNumEntries);
611 Teuchos::Array<scalar_type> LV(MaxNumEntries);
612 Teuchos::Array<scalar_type> UV(MaxNumEntries);
615 const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
620 L_->setAllToScalar (STS::zero ());
621 U_->setAllToScalar (STS::zero ());
624 D_->putScalar (STS::zero ());
625 auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
627 RCP<const map_type> rowMap = L_->getRowMap ();
637 Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getNodeElementList();
638 for (
size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
639 local_ordinal_type local_row = myRow;
643 A.getLocalRowCopy (local_row, InI, InV, NumIn);
651 for (
size_t j = 0; j < NumIn; ++j) {
652 const local_ordinal_type k = InI[j];
654 if (k == local_row) {
657 DV(local_row) += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
660 TEUCHOS_TEST_FOR_EXCEPTION(
661 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current " 662 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is " 663 "probably an artifact of the undocumented assumptions of the " 664 "original implementation (likely copied and pasted from Ifpack). " 665 "Nevertheless, the code I found here insisted on this being an error " 666 "state, so I will throw an exception here.");
668 else if (k < local_row) {
673 else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
685 DV(local_row) = Athresh_;
690 L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
694 L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
700 U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
704 U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
712 isInitialized_ =
true;
716 template<
class MatrixType>
721 using Teuchos::rcp_const_cast;
722 using Teuchos::rcp_dynamic_cast;
723 using Teuchos::Array;
724 using Teuchos::ArrayView;
725 const char prefix[] =
"Ifpack2::RILUK::compute: ";
730 TEUCHOS_TEST_FOR_EXCEPTION
731 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please " 732 "call setMatrix() with a nonnull input before calling this method.");
733 TEUCHOS_TEST_FOR_EXCEPTION
734 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not " 735 "fill complete. You may not invoke initialize() or compute() with this " 736 "matrix until the matrix is fill complete. If your matrix is a " 737 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and " 738 "range Maps, if appropriate) before calling this method.");
740 if (! isInitialized ()) {
744 Teuchos::Time timer (
"RILUK::compute");
747 Teuchos::TimeMonitor timeMon (timer);
748 double startTime = timer.wallTime();
752 if (!this->isKokkosKernelsSpiluk_) {
755 initAllValues (*A_local_);
761 const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
763 size_t NumIn, NumL, NumU;
766 const size_t MaxNumEntries =
767 L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
769 Teuchos::Array<local_ordinal_type> InI(MaxNumEntries);
770 Teuchos::Array<scalar_type> InV(MaxNumEntries);
771 size_t num_cols = U_->getColMap()->getNodeNumElements();
772 Teuchos::Array<int> colflag(num_cols);
774 auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
778 for (
size_t j = 0; j < num_cols; ++j) {
781 using IST =
typename row_matrix_type::impl_scalar_type;
782 for (
size_t i = 0; i < L_->getNodeNumRows (); ++i) {
786 local_inds_host_view_type UUI;
787 values_host_view_type UUV;
791 NumIn = MaxNumEntries;
792 nonconst_local_inds_host_view_type InI_v(InI.data(),MaxNumEntries);
793 nonconst_values_host_view_type InV_v(reinterpret_cast<IST*>(InV.data()),MaxNumEntries);
795 L_->getLocalRowCopy (local_row, InI_v , InV_v, NumL);
798 InI[NumL] = local_row;
800 nonconst_local_inds_host_view_type InI_sub(InI.data()+NumL+1,MaxNumEntries-NumL-1);
801 nonconst_values_host_view_type InV_sub(reinterpret_cast<IST*>(InV.data())+NumL+1,MaxNumEntries-NumL-1);
803 U_->getLocalRowCopy (local_row, InI_sub,InV_sub, NumU);
807 for (
size_t j = 0; j < NumIn; ++j) {
813 for (
size_t jj = 0; jj < NumL; ++jj) {
815 IST multiplier = InV[jj];
819 U_->getLocalRowView(j, UUI, UUV);
822 if (RelaxValue_ == STM::zero ()) {
823 for (
size_t k = 0; k < NumUU; ++k) {
824 const int kk = colflag[UUI[k]];
829 InV[kk] -=
static_cast<scalar_type>(multiplier * UUV[k]);
835 for (
size_t k = 0; k < NumUU; ++k) {
839 const int kk = colflag[UUI[k]];
841 InV[kk] -=
static_cast<scalar_type>(multiplier*UUV[k]);
844 diagmod -=
static_cast<scalar_type>(multiplier*UUV[k]);
852 L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
857 if (RelaxValue_ != STM::zero ()) {
858 DV(i) += RelaxValue_*diagmod;
861 if (STS::magnitude (DV(i)) > STS::magnitude (MaxDiagonalValue)) {
862 if (STS::real (DV(i)) < STM::zero ()) {
863 DV(i) = -MinDiagonalValue;
866 DV(i) = MinDiagonalValue;
873 for (
size_t j = 0; j < NumU; ++j) {
879 U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
883 for (
size_t j = 0; j < NumIn; ++j) {
884 colflag[InI[j]] = -1;
895 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
896 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
899 L_solver_->setMatrix (L_);
900 L_solver_->compute ();
901 U_solver_->setMatrix (U_);
902 U_solver_->compute ();
906 RCP<const crs_matrix_type> A_local_crs =
908 if (A_local_crs.is_null ()) {
910 Array<size_t> entriesPerRow(numRows);
912 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
914 RCP<crs_matrix_type> A_local_crs_nc =
916 A_local_->getColMap (),
919 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getNodeMaxNumRowEntries());
920 nonconst_values_host_view_type values(
"values",A_local_->getNodeMaxNumRowEntries());
922 size_t numEntries = 0;
923 A_local_->getLocalRowCopy(i, indices, values, numEntries);
924 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()),indices.data());
926 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
929 A_local_rowmap_ = A_local_crs->getLocalMatrixDevice().graph.row_map;
930 A_local_entries_ = A_local_crs->getLocalMatrixDevice().graph.entries;
931 A_local_values_ = A_local_crs->getLocalValuesView();
937 if (L_->isStaticGraph () || L_->isLocallyIndexed ()) {
938 L_->setAllToScalar (STS::zero ());
939 U_->setAllToScalar (STS::zero ());
942 using row_map_type =
typename crs_matrix_type::local_matrix_device_type::row_map_type;
944 row_map_type L_rowmap = L_->getLocalMatrixDevice().graph.row_map;
945 auto L_entries = L_->getLocalMatrixDevice().graph.entries;
946 auto L_values = L_->getLocalValuesView();
947 row_map_type U_rowmap = U_->getLocalMatrixDevice().graph.row_map;
948 auto U_entries = U_->getLocalMatrixDevice().graph.entries;
949 auto U_values = U_->getLocalValuesView();
951 KokkosSparse::Experimental::spiluk_numeric( KernelHandle_.getRawPtr(), LevelOfFill_,
952 A_local_rowmap_, A_local_entries_, A_local_values_,
953 L_rowmap, L_entries, L_values, U_rowmap, U_entries, U_values );
955 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
956 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
958 L_solver_->setMatrix (L_);
959 L_solver_->compute ();
960 U_solver_->setMatrix (U_);
961 U_solver_->compute ();
966 computeTime_ += (timer.wallTime() - startTime);
970 template<
class MatrixType>
973 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
974 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
975 Teuchos::ETransp mode,
980 using Teuchos::rcpFromRef;
982 TEUCHOS_TEST_FOR_EXCEPTION(
983 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::apply: The matrix is " 984 "null. Please call setMatrix() with a nonnull input, then initialize() " 985 "and compute(), before calling this method.");
986 TEUCHOS_TEST_FOR_EXCEPTION(
987 ! isComputed (), std::runtime_error,
988 "Ifpack2::RILUK::apply: If you have not yet called compute(), " 989 "you must call compute() before calling this method.");
990 TEUCHOS_TEST_FOR_EXCEPTION(
991 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
992 "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. " 993 "X.getNumVectors() = " << X.getNumVectors ()
994 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
995 TEUCHOS_TEST_FOR_EXCEPTION(
996 STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
997 "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for " 998 "complex Scalar type. Please talk to the Ifpack2 developers to get this " 999 "fixed. There is a FIXME in this file about this very issue.");
1000 #ifdef HAVE_IFPACK2_DEBUG 1003 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.");
1004 Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
1007 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
1008 if (STM::isnaninf (norms[j])) {
1013 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
1015 #endif // HAVE_IFPACK2_DEBUG 1020 Teuchos::Time timer (
"RILUK::apply");
1021 double startTime = timer.wallTime();
1023 Teuchos::TimeMonitor timeMon (timer);
1024 if (alpha == one && beta == zero) {
1025 if (mode == Teuchos::NO_TRANS) {
1027 L_solver_->apply (X, Y, mode);
1029 if (!this->isKokkosKernelsSpiluk_) {
1032 Y.elementWiseMultiply (one, *D_, Y, zero);
1035 U_solver_->apply (Y, Y, mode);
1039 U_solver_->apply (X, Y, mode);
1041 if (!this->isKokkosKernelsSpiluk_) {
1047 Y.elementWiseMultiply (one, *D_, Y, zero);
1050 L_solver_->apply (Y, Y, mode);
1054 if (alpha == zero) {
1064 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1065 apply (X, Y_tmp, mode);
1066 Y.update (alpha, Y_tmp, beta);
1071 #ifdef HAVE_IFPACK2_DEBUG 1073 Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
1076 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
1077 if (STM::isnaninf (norms[j])) {
1082 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
1084 #endif // HAVE_IFPACK2_DEBUG 1087 applyTime_ += (timer.wallTime() - startTime);
1124 template<
class MatrixType>
1127 std::ostringstream os;
1132 os <<
"\"Ifpack2::RILUK\": {";
1133 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 1134 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
1136 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
1138 if (A_.is_null ()) {
1139 os <<
"Matrix: null";
1142 os <<
"Global matrix dimensions: [" 1143 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]" 1147 if (! L_solver_.is_null ()) os <<
", " << L_solver_->description ();
1148 if (! U_solver_.is_null ()) os <<
", " << U_solver_->description ();
1156 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \ 1157 template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >; Ifpack2::Preconditioner< Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::scalar_type, Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::local_ordinal_type, Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::global_ordinal_type, Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::node_type >::magnitude_type Teuchos::ScalarTraits< Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Preconditioner.hpp:111
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:263
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:281
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_RILUK_def.hpp:137
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:245
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:284
Ifpack2 implementation details.
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:266
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:260
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:234
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:79
crs_matrix_type::impl_scalar_type impl_scalar_type
Scalar type stored in Kokkos::Views (CrsMatrix and MultiVector)
Definition: Ifpack2_RILUK_decl.hpp:293
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:100
Tpetra::global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_RILUK_decl.hpp:546
Definition: Ifpack2_Container_decl.hpp:576
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
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:253
Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:257