43 #ifndef IFPACK2_ILUT_DEF_HPP
44 #define IFPACK2_ILUT_DEF_HPP
46 #include <type_traits>
47 #include "Kokkos_StaticCrsGraph.hpp"
49 #include "Teuchos_StandardParameterEntryValidators.hpp"
51 #include "Tpetra_CrsMatrix.hpp"
52 #include "KokkosSparse_par_ilut.hpp"
54 #include "Ifpack2_Heap.hpp"
55 #include "Ifpack2_LocalFilter.hpp"
56 #include "Ifpack2_LocalSparseTriangularSolver.hpp"
57 #include "Ifpack2_Parameters.hpp"
58 #include "Ifpack2_Details_getParamTryingTypes.hpp"
72 type_strs[0] =
"serial";
73 type_strs[1] =
"par_ilut";
75 type_enums[0] = Serial;
76 type_enums[1] = PAR_ILUT;
105 template<
class ScalarType>
107 ilutDefaultDropTolerance () {
109 typedef typename STS::magnitudeType magnitude_type;
113 const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
118 return std::min (static_cast<magnitude_type> (1000) * STS::magnitude (STS::eps ()), oneHalf);
125 ilutDefaultDropTolerance<double> () {
132 template <
class MatrixType>
135 Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
136 Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
137 RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
139 DropTolerance_ (ilutDefaultDropTolerance<
scalar_type> ()),
140 par_ilut_options_{1, 0., -1, -1, 0.75,
false},
141 InitializeTime_ (0.0),
147 IsInitialized_ (
false),
149 useKokkosKernelsParILUT_(
false)
155 template<
class MatrixType>
156 void ILUT<MatrixType>::allocateSolvers ()
158 L_solver_ =
Teuchos::rcp (
new LocalSparseTriangularSolver<row_matrix_type> ());
159 L_solver_->setObjectLabel(
"lower");
160 U_solver_ =
Teuchos::rcp (
new LocalSparseTriangularSolver<row_matrix_type> ());
161 U_solver_->setObjectLabel(
"upper");
164 template <
class MatrixType>
167 using Ifpack2::Details::getParamTryingTypes;
168 const char prefix[] =
"Ifpack2::ILUT: ";
175 IlutImplType::Enum ilutimplType = IlutImplType::Serial;
177 static const char typeName[] =
"fact: type";
179 if ( ! params.
isType<std::string>(typeName))
break;
184 IlutImplType::loadPLTypeOption (ilutimplTypeStrs, ilutimplTypeEnums);
186 s2i(ilutimplTypeStrs (), ilutimplTypeEnums (), typeName,
false);
191 if (ilutimplType == IlutImplType::PAR_ILUT) {
192 this->useKokkosKernelsParILUT_ =
true;
195 this->useKokkosKernelsParILUT_ =
false;
201 double fillLevel = LevelOfFill_;
203 const std::string paramName (
"fact: ilut level-of-fill");
205 (params.
isParameter(paramName) && this->useKokkosKernelsParILUT_), std::runtime_error,
206 "Ifpack2::ILUT: Parameter " << paramName <<
" is meaningless for algorithm par_ilut.");
207 getParamTryingTypes<double, double, float>
208 (fillLevel, params, paramName, prefix);
210 (fillLevel < 1.0, std::runtime_error,
211 "Ifpack2::ILUT: The \"" << paramName <<
"\" parameter must be >= "
212 "1.0, but you set it to " << fillLevel <<
". For ILUT, the fill level "
213 "means something different than it does for ILU(k). ILU(0) produces "
214 "factors with the same sparsity structure as the input matrix A. For "
215 "ILUT, level-of-fill = 1.0 will produce factors with nonzeros matching "
216 "the sparsity structure of A. level-of-fill > 1.0 allows for additional "
220 magnitude_type absThresh = Athresh_;
222 const std::string paramName (
"fact: absolute threshold");
223 getParamTryingTypes<magnitude_type, magnitude_type, double>
224 (absThresh, params, paramName, prefix);
227 magnitude_type relThresh = Rthresh_;
229 const std::string paramName (
"fact: relative threshold");
230 getParamTryingTypes<magnitude_type, magnitude_type, double>
231 (relThresh, params, paramName, prefix);
234 magnitude_type relaxValue = RelaxValue_;
236 const std::string paramName (
"fact: relax value");
237 getParamTryingTypes<magnitude_type, magnitude_type, double>
238 (relaxValue, params, paramName, prefix);
241 magnitude_type dropTol = DropTolerance_;
243 const std::string paramName (
"fact: drop tolerance");
244 getParamTryingTypes<magnitude_type, magnitude_type, double>
245 (dropTol, params, paramName, prefix);
248 int par_ilut_max_iter=20;
249 magnitude_type par_ilut_residual_norm_delta_stop=1e-2;
250 int par_ilut_team_size=0;
251 int par_ilut_vector_size=0;
252 float par_ilut_fill_in_limit=0.75;
253 bool par_ilut_verbose=
false;
254 if (this->useKokkosKernelsParILUT_) {
255 par_ilut_max_iter = par_ilut_options_.max_iter;
256 par_ilut_residual_norm_delta_stop = par_ilut_options_.residual_norm_delta_stop;
257 par_ilut_team_size = par_ilut_options_.team_size;
258 par_ilut_vector_size = par_ilut_options_.vector_size;
259 par_ilut_fill_in_limit = par_ilut_options_.fill_in_limit;
260 par_ilut_verbose = par_ilut_options_.verbose;
262 std::string par_ilut_plist_name(
"parallel ILUT options");
263 if (params.
isSublist(par_ilut_plist_name)) {
266 std::string paramName(
"maximum iterations");
267 getParamTryingTypes<int, int>(par_ilut_max_iter, par_ilut_plist, paramName, prefix);
269 paramName =
"residual norm delta stop";
270 getParamTryingTypes<magnitude_type, magnitude_type, double>(par_ilut_residual_norm_delta_stop, par_ilut_plist, paramName, prefix);
272 paramName =
"team size";
273 getParamTryingTypes<int, int>(par_ilut_team_size, par_ilut_plist, paramName, prefix);
275 paramName =
"vector size";
276 getParamTryingTypes<int, int>(par_ilut_vector_size, par_ilut_plist, paramName, prefix);
278 paramName =
"fill in limit";
279 getParamTryingTypes<float, float, double>(par_ilut_fill_in_limit, par_ilut_plist, paramName, prefix);
281 paramName =
"verbose";
282 getParamTryingTypes<bool, bool>(par_ilut_verbose, par_ilut_plist, paramName, prefix);
286 par_ilut_options_.max_iter = par_ilut_max_iter;
287 par_ilut_options_.residual_norm_delta_stop = par_ilut_residual_norm_delta_stop;
288 par_ilut_options_.team_size = par_ilut_team_size;
289 par_ilut_options_.vector_size = par_ilut_vector_size;
290 par_ilut_options_.fill_in_limit = par_ilut_fill_in_limit;
291 par_ilut_options_.verbose = par_ilut_verbose;
296 L_solver_->setParameters(params);
297 U_solver_->setParameters(params);
299 LevelOfFill_ = fillLevel;
300 Athresh_ = absThresh;
301 Rthresh_ = relThresh;
302 RelaxValue_ = relaxValue;
303 DropTolerance_ = dropTol;
307 template <
class MatrixType>
311 A_.
is_null (), std::runtime_error,
"Ifpack2::ILUT::getComm: "
312 "The matrix is null. Please call setMatrix() with a nonnull input "
313 "before calling this method.");
314 return A_->getComm ();
318 template <
class MatrixType>
325 template <
class MatrixType>
330 A_.
is_null (), std::runtime_error,
"Ifpack2::ILUT::getDomainMap: "
331 "The matrix is null. Please call setMatrix() with a nonnull input "
332 "before calling this method.");
333 return A_->getDomainMap ();
337 template <
class MatrixType>
342 A_.
is_null (), std::runtime_error,
"Ifpack2::ILUT::getRangeMap: "
343 "The matrix is null. Please call setMatrix() with a nonnull input "
344 "before calling this method.");
345 return A_->getRangeMap ();
349 template <
class MatrixType>
355 template <
class MatrixType>
357 return NumInitialize_;
361 template <
class MatrixType>
367 template <
class MatrixType>
373 template <
class MatrixType>
375 return InitializeTime_;
379 template<
class MatrixType>
385 template<
class MatrixType>
391 template<
class MatrixType>
394 A_.
is_null (), std::runtime_error,
"Ifpack2::ILUT::getNodeSmootherComplexity: "
395 "The input matrix A is null. Please call setMatrix() with a nonnull "
396 "input matrix, then call compute(), before calling this method.");
398 return A_->getLocalNumEntries() + getLocalNumEntries();
402 template<
class MatrixType>
404 return L_->getGlobalNumEntries () + U_->getGlobalNumEntries ();
408 template<
class MatrixType>
410 return L_->getLocalNumEntries () + U_->getLocalNumEntries ();
414 template<
class MatrixType>
420 ! A.
is_null () && A->getComm ()->getSize () == 1 &&
421 A->getLocalNumRows () != A->getLocalNumCols (),
422 std::runtime_error,
"Ifpack2::ILUT::setMatrix: If A's communicator only "
423 "contains one process, then A must be square. Instead, you provided a "
424 "matrix A with " << A->getLocalNumRows () <<
" rows and "
425 << A->getLocalNumCols () <<
" columns.");
431 IsInitialized_ =
false;
433 A_local_ = Teuchos::null;
440 if (! L_solver_.is_null ()) {
441 L_solver_->setMatrix (Teuchos::null);
443 if (! U_solver_.is_null ()) {
444 U_solver_->setMatrix (Teuchos::null);
453 template <
class MatrixType>
459 using Teuchos::rcp_dynamic_cast;
460 using Teuchos::rcp_implicit_cast;
465 if (A->getRowMap ()->getComm ()->getSize () == 1 ||
466 A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
473 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
474 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> > (A);
475 if (! A_lf_r.is_null ()) {
476 return rcp_implicit_cast<
const row_matrix_type> (A_lf_r);
482 return rcp (
new LocalFilter<row_matrix_type> (A));
487 template<
class MatrixType>
492 using Teuchos::rcp_const_cast;
494 double startTime = timer.
wallTime();
500 A_.
is_null (), std::runtime_error,
"Ifpack2::ILUT::initialize: "
501 "The matrix to precondition is null. Please call setMatrix() with a "
502 "nonnull input before calling this method.");
505 IsInitialized_ =
false;
507 A_local_ = Teuchos::null;
511 A_local_ = makeLocalFilter(A_);
513 A_local_.is_null(), std::logic_error,
"Ifpack2::RILUT::initialize: "
514 "makeLocalFilter returned null; it failed to compute A_local. "
515 "Please report this bug to the Ifpack2 developers.");
517 if (this->useKokkosKernelsParILUT_) {
518 this->KernelHandle_ =
Teuchos::rcp(
new kk_handle_type());
519 KernelHandle_->create_par_ilut_handle();
520 auto par_ilut_handle = KernelHandle_->get_par_ilut_handle();
521 par_ilut_handle->set_residual_norm_delta_stop(par_ilut_options_.residual_norm_delta_stop);
522 par_ilut_handle->set_team_size(par_ilut_options_.team_size);
523 par_ilut_handle->set_vector_size(par_ilut_options_.vector_size);
524 par_ilut_handle->set_fill_in_limit(par_ilut_options_.fill_in_limit);
525 par_ilut_handle->set_verbose(par_ilut_options_.verbose);
526 par_ilut_handle->set_async_update(
false);
528 RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_local_);
529 if (A_local_crs.is_null()) {
532 Array<size_t> entriesPerRow(numRows);
534 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
536 RCP<crs_matrix_type> A_local_crs_nc =
538 A_local_->getColMap (),
541 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getLocalMaxNumRowEntries());
542 nonconst_values_host_view_type values(
"values",A_local_->getLocalMaxNumRowEntries());
544 size_t numEntries = 0;
545 A_local_->getLocalRowCopy(i, indices, values, numEntries);
546 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
548 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
551 auto A_local_crs_device = A_local_crs->getLocalMatrixDevice();
554 typedef typename Kokkos::View<usize_type*, array_layout, device_type> ulno_row_view_t;
555 const int NumMyRows = A_local_crs->getRowMap()->getLocalNumElements();
556 L_rowmap_ = ulno_row_view_t(
"L_row_map", NumMyRows + 1);
557 U_rowmap_ = ulno_row_view_t(
"U_row_map", NumMyRows + 1);
559 KokkosSparse::Experimental::par_ilut_symbolic(KernelHandle_.getRawPtr(),
560 A_local_crs_device.graph.row_map, A_local_crs_device.graph.entries,
565 IsInitialized_ =
true;
568 InitializeTime_ += (timer.
wallTime() - startTime);
572 template<
typename ScalarType>
574 scalar_mag (
const ScalarType& s)
580 template<
class MatrixType>
588 using Teuchos::reduceAll;
590 using Teuchos::rcp_const_cast;
593 if (! isInitialized ()) {
598 double startTime = timer.
wallTime();
602 if (!this->useKokkosKernelsParILUT_)
635 #ifdef IFPACK2_WRITE_ILUT_FACTORS
636 std::ofstream ofsL(
"L.ifpack2_ilut.mtx", std::ios::out);
637 std::ofstream ofsU(
"U.ifpack2_ilut.mtx", std::ios::out);
642 double local_nnz =
static_cast<double> (A_local_->getLocalNumEntries ());
643 double fill = ((getLevelOfFill () - 1.0) * local_nnz) / (2 * myNumRows);
648 double fill_ceil=std::ceil(fill);
652 size_type fillL =
static_cast<size_type
>(fill_ceil);
653 size_type fillU =
static_cast<size_type
>(fill_ceil);
655 Array<scalar_type> InvDiagU (myNumRows, zero);
657 Array<Array<local_ordinal_type> > L_tmp_idx(myNumRows);
658 Array<Array<scalar_type> > L_tmpv(myNumRows);
659 Array<Array<local_ordinal_type> > U_tmp_idx(myNumRows);
660 Array<Array<scalar_type> > U_tmpv(myNumRows);
662 enum { UNUSED, ORIG, FILL };
665 Array<int> pattern(max_col, UNUSED);
666 Array<scalar_type> cur_row(max_col, zero);
667 Array<magnitude_type> unorm(max_col);
668 magnitude_type rownorm;
669 Array<local_ordinal_type> L_cols_heap;
670 Array<local_ordinal_type> U_cols;
671 Array<local_ordinal_type> L_vals_heap;
672 Array<local_ordinal_type> U_vals_heap;
677 greater_indirect<scalar_type,local_ordinal_type> vals_comp(cur_row);
682 nonconst_local_inds_host_view_type ColIndicesARCP;
683 nonconst_values_host_view_type ColValuesARCP;
684 if (! A_local_->supportsRowViews ()) {
685 const size_t maxnz = A_local_->getLocalMaxNumRowEntries ();
686 Kokkos::resize(ColIndicesARCP,maxnz);
687 Kokkos::resize(ColValuesARCP,maxnz);
691 local_inds_host_view_type ColIndicesA;
692 values_host_view_type ColValuesA;
695 if (A_local_->supportsRowViews ()) {
696 A_local_->getLocalRowView (row_i, ColIndicesA, ColValuesA);
697 RowNnz = ColIndicesA.size ();
700 A_local_->getLocalRowCopy (row_i, ColIndicesARCP, ColValuesARCP, RowNnz);
701 ColIndicesA = Kokkos::subview(ColIndicesARCP,std::make_pair((
size_t)0, RowNnz));
702 ColValuesA = Kokkos::subview(ColValuesARCP,std::make_pair((
size_t)0, RowNnz));
707 U_cols.push_back(row_i);
708 cur_row[row_i] = zero;
709 pattern[row_i] = ORIG;
711 size_type L_cols_heaplen = 0;
712 rownorm = STM::zero ();
713 for (
size_t i = 0; i < RowNnz; ++i) {
714 if (ColIndicesA[i] < myNumRows) {
715 if (ColIndicesA[i] < row_i) {
716 add_to_heap(ColIndicesA[i], L_cols_heap, L_cols_heaplen);
718 else if (ColIndicesA[i] > row_i) {
719 U_cols.push_back(ColIndicesA[i]);
722 cur_row[ColIndicesA[i]] = ColValuesA[i];
723 pattern[ColIndicesA[i]] = ORIG;
724 rownorm += scalar_mag(ColValuesA[i]);
731 const magnitude_type rthresh = getRelativeThreshold();
733 cur_row[row_i] = as<scalar_type> (getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh*v;
735 size_type orig_U_len = U_cols.size();
736 RowNnz = L_cols_heap.size() + orig_U_len;
737 rownorm = getDropTolerance() * rownorm/RowNnz;
740 size_type L_vals_heaplen = 0;
741 while (L_cols_heaplen > 0) {
744 scalar_type multiplier = cur_row[row_k] * InvDiagU[row_k];
745 cur_row[row_k] = multiplier;
746 magnitude_type mag_mult = scalar_mag(multiplier);
747 if (mag_mult*unorm[row_k] < rownorm) {
748 pattern[row_k] = UNUSED;
752 if (pattern[row_k] != ORIG) {
753 if (L_vals_heaplen < fillL) {
754 add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
756 else if (L_vals_heaplen==0 ||
757 mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
758 pattern[row_k] = UNUSED;
763 pattern[L_vals_heap.front()] = UNUSED;
765 add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
771 ArrayView<local_ordinal_type> ColIndicesU = U_tmp_idx[row_k]();
772 ArrayView<scalar_type> ColValuesU = U_tmpv[row_k]();
773 size_type ColNnzU = ColIndicesU.size();
775 for(size_type j=0; j<ColNnzU; ++j) {
776 if (ColIndicesU[j] > row_k) {
779 if (pattern[col_j] != UNUSED) {
780 cur_row[col_j] -= tmp;
782 else if (scalar_mag(tmp) > rownorm) {
783 cur_row[col_j] = -tmp;
784 pattern[col_j] = FILL;
786 U_cols.push_back(col_j);
802 for (size_type i = 0; i < (size_type)ColIndicesA.size (); ++i) {
803 if (ColIndicesA[i] < row_i) {
804 L_tmp_idx[row_i].push_back(ColIndicesA[i]);
805 L_tmpv[row_i].push_back(cur_row[ColIndicesA[i]]);
806 pattern[ColIndicesA[i]] = UNUSED;
811 for (size_type j = 0; j < L_vals_heaplen; ++j) {
812 L_tmp_idx[row_i].push_back(L_vals_heap[j]);
813 L_tmpv[row_i].push_back(cur_row[L_vals_heap[j]]);
814 pattern[L_vals_heap[j]] = UNUSED;
822 #ifdef IFPACK2_WRITE_ILUT_FACTORS
823 for (size_type ii = 0; ii < L_tmp_idx[row_i].size (); ++ii) {
824 ofsL << row_i <<
" " << L_tmp_idx[row_i][ii] <<
" "
825 << L_tmpv[row_i][ii] << std::endl;
831 if (cur_row[row_i] == zero) {
832 std::cerr <<
"Ifpack2::ILUT::Compute: zero pivot encountered! "
833 <<
"Replacing with rownorm and continuing..."
834 <<
"(You may need to set the parameter "
835 <<
"'fact: absolute threshold'.)" << std::endl;
836 cur_row[row_i] = rownorm;
838 InvDiagU[row_i] = one / cur_row[row_i];
841 U_tmp_idx[row_i].push_back(row_i);
842 U_tmpv[row_i].push_back(cur_row[row_i]);
843 unorm[row_i] = scalar_mag(cur_row[row_i]);
844 pattern[row_i] = UNUSED;
850 size_type U_vals_heaplen = 0;
851 for(size_type j=1; j<U_cols.size(); ++j) {
853 if (pattern[col] != ORIG) {
854 if (U_vals_heaplen < fillU) {
855 add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
857 else if (U_vals_heaplen!=0 && scalar_mag(cur_row[col]) >
858 scalar_mag(cur_row[U_vals_heap.front()])) {
860 add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
864 U_tmp_idx[row_i].push_back(col);
865 U_tmpv[row_i].push_back(cur_row[col]);
866 unorm[row_i] += scalar_mag(cur_row[col]);
868 pattern[col] = UNUSED;
871 for(size_type j=0; j<U_vals_heaplen; ++j) {
872 U_tmp_idx[row_i].push_back(U_vals_heap[j]);
873 U_tmpv[row_i].push_back(cur_row[U_vals_heap[j]]);
874 unorm[row_i] += scalar_mag(cur_row[U_vals_heap[j]]);
877 unorm[row_i] /= (orig_U_len + U_vals_heaplen);
879 #ifdef IFPACK2_WRITE_ILUT_FACTORS
880 for(
int ii=0; ii<U_tmp_idx[row_i].size(); ++ii) {
881 ofsU <<row_i<<
" " <<U_tmp_idx[row_i][ii]<<
" "
882 <<U_tmpv[row_i][ii]<< std::endl;
893 Array<size_t> nnzPerRow(myNumRows);
899 L_solver_->setMatrix(Teuchos::null);
900 U_solver_->setMatrix(Teuchos::null);
903 nnzPerRow[row_i] = L_tmp_idx[row_i].size();
910 L_->insertLocalValues (row_i, L_tmp_idx[row_i](), L_tmpv[row_i]());
916 nnzPerRow[row_i] = U_tmp_idx[row_i].size();
923 U_->insertLocalValues (row_i, U_tmp_idx[row_i](), U_tmpv[row_i]());
928 L_solver_->setMatrix(L_);
929 L_solver_->initialize ();
930 L_solver_->compute ();
932 U_solver_->setMatrix(U_);
933 U_solver_->initialize ();
934 U_solver_->compute ();
938 RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_local_);
940 if(A_local_crs.is_null()) {
942 Array<size_t> entriesPerRow(numRows);
944 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
946 RCP<crs_matrix_type> A_local_crs_nc =
948 A_local_->getColMap (),
951 nonconst_local_inds_host_view_type indices(
"indices",A_local_->getLocalMaxNumRowEntries());
952 nonconst_values_host_view_type values(
"values",A_local_->getLocalMaxNumRowEntries());
954 size_t numEntries = 0;
955 A_local_->getLocalRowCopy(i, indices, values, numEntries);
956 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()),indices.data());
958 A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
961 auto lclMtx = A_local_crs->getLocalMatrixDevice();
962 A_local_rowmap_ = lclMtx.graph.row_map;
963 A_local_entries_ = lclMtx.graph.entries;
964 A_local_values_ = lclMtx.values;
968 auto par_ilut_handle = KernelHandle_->get_par_ilut_handle();
969 auto nnzL = par_ilut_handle->get_nnzL();
970 static_graph_entries_t L_entries_ = static_graph_entries_t(
"L_entries", nnzL);
971 local_matrix_values_t L_values_ = local_matrix_values_t(
"L_values", nnzL);
973 auto nnzU = par_ilut_handle->get_nnzU();
974 static_graph_entries_t U_entries_ = static_graph_entries_t(
"U_entries", nnzU);
975 local_matrix_values_t U_values_ = local_matrix_values_t(
"U_values", nnzU);
977 KokkosSparse::Experimental::par_ilut_numeric(KernelHandle_.getRawPtr(),
978 A_local_rowmap_, A_local_entries_, A_local_values_,
979 L_rowmap_, L_entries_, L_values_, U_rowmap_, U_entries_, U_values_);
981 auto L_kokkosCrsGraph = local_graph_device_type(L_entries_, L_rowmap_);
982 auto U_kokkosCrsGraph = local_graph_device_type(U_entries_, U_rowmap_);
984 local_matrix_device_type L_localCrsMatrix_device;
985 L_localCrsMatrix_device = local_matrix_device_type(
"L_Factor_localmatrix",
986 A_local_->getLocalNumRows(),
991 A_local_crs->getRowMap(),
992 A_local_crs->getColMap(),
993 A_local_crs->getDomainMap(),
994 A_local_crs->getRangeMap(),
995 A_local_crs->getGraph()->getImporter(),
996 A_local_crs->getGraph()->getExporter()));
998 local_matrix_device_type U_localCrsMatrix_device;
999 U_localCrsMatrix_device = local_matrix_device_type(
"U_Factor_localmatrix",
1000 A_local_->getLocalNumRows(),
1005 A_local_crs->getRowMap(),
1006 A_local_crs->getColMap(),
1007 A_local_crs->getDomainMap(),
1008 A_local_crs->getRangeMap(),
1009 A_local_crs->getGraph()->getImporter(),
1010 A_local_crs->getGraph()->getExporter()));
1012 L_solver_->setMatrix (L_);
1013 L_solver_->compute ();
1014 U_solver_->setMatrix (U_);
1015 U_solver_->compute ();
1020 ComputeTime_ += (timer.
wallTime() - startTime);
1026 template <
class MatrixType>
1028 apply (
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
1029 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
1036 using Teuchos::rcpFromRef;
1039 ! isComputed (), std::runtime_error,
1040 "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
1041 "factorization, before calling apply().");
1044 X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
1045 "Ifpack2::ILUT::apply: X and Y must have the same number of columns. "
1046 "X has " << X.getNumVectors () <<
" columns, but Y has "
1047 << Y.getNumVectors () <<
" columns.");
1053 double startTime = timer.
wallTime();
1057 if (alpha == one && beta == zero) {
1060 L_solver_->apply (X, Y, mode);
1063 U_solver_->apply (Y, Y, mode);
1068 U_solver_->apply (X, Y, mode);
1071 L_solver_->apply (Y, Y, mode);
1075 if (alpha == zero) {
1085 MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1086 apply (X, Y_tmp, mode);
1087 Y.update (alpha, Y_tmp, beta);
1093 ApplyTime_ += (timer.
wallTime() - startTime);
1097 template <
class MatrixType>
1100 std::ostringstream os;
1105 os <<
"\"Ifpack2::ILUT\": {";
1106 os <<
"Initialized: " << (isInitialized () ?
"true" :
"false") <<
", "
1107 <<
"Computed: " << (isComputed () ?
"true" :
"false") <<
", ";
1109 os <<
"Level-of-fill: " << getLevelOfFill() <<
", "
1110 <<
"absolute threshold: " << getAbsoluteThreshold() <<
", "
1111 <<
"relative threshold: " << getRelativeThreshold() <<
", "
1112 <<
"relaxation value: " << getRelaxValue() <<
", ";
1115 os <<
"Matrix: null";
1118 os <<
"Global matrix dimensions: ["
1119 << A_->getGlobalNumRows () <<
", " << A_->getGlobalNumCols () <<
"]"
1120 <<
", Global nnz: " << A_->getGlobalNumEntries();
1128 template <
class MatrixType>
1151 out <<
"\"Ifpack2::ILUT\":" << endl;
1153 out <<
"MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
1154 if (this->getObjectLabel () !=
"") {
1155 out <<
"Label: \"" << this->getObjectLabel () <<
"\"" << endl;
1157 out <<
"Initialized: " << (isInitialized () ?
"true" :
"false")
1159 <<
"Computed: " << (isComputed () ?
"true" :
"false")
1161 <<
"Level of fill: " << getLevelOfFill () << endl
1162 <<
"Absolute threshold: " << getAbsoluteThreshold () << endl
1163 <<
"Relative threshold: " << getRelativeThreshold () << endl
1164 <<
"Relax value: " << getRelaxValue () << endl;
1167 const double fillFraction =
1168 (double) getGlobalNumEntries () / (double) A_->getGlobalNumEntries ();
1169 const double nnzToRows =
1170 (double) getGlobalNumEntries () / (double) U_->getGlobalNumRows ();
1172 out <<
"Dimensions of L: [" << L_->getGlobalNumRows () <<
", "
1173 << L_->getGlobalNumRows () <<
"]" << endl
1174 <<
"Dimensions of U: [" << U_->getGlobalNumRows () <<
", "
1175 << U_->getGlobalNumRows () <<
"]" << endl
1176 <<
"Number of nonzeros in factors: " << getGlobalNumEntries () << endl
1177 <<
"Fill fraction of factors over A: " << fillFraction << endl
1178 <<
"Ratio of nonzeros to rows: " << nnzToRows << endl;
1181 out <<
"Number of initialize calls: " << getNumInitialize () << endl
1182 <<
"Number of compute calls: " << getNumCompute () << endl
1183 <<
"Number of apply calls: " << getNumApply () << endl
1184 <<
"Total time in seconds for initialize: " << getInitializeTime () << endl
1185 <<
"Total time in seconds for compute: " << getComputeTime () << endl
1186 <<
"Total time in seconds for apply: " << getApplyTime () << endl;
1188 out <<
"Local matrix:" << endl;
1189 A_local_->describe (out, vl);
1201 #define IFPACK2_ILUT_INSTANT(S,LO,GO,N) \
1202 template class Ifpack2::ILUT< Tpetra::RowMatrix<S, LO, GO, N> >;
ILUT(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_ILUT_def.hpp:133
bool hasTransposeApply() const
Whether this object's apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_ILUT_def.hpp:350
basic_OSTab< char > OSTab
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_ILUT_def.hpp:403
void initialize()
Clear any previously computed factors, and potentially compute sparsity patterns of factors...
Definition: Ifpack2_ILUT_def.hpp:488
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_ILUT_def.hpp:1131
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_ILUT_def.hpp:1098
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:93
void rm_heap_root(Teuchos::Array< Ordinal > &heap, SizeType &heap_len)
Definition: Ifpack2_Heap.hpp:92
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isSublist(const std::string &name) const
Teuchos::RCP< const map_type > getRangeMap() const
Tpetra::Map representing the range of this operator.
Definition: Ifpack2_ILUT_def.hpp:339
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_ILUT_decl.hpp:108
void resize(size_type new_size, const value_type &x=value_type())
void compute()
Compute factors L and U using the specified diagonal perturbation thresholds and relaxation parameter...
Definition: Ifpack2_ILUT_def.hpp:581
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:415
IntegralType getIntegralValue(const std::string &str, const std::string ¶mName="", const std::string &sublistName="") const
static magnitudeType magnitude(T a)
int getNumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack2_ILUT_def.hpp:356
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_ILUT_def.hpp:386
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the input matrix's communicator.
Definition: Ifpack2_ILUT_def.hpp:309
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Type of the Tpetra::CrsMatrix specialization that this class uses for the L and U factors...
Definition: Ifpack2_ILUT_decl.hpp:136
TypeTo as(const TypeFrom &t)
Teuchos::RCP< const map_type > getDomainMap() const
Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_ILUT_def.hpp:327
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
void add_to_heap(const Ordinal &idx, Teuchos::Array< Ordinal > &heap, SizeType &heap_len)
Definition: Ifpack2_Heap.hpp:70
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_ILUT_def.hpp:392
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_ILUT_def.hpp:374
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_ILUT_def.hpp:368
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 ILUT preconditioner to X, resulting in Y.
Definition: Ifpack2_ILUT_def.hpp:1028
size_t getLocalNumEntries() const
Returns the number of nonzero entries in the local graph.
Definition: Ifpack2_ILUT_def.hpp:409
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_ILUT_def.hpp:380
void setParameters(const Teuchos::ParameterList ¶ms)
Set preconditioner parameters.
Definition: Ifpack2_ILUT_def.hpp:165
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_ILUT_def.hpp:362
std::string typeName(const T &t)
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_ILUT_decl.hpp:111
Teuchos::RCP< const row_matrix_type > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:320