43 #ifndef IFPACK2_CRSRILUK_DEF_HPP 44 #define IFPACK2_CRSRILUK_DEF_HPP 46 #include "Ifpack2_LocalFilter.hpp" 47 #include "Tpetra_CrsMatrix.hpp" 48 #include "Ifpack2_LocalSparseTriangularSolver.hpp" 50 #ifdef IFPACK2_ILUK_EXPERIMENTAL 51 #include <shylubasker_def.hpp> 52 #include <Kokkos_CrsMatrix.hpp> 53 # ifdef IFPACK2_HTS_EXPERIMENTAL 54 # include <shylu_hts.hpp> 60 template<
class MatrixType>
65 isInitialized_ (false),
67 isExperimental_ (false),
71 initializeTime_ (0.0),
80 template<
class MatrixType>
85 isInitialized_ (false),
87 isExperimental_ (false),
91 initializeTime_ (0.0),
100 template<
class MatrixType>
104 template<
class MatrixType>
112 if (A.getRawPtr () != A_.getRawPtr ()) {
113 isAllocated_ =
false;
114 isInitialized_ =
false;
116 A_local_ = Teuchos::null;
117 Graph_ = Teuchos::null;
124 if (! L_solver_.is_null ()) {
127 if (! U_solver_.is_null ()) {
128 U_solver_->setMatrix (Teuchos::null);
140 template<
class MatrixType>
144 TEUCHOS_TEST_FOR_EXCEPTION(
145 L_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getL: The L factor " 146 "is null. Please call initialize() (and preferably also compute()) " 147 "before calling this method. If the input matrix has not yet been set, " 148 "you must first call setMatrix() with a nonnull input matrix before you " 149 "may call initialize() or compute().");
154 template<
class MatrixType>
155 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
161 TEUCHOS_TEST_FOR_EXCEPTION(
162 D_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getD: The D factor " 163 "(of diagonal entries) is null. Please call initialize() (and " 164 "preferably also compute()) before calling this method. If the input " 165 "matrix has not yet been set, you must first call setMatrix() with a " 166 "nonnull input matrix before you may call initialize() or compute().");
171 template<
class MatrixType>
175 TEUCHOS_TEST_FOR_EXCEPTION(
176 U_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getU: The U 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 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
190 TEUCHOS_TEST_FOR_EXCEPTION(
191 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: " 192 "The matrix is null. Please call setMatrix() with a nonnull input " 193 "before calling this method.");
196 TEUCHOS_TEST_FOR_EXCEPTION(
197 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getDomainMap: " 198 "The computed graph is null. Please call initialize() before calling " 204 template<
class MatrixType>
205 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
209 TEUCHOS_TEST_FOR_EXCEPTION(
210 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: " 211 "The matrix is null. Please call setMatrix() with a nonnull input " 212 "before calling this method.");
215 TEUCHOS_TEST_FOR_EXCEPTION(
216 Graph_.is_null (), std::runtime_error,
"Ifpack2::RILUK::getRangeMap: " 217 "The computed graph is null. Please call initialize() before calling " 223 template<
class MatrixType>
229 if (! isAllocated_) {
238 L_ = rcp (
new crs_matrix_type (Graph_->getL_Graph ()));
239 U_ = rcp (
new crs_matrix_type (Graph_->getU_Graph ()));
240 L_->setAllToScalar (STS::zero ());
241 U_->setAllToScalar (STS::zero ());
246 D_ = rcp (
new vec_type (Graph_->getL_Graph ()->getRowMap ()));
252 template<
class MatrixType>
255 setParameters (
const Teuchos::ParameterList& params)
258 using Teuchos::Exceptions::InvalidParameterName;
259 using Teuchos::Exceptions::InvalidParameterType;
275 bool gotFillLevel =
false;
277 fillLevel = params.get<
int> (
"fact: iluk level-of-fill");
280 catch (InvalidParameterType&) {
285 catch (InvalidParameterName&) {
289 if (! gotFillLevel) {
295 catch (InvalidParameterType&) {
301 if (! gotFillLevel) {
305 fillLevel = as<int> (params.get<
magnitude_type> (
"fact: iluk level-of-fill"));
308 catch (InvalidParameterType&) {
314 if (! gotFillLevel) {
318 fillLevel = as<int> (params.get<
double> (
"fact: iluk level-of-fill"));
321 catch (InvalidParameterType& e) {
331 TEUCHOS_TEST_FOR_EXCEPTION(
334 "Ifpack2::RILUK::setParameters: We should never get here! " 335 "The method should either have read the \"fact: iluk level-of-fill\" " 336 "parameter by this point, or have thrown an exception. " 337 "Please let the Ifpack2 developers know about this bug.");
345 absThresh = params.get<
magnitude_type> (
"fact: absolute threshold");
347 catch (InvalidParameterType&) {
350 absThresh = as<magnitude_type> (params.get<
double> (
"fact: absolute threshold"));
352 catch (InvalidParameterName&) {
357 relThresh = params.get<
magnitude_type> (
"fact: relative threshold");
359 catch (InvalidParameterType&) {
362 relThresh = as<magnitude_type> (params.get<
double> (
"fact: relative threshold"));
364 catch (InvalidParameterName&) {
371 catch (InvalidParameterType&) {
374 relaxValue = as<magnitude_type> (params.get<
double> (
"fact: relax value"));
376 catch (InvalidParameterName&) {
383 #ifdef IFPACK2_ILUK_EXPERIMENTAL 385 isExperimental_ = params.get<
bool> (
"fact: iluk experimental basker");
387 catch (InvalidParameterType&) {
390 catch (InvalidParameterName&) {
395 basker_threads = params.get<
int> (
"fact: iluk experimental basker threads");
397 catch (InvalidParameterType&) {
400 catch (InvalidParameterName&) {
405 basker_user_fill = (
scalar_type) params.get<
double>(
"fact: iluk experimental basker user_fill");
407 catch (InvalidParameterType&) {
410 catch (InvalidParameterName&) {
414 # ifdef IFPACK2_HTS_EXPERIMENTAL 416 hts_nthreads_ = basker_threads;
418 const std::string name(
"fact: iluk experimental HTS");
419 if (params.isType<
bool> (name))
420 use_hts_ = params.get<
bool> (name);
423 const std::string name(
"fact: iluk experimental HTS threads");
424 if (params.isType<
int> (name))
425 hts_nthreads_ = params.get<
int> (name);
434 LevelOfFill_ = fillLevel;
441 if (absThresh != -STM::one ()) {
442 Athresh_ = absThresh;
444 if (relThresh != -STM::one ()) {
445 Rthresh_ = relThresh;
447 if (relaxValue != -STM::one ()) {
448 RelaxValue_ = relaxValue;
453 template<
class MatrixType>
454 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
460 template<
class MatrixType>
461 Teuchos::RCP<const typename RILUK<MatrixType>::crs_matrix_type>
467 template<
class MatrixType>
468 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
473 using Teuchos::rcp_dynamic_cast;
474 using Teuchos::rcp_implicit_cast;
479 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
480 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
487 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
488 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
489 if (! A_lf_r.is_null ()) {
490 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
496 return rcp (
new LocalFilter<row_matrix_type> (A));
501 template<
class MatrixType>
506 using Teuchos::rcp_const_cast;
507 using Teuchos::rcp_dynamic_cast;
508 using Teuchos::rcp_implicit_cast;
512 const char prefix[] =
"Ifpack2::RILUK::initialize: ";
514 TEUCHOS_TEST_FOR_EXCEPTION
515 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please " 516 "call setMatrix() with a nonnull input before calling this method.");
517 TEUCHOS_TEST_FOR_EXCEPTION
518 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not " 519 "fill complete. You may not invoke initialize() or compute() with this " 520 "matrix until the matrix is fill complete. If your matrix is a " 521 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and " 522 "range Maps, if appropriate) before calling this method.");
524 Teuchos::Time timer (
"RILUK::initialize");
526 Teuchos::TimeMonitor timeMon (timer);
535 isInitialized_ =
false;
536 isAllocated_ =
false;
538 Graph_ = Teuchos::null;
540 A_local_ = makeLocalFilter (A_);
541 TEUCHOS_TEST_FOR_EXCEPTION(
542 A_local_.is_null (), std::logic_error,
"Ifpack2::RILUK::initialize: " 543 "makeLocalFilter returned null; it failed to compute A_local. " 544 "Please report this bug to the Ifpack2 developers.");
552 RCP<const crs_matrix_type> A_local_crs =
554 if (A_local_crs.is_null ()) {
559 RCP<crs_matrix_type> A_local_crs_nc =
561 A_local_->getColMap (), 0));
570 import_type
import (A_local_->getRowMap (), A_local_->getRowMap ());
571 A_local_crs_nc->doImport (*A_local_,
import, Tpetra::REPLACE);
572 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
577 #ifdef IFPACK2_ILUK_EXPERIMENTAL 578 if (isExperimental_) A_local_crs_ = A_local_crs;
584 Graph_->initialize ();
586 checkOrderingConsistency (*A_local_);
593 #ifdef IFPACK2_ILUK_EXPERIMENTAL 594 typedef typename node_type::device_type kokkos_device;
595 typedef typename kokkos_device::execution_space kokkos_exe;
596 typedef typename crs_matrix_type::local_matrix_type kokkos_csr_matrix;
598 static_assert( std::is_same< kokkos_exe,
599 Kokkos::OpenMP>::value,
600 "Kokkos node type not supported by experimental thread basker RILUK");
603 myBasker = rcp(
new BaskerNS::Basker<local_ordinal_type, scalar_type, Kokkos::OpenMP>);
604 myBasker->Options.no_pivot =
true;
605 myBasker->Options.transpose =
true;
606 myBasker->Options.symmetric =
false;
607 myBasker->Options.realloc =
true;
608 myBasker->Options.btf =
false;
609 myBasker->Options.incomplete =
true;
610 myBasker->Options.inc_lvl = LevelOfFill_;
611 myBasker->Options.user_fill = basker_user_fill;
612 myBasker->Options.same_pattern =
false;
613 myBasker->SetThreads(basker_threads);
614 basker_reuse =
false;
616 kokkos_csr_matrix kcsr = A_local_crs_->getLocalMatrix();
624 bi < A_local_crs_->getNodeNumRows()+1; bi++)
635 &(kcsr.graph.entries(0)),
638 TEUCHOS_TEST_FOR_EXCEPTION(
639 basker_error != 0, std::logic_error,
"Ifpack2::RILUK::initialize:" 640 "Error in basker Symbolic");
644 TEUCHOS_TEST_FOR_EXCEPTION(
645 0==0, std::logic_error,
"Ifpack2::RILUK::initialize: " 646 "Using experimental ILUK without compiling experimental " 647 "Try again with -DIFPACK2_ILUK_EXPERIMENAL.");
653 isInitialized_ =
true;
655 initializeTime_ += timer.totalElapsedTime ();
658 #ifdef IFPACK2_ILUK_EXPERIMENTAL 660 template<
class MatrixType>
667 using Teuchos::rcp_dynamic_cast;
668 using Teuchos::rcp_implicit_cast;
669 using Teuchos::rcp_const_cast;
671 RCP<crs_matrix_type> A_local_crs_nc = rcp_const_cast<crs_matrix_type> (A_local_crs_);
672 A_local_crs_ = rcp_dynamic_cast<
const crs_matrix_type> (A_local_);
673 if (A_local_crs_.is_null ()) {
674 Teuchos::Array<local_ordinal_type> lclColInds;
675 Teuchos::Array<scalar_type> vals;
676 A_local_crs_nc->resumeFill();
677 for (
size_t i = 0; i < A_local_crs_nc->getNodeNumRows(); ++i) {
679 const auto nc = A_local_->getNumEntriesInLocalRow(i);
680 if (static_cast<size_t>(lclColInds.size()) < nc) {
681 lclColInds.resize(nc);
684 A_local_->getLocalRowCopy(i, lclColInds(), vals(), numEnt);
685 A_local_crs_nc->replaceLocalValues(i, lclColInds.view(0, numEnt), vals.view(0, numEnt));
687 A_local_crs_nc->fillComplete();
689 A_local_crs_ = rcp_const_cast<
const crs_matrix_type> (A_local_crs_nc);
693 template<
class MatrixType>
696 checkOrderingConsistency (
const row_matrix_type& A)
701 Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getNodeElementList();
702 Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getNodeElementList();
703 bool gidsAreConsistentlyOrdered=
true;
704 global_ordinal_type indexOfInconsistentGID=0;
705 for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
706 if (rowGIDs[i] != colGIDs[i]) {
707 gidsAreConsistentlyOrdered=
false;
708 indexOfInconsistentGID=i;
712 TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==
false, std::runtime_error,
713 "The ordering of the local GIDs in the row and column maps is not the same" 714 << std::endl <<
"at index " << indexOfInconsistentGID
715 <<
". Consistency is required, as all calculations are done with" 716 << std::endl <<
"local indexing.");
719 template<
class MatrixType>
722 initAllValues (
const row_matrix_type& A)
724 using Teuchos::ArrayRCP;
728 using Teuchos::REDUCE_SUM;
729 using Teuchos::reduceAll;
730 typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
732 size_t NumIn = 0, NumL = 0, NumU = 0;
733 bool DiagFound =
false;
734 size_t NumNonzeroDiags = 0;
735 size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
739 Teuchos::Array<local_ordinal_type> InI(MaxNumEntries);
740 Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
741 Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
742 Teuchos::Array<scalar_type> InV(MaxNumEntries);
743 Teuchos::Array<scalar_type> LV(MaxNumEntries);
744 Teuchos::Array<scalar_type> UV(MaxNumEntries);
747 const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
752 L_->setAllToScalar (STS::zero ());
753 U_->setAllToScalar (STS::zero ());
756 D_->putScalar (STS::zero ());
757 ArrayRCP<scalar_type> DV = D_->get1dViewNonConst ();
759 RCP<const map_type> rowMap = L_->getRowMap ();
769 Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getNodeElementList();
770 for (
size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
771 local_ordinal_type local_row = myRow;
775 A.getLocalRowCopy (local_row, InI(), InV(), NumIn);
783 for (
size_t j = 0; j < NumIn; ++j) {
784 const local_ordinal_type k = InI[j];
786 if (k == local_row) {
789 DV[local_row] += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
792 TEUCHOS_TEST_FOR_EXCEPTION(
793 true, std::runtime_error,
"Ifpack2::RILUK::initAllValues: current " 794 "GID k = " << k <<
" < 0. I'm not sure why this is an error; it is " 795 "probably an artifact of the undocumented assumptions of the " 796 "original implementation (likely copied and pasted from Ifpack). " 797 "Nevertheless, the code I found here insisted on this being an error " 798 "state, so I will throw an exception here.");
800 else if (k < local_row) {
805 else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
817 DV[local_row] = Athresh_;
822 L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
826 L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
832 U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
836 U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
844 isInitialized_ =
true;
848 template<
class MatrixType>
852 const char prefix[] =
"Ifpack2::RILUK::compute: ";
857 TEUCHOS_TEST_FOR_EXCEPTION
858 (A_.is_null (), std::runtime_error, prefix <<
"The matrix is null. Please " 859 "call setMatrix() with a nonnull input before calling this method.");
860 TEUCHOS_TEST_FOR_EXCEPTION
861 (! A_->isFillComplete (), std::runtime_error, prefix <<
"The matrix is not " 862 "fill complete. You may not invoke initialize() or compute() with this " 863 "matrix until the matrix is fill complete. If your matrix is a " 864 "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and " 865 "range Maps, if appropriate) before calling this method.");
867 if (! isInitialized ()) {
871 Teuchos::Time timer (
"RILUK::compute");
872 if ( ! isExperimental_) {
874 Teuchos::TimeMonitor timeMon (timer);
880 initAllValues (*A_local_);
886 const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
888 size_t NumIn, NumL, NumU;
891 const size_t MaxNumEntries =
892 L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
894 Teuchos::Array<local_ordinal_type> InI(MaxNumEntries);
895 Teuchos::Array<scalar_type> InV(MaxNumEntries);
896 size_t num_cols = U_->getColMap()->getNodeNumElements();
897 Teuchos::Array<int> colflag(num_cols);
899 Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst();
905 Teuchos::ArrayView<const local_ordinal_type> UUI;
906 Teuchos::ArrayView<const scalar_type> UUV;
907 for (
size_t j = 0; j < num_cols; ++j) {
911 for (
size_t i = 0; i < L_->getNodeNumRows (); ++i) {
916 NumIn = MaxNumEntries;
917 L_->getLocalRowCopy (local_row, InI (), InV (), NumL);
920 InI[NumL] = local_row;
922 U_->getLocalRowCopy (local_row, InI (NumL+1, MaxNumEntries-NumL-1),
923 InV (NumL+1, MaxNumEntries-NumL-1), NumU);
927 for (
size_t j = 0; j < NumIn; ++j) {
933 for (
size_t jj = 0; jj < NumL; ++jj) {
939 U_->getLocalRowView(j, UUI, UUV);
942 if (RelaxValue_ == STM::zero ()) {
943 for (
size_t k = 0; k < NumUU; ++k) {
944 const int kk = colflag[UUI[k]];
949 InV[kk] -= multiplier * UUV[k];
954 for (
size_t k = 0; k < NumUU; ++k) {
958 const int kk = colflag[UUI[k]];
960 InV[kk] -= multiplier*UUV[k];
963 diagmod -= multiplier*UUV[k];
970 L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
975 if (RelaxValue_ != STM::zero ()) {
976 DV[i] += RelaxValue_*diagmod;
979 if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
980 if (STS::real (DV[i]) < STM::zero ()) {
981 DV[i] = -MinDiagonalValue;
984 DV[i] = MinDiagonalValue;
988 DV[i] = STS::one () / DV[i];
991 for (
size_t j = 0; j < NumU; ++j) {
992 InV[NumL+1+j] *= DV[i];
997 U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
1001 for (
size_t j = 0; j < NumIn; ++j) {
1002 colflag[InI[j]] = -1;
1013 L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
1014 U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
1019 TEUCHOS_TEST_FOR_EXCEPTION(
1020 0 < L_->getNodeNumRows() &&
1021 ! L_->isLowerTriangular (), std::runtime_error,
1022 "Ifpack2::RILUK::compute: L isn't lower triangular.");
1023 TEUCHOS_TEST_FOR_EXCEPTION(
1024 0 < U_->getNodeNumRows() &&
1025 ! U_->isUpperTriangular (), std::runtime_error,
1026 "Ifpack2::RILUK::compute: U isn't lower triangular.");
1029 L_solver_->initialize ();
1030 L_solver_->compute ();
1033 U_solver_->initialize ();
1034 U_solver_->compute ();
1037 #ifdef IFPACK2_ILUK_EXPERIMENTAL 1038 Teuchos::Time timerbasker (
"RILUK::basercompute");
1041 Teuchos::TimeMonitor timeMon (timer);
1043 if(basker_reuse ==
false)
1048 myBasker->Factor_Inc(0);
1049 basker_reuse =
true;
1055 myBasker->Options.same_pattern =
true;
1057 typedef typename crs_matrix_type::local_matrix_type kokkos_csr_matrix;
1058 kokkos_csr_matrix kcsr = A_local_crs_->getLocalMatrix();
1076 &(kcsr.graph.entries(0)),
1081 TEUCHOS_TEST_FOR_EXCEPTION(
1082 basker_error != 0, std::logic_error,
"Ifpack2::RILUK::initialize:" 1083 "Error in basker compute");
1088 # ifdef IFPACK2_HTS_EXPERIMENTAL 1091 Teuchos::Time basker_getL(
"basker_getL");
1092 Teuchos::Time hts_buildL (
"hts_buildL");
1093 Teuchos::Time basker_getU(
"basker_getU");
1094 Teuchos::Time hts_buildU (
"hts_buildU");
1096 hts_mgr_ = Teuchos::rcp(
new HTSManager());
1100 myBasker->GetPerm(&p, &q);
1103 myBasker->GetL(d.n, nnz, &d.ir, &d.jc, &d.v);
1105 typename HTST::CrsMatrix* T = HTST::make_CrsMatrix(d.n, d.ir, d.jc, d.v,
true);
1106 hts_mgr_->Limpl = HTST::preprocess(T, 1, hts_nthreads_,
true, p, 0);
1107 HTST::delete_CrsMatrix(T);
1112 myBasker->GetU(d.n, nnz, &d.ir, &d.jc, &d.v);
1114 typename HTST::CrsMatrix* T = HTST::make_CrsMatrix(d.n, d.ir, d.jc, d.v,
true);
1115 hts_mgr_->Uimpl = HTST::preprocess(T, 1, hts_nthreads_,
true, 0, q);
1116 HTST::delete_CrsMatrix(T);
1123 TEUCHOS_TEST_FOR_EXCEPTION(
1124 false, std::runtime_error,
1125 "Ifpack2::RILUK::compute: experimental not enabled");
1131 computeTime_ += timer.totalElapsedTime ();
1135 template<
class MatrixType>
1138 apply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1139 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1140 Teuchos::ETransp mode,
1145 using Teuchos::rcpFromRef;
1147 TEUCHOS_TEST_FOR_EXCEPTION(
1148 A_.is_null (), std::runtime_error,
"Ifpack2::RILUK::apply: The matrix is " 1149 "null. Please call setMatrix() with a nonnull input, then initialize() " 1150 "and compute(), before calling this method.");
1151 TEUCHOS_TEST_FOR_EXCEPTION(
1152 ! isComputed (), std::runtime_error,
1153 "Ifpack2::RILUK::apply: If you have not yet called compute(), " 1154 "you must call compute() before calling this method.");
1155 TEUCHOS_TEST_FOR_EXCEPTION(
1156 X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
1157 "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. " 1158 "X.getNumVectors() = " << X.getNumVectors ()
1159 <<
" != Y.getNumVectors() = " << Y.getNumVectors () <<
".");
1160 TEUCHOS_TEST_FOR_EXCEPTION(
1161 STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
1162 "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for " 1163 "complex Scalar type. Please talk to the Ifpack2 developers to get this " 1164 "fixed. There is a FIXME in this file about this very issue.");
1165 #ifdef HAVE_IFPACK2_DEBUG 1168 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.");
1169 Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
1172 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
1173 if (STM::isnaninf (norms[j])) {
1178 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
1180 #endif // HAVE_IFPACK2_DEBUG 1185 Teuchos::Time timer (
"RILUK::apply");
1187 if(!isExperimental_){
1188 Teuchos::TimeMonitor timeMon (timer);
1189 if (alpha == one && beta ==
zero) {
1190 if (mode == Teuchos::NO_TRANS) {
1197 MV C (D_->getMap (), X.getNumVectors ());
1198 L_solver_->apply (X, C, mode);
1203 C.elementWiseMultiply (one, *D_, C,
zero);
1205 U_solver_->apply (C, Y, mode);
1215 MV C (D_->getMap (), X.getNumVectors ());
1216 U_solver_->apply (X, C, mode);
1225 C.elementWiseMultiply (one, *D_, C,
zero);
1227 L_solver_->apply (C, Y, mode);
1231 if (alpha ==
zero) {
1238 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1239 apply (X, Y_tmp, mode);
1240 Y.update (alpha, Y_tmp, beta);
1246 #ifdef IFPACK2_ILUK_EXPERIMENTAL 1247 Teuchos::ArrayRCP<const scalar_type> XX;
1248 Teuchos::ArrayRCP<const scalar_type> YY;
1252 # ifdef IFPACK2_HTS_EXPERIMENTAL 1254 HTST::solve_omp(hts_mgr_->Limpl, const_cast<scalar_type*>(XX.getRawPtr()),
1255 X.getNumVectors(),
const_cast<scalar_type*
>(YY.getRawPtr()));
1256 HTST::solve_omp(hts_mgr_->Uimpl, const_cast<scalar_type*>(YY.getRawPtr()),
1262 (const_cast<scalar_type *>(XX.getRawPtr())),
1266 TEUCHOS_TEST_FOR_EXCEPTION(
1267 0==1, std::runtime_error,
1268 "Ifpack2::RILUK::apply: Experimental no enabled");
1273 #ifdef HAVE_IFPACK2_DEBUG 1275 Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
1278 for (
size_t j = 0; j < Y.getNumVectors (); ++j) {
1279 if (STM::isnaninf (norms[j])) {
1284 TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error,
"Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
1286 #endif // HAVE_IFPACK2_DEBUG 1289 applyTime_ = timer.totalElapsedTime ();
1293 template<
class MatrixType>
1295 multiply (
const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1296 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1297 const Teuchos::ETransp mode)
const 1299 const scalar_type
zero = STS::zero ();
1300 const scalar_type one = STS::one ();
1302 if (mode != Teuchos::NO_TRANS) {
1303 U_->apply (X, Y, mode);
1304 Y.update (one, X, one);
1309 Y.elementWiseMultiply (one, *D_, Y, zero);
1311 MV Y_tmp (Y, Teuchos::Copy);
1312 L_->apply (Y_tmp, Y, mode);
1313 Y.update (one, Y_tmp, one);
1316 L_->apply (X, Y, mode);
1317 Y.update (one, X, one);
1318 Y.elementWiseMultiply (one, *D_, Y, zero);
1319 MV Y_tmp (Y, Teuchos::Copy);
1320 U_->apply (Y_tmp, Y, mode);
1321 Y.update (one, Y_tmp, one);
1326 template<
class MatrixType>
1329 std::ostringstream os;
1334 os <<
"\"Ifpack2::RILUK\": {";
1335 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", " 1336 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
1338 os <<
"Level-of-fill: " << getLevelOfFill() <<
", ";
1340 if (A_.is_null ()) {
1341 os <<
"Matrix: null";
1344 os <<
"Global matrix dimensions: [" 1345 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]" 1353 #if defined IFPACK2_ILUK_EXPERIMENTAL && defined IFPACK2_HTS_EXPERIMENTAL 1354 template<
class MatrixType>
1356 : Limpl(0), Uimpl(0)
1359 template<
class MatrixType>
1360 RILUK<MatrixType>::HTSManager::~HTSManager ()
1362 HTST::delete_Impl(Limpl);
1363 HTST::delete_Impl(Uimpl);
1366 template<
class MatrixType>
1367 RILUK<MatrixType>::HTSData::HTSData ()
1368 : jc(0), ir(0), v(0)
1371 template<
class MatrixType>
1372 RILUK<MatrixType>::HTSData::~HTSData ()
1374 free_CrsMatrix_data();
1377 template<
class MatrixType>
1378 void RILUK<MatrixType>::HTSData::free_CrsMatrix_data ()
1380 if (jc)
delete[] jc;
1381 if (ir)
delete[] ir;
1387 template<
class MatrixType>
1388 void RILUK<MatrixType>::HTSData::sort ()
1390 if ( ! ir || ! jc)
return;
1391 std::vector<Entry> es;
1392 for (local_ordinal_type i = 0; i < n; ++i) {
1393 es.resize(ir[i+1] - ir[i]);
1394 const local_ordinal_type os = ir[i];
1395 for (local_ordinal_type j = 0; j < ir[i+1] - os; ++j) {
1399 std::sort(es.begin(), es.end());
1400 for (local_ordinal_type j = 0; j < ir[i+1] - os; ++j) {
1410 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \ 1411 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:272
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:284
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_RILUK_def.hpp:106
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:254
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:287
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:275
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:269
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:189
LinearOp zero(const VectorSpace &vs)
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:83
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:97
Tpetra::global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_RILUK_decl.hpp:524
Definition: Ifpack2_Container.hpp:761
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
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:208
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:266