10 #ifndef IFPACK2_ILUT_DEF_HPP
11 #define IFPACK2_ILUT_DEF_HPP
13 #include <type_traits>
15 #include "Teuchos_StandardParameterEntryValidators.hpp"
17 #include "Tpetra_CrsMatrix.hpp"
18 #include "KokkosSparse_par_ilut.hpp"
20 #include "Ifpack2_Heap.hpp"
21 #include "Ifpack2_LocalFilter.hpp"
22 #include "Ifpack2_LocalSparseTriangularSolver.hpp"
23 #include "Ifpack2_Parameters.hpp"
24 #include "Ifpack2_Details_getParamTryingTypes.hpp"
38 type_strs[0] =
"serial";
39 type_strs[1] =
"par_ilut";
41 type_enums[0] = Serial;
42 type_enums[1] = PAR_ILUT;
70 template <
class ScalarType>
72 ilutDefaultDropTolerance() {
74 typedef typename STS::magnitudeType magnitude_type;
78 const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
83 return std::min(static_cast<magnitude_type>(1000) * STS::magnitude(STS::eps()), oneHalf);
90 ilutDefaultDropTolerance<double>() {
96 template <
class MatrixType>
99 , Athresh_(Teuchos::ScalarTraits<magnitude_type>::zero())
100 , Rthresh_(Teuchos::ScalarTraits<magnitude_type>::one())
101 , RelaxValue_(Teuchos::ScalarTraits<magnitude_type>::zero())
103 , DropTolerance_(ilutDefaultDropTolerance<
scalar_type>())
104 , par_ilut_options_{1, 0., -1, -1, 0.75,
false}
105 , InitializeTime_(0.0)
111 , IsInitialized_(
false)
113 , useKokkosKernelsParILUT_(
false)
119 template <
class MatrixType>
120 void ILUT<MatrixType>::allocateSolvers() {
121 L_solver_ =
Teuchos::rcp(
new LocalSparseTriangularSolver<row_matrix_type>());
122 L_solver_->setObjectLabel(
"lower");
123 U_solver_ =
Teuchos::rcp(
new LocalSparseTriangularSolver<row_matrix_type>());
124 U_solver_->setObjectLabel(
"upper");
127 template <
class MatrixType>
129 using Ifpack2::Details::getParamTryingTypes;
130 const char prefix[] =
"Ifpack2::ILUT: ";
137 IlutImplType::Enum ilutimplType = IlutImplType::Serial;
139 static const char typeName[] =
"fact: type";
141 if (!params.
isType<std::string>(typeName))
break;
146 IlutImplType::loadPLTypeOption(ilutimplTypeStrs, ilutimplTypeEnums);
148 s2i(ilutimplTypeStrs(), ilutimplTypeEnums(), typeName,
false);
153 if (ilutimplType == IlutImplType::PAR_ILUT) {
154 this->useKokkosKernelsParILUT_ =
true;
156 this->useKokkosKernelsParILUT_ =
false;
162 double fillLevel = LevelOfFill_;
164 const std::string paramName(
"fact: ilut level-of-fill");
166 (params.
isParameter(paramName) && this->useKokkosKernelsParILUT_), std::runtime_error,
167 "Ifpack2::ILUT: Parameter " << paramName <<
" is meaningless for algorithm par_ilut.");
168 getParamTryingTypes<double, double, float>(fillLevel, params, paramName, prefix);
170 "Ifpack2::ILUT: The \"" << paramName <<
"\" parameter must be >= "
171 "1.0, but you set it to "
172 << fillLevel <<
". For ILUT, the fill level "
173 "means something different than it does for ILU(k). ILU(0) produces "
174 "factors with the same sparsity structure as the input matrix A. For "
175 "ILUT, level-of-fill = 1.0 will produce factors with nonzeros matching "
176 "the sparsity structure of A. level-of-fill > 1.0 allows for additional "
180 magnitude_type absThresh = Athresh_;
182 const std::string paramName(
"fact: absolute threshold");
183 getParamTryingTypes<magnitude_type, magnitude_type, double>(absThresh, params, paramName, prefix);
186 magnitude_type relThresh = Rthresh_;
188 const std::string paramName(
"fact: relative threshold");
189 getParamTryingTypes<magnitude_type, magnitude_type, double>(relThresh, params, paramName, prefix);
192 magnitude_type relaxValue = RelaxValue_;
194 const std::string paramName(
"fact: relax value");
195 getParamTryingTypes<magnitude_type, magnitude_type, double>(relaxValue, params, paramName, prefix);
198 magnitude_type dropTol = DropTolerance_;
200 const std::string paramName(
"fact: drop tolerance");
201 getParamTryingTypes<magnitude_type, magnitude_type, double>(dropTol, params, paramName, prefix);
204 int par_ilut_max_iter = 20;
205 magnitude_type par_ilut_residual_norm_delta_stop = 1e-2;
206 int par_ilut_team_size = 0;
207 int par_ilut_vector_size = 0;
208 float par_ilut_fill_in_limit = 0.75;
209 bool par_ilut_verbose =
false;
210 if (this->useKokkosKernelsParILUT_) {
211 par_ilut_max_iter = par_ilut_options_.max_iter;
212 par_ilut_residual_norm_delta_stop = par_ilut_options_.residual_norm_delta_stop;
213 par_ilut_team_size = par_ilut_options_.team_size;
214 par_ilut_vector_size = par_ilut_options_.vector_size;
215 par_ilut_fill_in_limit = par_ilut_options_.fill_in_limit;
216 par_ilut_verbose = par_ilut_options_.verbose;
218 std::string par_ilut_plist_name(
"parallel ILUT options");
219 if (params.
isSublist(par_ilut_plist_name)) {
222 std::string paramName(
"maximum iterations");
223 getParamTryingTypes<int, int>(par_ilut_max_iter, par_ilut_plist, paramName, prefix);
225 paramName =
"residual norm delta stop";
226 getParamTryingTypes<magnitude_type, magnitude_type, double>(par_ilut_residual_norm_delta_stop, par_ilut_plist, paramName, prefix);
228 paramName =
"team size";
229 getParamTryingTypes<int, int>(par_ilut_team_size, par_ilut_plist, paramName, prefix);
231 paramName =
"vector size";
232 getParamTryingTypes<int, int>(par_ilut_vector_size, par_ilut_plist, paramName, prefix);
234 paramName =
"fill in limit";
235 getParamTryingTypes<float, float, double>(par_ilut_fill_in_limit, par_ilut_plist, paramName, prefix);
237 paramName =
"verbose";
238 getParamTryingTypes<bool, bool>(par_ilut_verbose, par_ilut_plist, paramName, prefix);
242 par_ilut_options_.max_iter = par_ilut_max_iter;
243 par_ilut_options_.residual_norm_delta_stop = par_ilut_residual_norm_delta_stop;
244 par_ilut_options_.team_size = par_ilut_team_size;
245 par_ilut_options_.vector_size = par_ilut_vector_size;
246 par_ilut_options_.fill_in_limit = par_ilut_fill_in_limit;
247 par_ilut_options_.verbose = par_ilut_verbose;
252 L_solver_->setParameters(params);
253 U_solver_->setParameters(params);
255 LevelOfFill_ = fillLevel;
256 Athresh_ = absThresh;
257 Rthresh_ = relThresh;
258 RelaxValue_ = relaxValue;
259 DropTolerance_ = dropTol;
262 template <
class MatrixType>
266 A_.
is_null(), std::runtime_error,
267 "Ifpack2::ILUT::getComm: "
268 "The matrix is null. Please call setMatrix() with a nonnull input "
269 "before calling this method.");
270 return A_->getComm();
273 template <
class MatrixType>
279 template <
class MatrixType>
283 A_.
is_null(), std::runtime_error,
284 "Ifpack2::ILUT::getDomainMap: "
285 "The matrix is null. Please call setMatrix() with a nonnull input "
286 "before calling this method.");
287 return A_->getDomainMap();
290 template <
class MatrixType>
294 A_.
is_null(), std::runtime_error,
295 "Ifpack2::ILUT::getRangeMap: "
296 "The matrix is null. Please call setMatrix() with a nonnull input "
297 "before calling this method.");
298 return A_->getRangeMap();
301 template <
class MatrixType>
306 template <
class MatrixType>
308 return NumInitialize_;
311 template <
class MatrixType>
316 template <
class MatrixType>
321 template <
class MatrixType>
323 return InitializeTime_;
326 template <
class MatrixType>
331 template <
class MatrixType>
336 template <
class MatrixType>
339 A_.
is_null(), std::runtime_error,
340 "Ifpack2::ILUT::getNodeSmootherComplexity: "
341 "The input matrix A is null. Please call setMatrix() with a nonnull "
342 "input matrix, then call compute(), before calling this method.");
344 return A_->getLocalNumEntries() + getLocalNumEntries();
347 template <
class MatrixType>
349 return L_->getGlobalNumEntries() + U_->getGlobalNumEntries();
352 template <
class MatrixType>
354 return L_->getLocalNumEntries() + U_->getLocalNumEntries();
357 template <
class MatrixType>
362 !A.
is_null() && A->getComm()->getSize() == 1 &&
363 A->getLocalNumRows() != A->getLocalNumCols(),
365 "Ifpack2::ILUT::setMatrix: If A's communicator only "
366 "contains one process, then A must be square. Instead, you provided a "
368 << A->getLocalNumRows() <<
" rows and "
369 << A->getLocalNumCols() <<
" columns.");
375 IsInitialized_ =
false;
377 A_local_ = Teuchos::null;
384 if (!L_solver_.is_null()) {
385 L_solver_->setMatrix(Teuchos::null);
387 if (!U_solver_.is_null()) {
388 U_solver_->setMatrix(Teuchos::null);
397 template <
class MatrixType>
402 using Teuchos::rcp_dynamic_cast;
403 using Teuchos::rcp_implicit_cast;
408 if (A->getRowMap()->getComm()->getSize() == 1 ||
409 A->getRowMap()->isSameAs(*(A->getColMap()))) {
416 RCP<const LocalFilter<row_matrix_type> > A_lf_r =
417 rcp_dynamic_cast<
const LocalFilter<row_matrix_type> >(A);
418 if (!A_lf_r.is_null()) {
419 return rcp_implicit_cast<
const row_matrix_type>(A_lf_r);
424 return rcp(
new LocalFilter<row_matrix_type>(A));
428 template <
class MatrixType>
432 using Teuchos::rcp_const_cast;
434 double startTime = timer.
wallTime();
440 A_.
is_null(), std::runtime_error,
441 "Ifpack2::ILUT::initialize: "
442 "The matrix to precondition is null. Please call setMatrix() with a "
443 "nonnull input before calling this method.");
446 IsInitialized_ =
false;
448 A_local_ = Teuchos::null;
452 A_local_ = makeLocalFilter(A_);
454 A_local_.is_null(), std::logic_error,
455 "Ifpack2::RILUT::initialize: "
456 "makeLocalFilter returned null; it failed to compute A_local. "
457 "Please report this bug to the Ifpack2 developers.");
459 if (this->useKokkosKernelsParILUT_) {
460 this->KernelHandle_ =
Teuchos::rcp(
new kk_handle_type());
461 KernelHandle_->create_par_ilut_handle();
462 auto par_ilut_handle = KernelHandle_->get_par_ilut_handle();
463 par_ilut_handle->set_residual_norm_delta_stop(par_ilut_options_.residual_norm_delta_stop);
464 par_ilut_handle->set_team_size(par_ilut_options_.team_size);
465 par_ilut_handle->set_vector_size(par_ilut_options_.vector_size);
466 par_ilut_handle->set_fill_in_limit(par_ilut_options_.fill_in_limit);
467 par_ilut_handle->set_verbose(par_ilut_options_.verbose);
468 par_ilut_handle->set_async_update(
false);
470 RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_local_);
471 if (A_local_crs.is_null()) {
474 Array<size_t> entriesPerRow(numRows);
476 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
478 RCP<crs_matrix_type> A_local_crs_nc =
480 A_local_->getColMap(),
483 nonconst_local_inds_host_view_type indices(
"indices", A_local_->getLocalMaxNumRowEntries());
484 nonconst_values_host_view_type values(
"values", A_local_->getLocalMaxNumRowEntries());
486 size_t numEntries = 0;
487 A_local_->getLocalRowCopy(i, indices, values, numEntries);
488 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
490 A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
493 auto A_local_crs_device = A_local_crs->getLocalMatrixDevice();
496 typedef typename Kokkos::View<usize_type*, array_layout, device_type> ulno_row_view_t;
497 const int NumMyRows = A_local_crs->getRowMap()->getLocalNumElements();
498 L_rowmap_ = ulno_row_view_t(
"L_row_map", NumMyRows + 1);
499 U_rowmap_ = ulno_row_view_t(
"U_row_map", NumMyRows + 1);
500 L_rowmap_orig_ = ulno_row_view_t(
"L_row_map_orig", NumMyRows + 1);
501 U_rowmap_orig_ = ulno_row_view_t(
"U_row_map_orig", NumMyRows + 1);
503 KokkosSparse::Experimental::par_ilut_symbolic(KernelHandle_.getRawPtr(),
504 A_local_crs_device.graph.row_map, A_local_crs_device.graph.entries,
508 Kokkos::deep_copy(L_rowmap_orig_, L_rowmap_);
509 Kokkos::deep_copy(U_rowmap_orig_, U_rowmap_);
512 IsInitialized_ =
true;
515 InitializeTime_ += (timer.
wallTime() - startTime);
518 template <
typename ScalarType>
520 scalar_mag(
const ScalarType& s) {
524 template <
class MatrixType>
532 using Teuchos::rcp_const_cast;
533 using Teuchos::reduceAll;
536 if (!isInitialized()) {
541 double startTime = timer.
wallTime();
545 if (!this->useKokkosKernelsParILUT_) {
577 #ifdef IFPACK2_WRITE_ILUT_FACTORS
578 std::ofstream ofsL(
"L.ifpack2_ilut.mtx", std::ios::out);
579 std::ofstream ofsU(
"U.ifpack2_ilut.mtx", std::ios::out);
584 double local_nnz =
static_cast<double>(A_local_->getLocalNumEntries());
585 double fill = ((getLevelOfFill() - 1.0) * local_nnz) / (2 * myNumRows);
590 double fill_ceil = std::ceil(fill);
594 size_type fillL =
static_cast<size_type
>(fill_ceil);
595 size_type fillU =
static_cast<size_type
>(fill_ceil);
597 Array<scalar_type> InvDiagU(myNumRows, zero);
599 Array<Array<local_ordinal_type> > L_tmp_idx(myNumRows);
600 Array<Array<scalar_type> > L_tmpv(myNumRows);
601 Array<Array<local_ordinal_type> > U_tmp_idx(myNumRows);
602 Array<Array<scalar_type> > U_tmpv(myNumRows);
609 Array<int> pattern(max_col, UNUSED);
610 Array<scalar_type> cur_row(max_col, zero);
611 Array<magnitude_type> unorm(max_col);
612 magnitude_type rownorm;
613 Array<local_ordinal_type> L_cols_heap;
614 Array<local_ordinal_type> U_cols;
615 Array<local_ordinal_type> L_vals_heap;
616 Array<local_ordinal_type> U_vals_heap;
621 greater_indirect<scalar_type, local_ordinal_type> vals_comp(cur_row);
626 nonconst_local_inds_host_view_type ColIndicesARCP;
627 nonconst_values_host_view_type ColValuesARCP;
628 if (!A_local_->supportsRowViews()) {
629 const size_t maxnz = A_local_->getLocalMaxNumRowEntries();
630 Kokkos::resize(ColIndicesARCP, maxnz);
631 Kokkos::resize(ColValuesARCP, maxnz);
635 local_inds_host_view_type ColIndicesA;
636 values_host_view_type ColValuesA;
639 if (A_local_->supportsRowViews()) {
640 A_local_->getLocalRowView(row_i, ColIndicesA, ColValuesA);
641 RowNnz = ColIndicesA.size();
643 A_local_->getLocalRowCopy(row_i, ColIndicesARCP, ColValuesARCP, RowNnz);
644 ColIndicesA = Kokkos::subview(ColIndicesARCP, std::make_pair((
size_t)0, RowNnz));
645 ColValuesA = Kokkos::subview(ColValuesARCP, std::make_pair((
size_t)0, RowNnz));
650 U_cols.push_back(row_i);
651 cur_row[row_i] = zero;
652 pattern[row_i] = ORIG;
654 size_type L_cols_heaplen = 0;
655 rownorm = STM::zero();
656 for (
size_t i = 0; i < RowNnz; ++i) {
657 if (ColIndicesA[i] < myNumRows) {
658 if (ColIndicesA[i] < row_i) {
659 add_to_heap(ColIndicesA[i], L_cols_heap, L_cols_heaplen);
660 }
else if (ColIndicesA[i] > row_i) {
661 U_cols.push_back(ColIndicesA[i]);
664 cur_row[ColIndicesA[i]] = ColValuesA[i];
665 pattern[ColIndicesA[i]] = ORIG;
666 rownorm += scalar_mag(ColValuesA[i]);
673 const magnitude_type rthresh = getRelativeThreshold();
675 cur_row[row_i] = as<scalar_type>(getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh * v;
677 size_type orig_U_len = U_cols.size();
678 RowNnz = L_cols_heap.size() + orig_U_len;
679 rownorm = getDropTolerance() * rownorm / RowNnz;
682 size_type L_vals_heaplen = 0;
683 while (L_cols_heaplen > 0) {
686 scalar_type multiplier = cur_row[row_k] * InvDiagU[row_k];
687 cur_row[row_k] = multiplier;
688 magnitude_type mag_mult = scalar_mag(multiplier);
689 if (mag_mult * unorm[row_k] < rownorm) {
690 pattern[row_k] = UNUSED;
694 if (pattern[row_k] != ORIG) {
695 if (L_vals_heaplen < fillL) {
696 add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
697 }
else if (L_vals_heaplen == 0 ||
698 mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
699 pattern[row_k] = UNUSED;
703 pattern[L_vals_heap.front()] = UNUSED;
705 add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
711 ArrayView<local_ordinal_type> ColIndicesU = U_tmp_idx[row_k]();
712 ArrayView<scalar_type> ColValuesU = U_tmpv[row_k]();
713 size_type ColNnzU = ColIndicesU.size();
715 for (size_type j = 0; j < ColNnzU; ++j) {
716 if (ColIndicesU[j] > row_k) {
719 if (pattern[col_j] != UNUSED) {
720 cur_row[col_j] -= tmp;
721 }
else if (scalar_mag(tmp) > rownorm) {
722 cur_row[col_j] = -tmp;
723 pattern[col_j] = FILL;
725 U_cols.push_back(col_j);
739 for (size_type i = 0; i < (size_type)ColIndicesA.size(); ++i) {
740 if (ColIndicesA[i] < row_i) {
741 L_tmp_idx[row_i].push_back(ColIndicesA[i]);
742 L_tmpv[row_i].push_back(cur_row[ColIndicesA[i]]);
743 pattern[ColIndicesA[i]] = UNUSED;
748 for (size_type j = 0; j < L_vals_heaplen; ++j) {
749 L_tmp_idx[row_i].push_back(L_vals_heap[j]);
750 L_tmpv[row_i].push_back(cur_row[L_vals_heap[j]]);
751 pattern[L_vals_heap[j]] = UNUSED;
759 #ifdef IFPACK2_WRITE_ILUT_FACTORS
760 for (size_type ii = 0; ii < L_tmp_idx[row_i].size(); ++ii) {
761 ofsL << row_i <<
" " << L_tmp_idx[row_i][ii] <<
" "
762 << L_tmpv[row_i][ii] << std::endl;
767 if (cur_row[row_i] == zero) {
768 std::cerr <<
"Ifpack2::ILUT::Compute: zero pivot encountered! "
769 <<
"Replacing with rownorm and continuing..."
770 <<
"(You may need to set the parameter "
771 <<
"'fact: absolute threshold'.)" << std::endl;
772 cur_row[row_i] = rownorm;
774 InvDiagU[row_i] = one / cur_row[row_i];
777 U_tmp_idx[row_i].push_back(row_i);
778 U_tmpv[row_i].push_back(cur_row[row_i]);
779 unorm[row_i] = scalar_mag(cur_row[row_i]);
780 pattern[row_i] = UNUSED;
786 size_type U_vals_heaplen = 0;
787 for (size_type j = 1; j < U_cols.size(); ++j) {
789 if (pattern[col] != ORIG) {
790 if (U_vals_heaplen < fillU) {
791 add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
792 }
else if (U_vals_heaplen != 0 && scalar_mag(cur_row[col]) >
793 scalar_mag(cur_row[U_vals_heap.front()])) {
795 add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
798 U_tmp_idx[row_i].push_back(col);
799 U_tmpv[row_i].push_back(cur_row[col]);
800 unorm[row_i] += scalar_mag(cur_row[col]);
802 pattern[col] = UNUSED;
805 for (size_type j = 0; j < U_vals_heaplen; ++j) {
806 U_tmp_idx[row_i].push_back(U_vals_heap[j]);
807 U_tmpv[row_i].push_back(cur_row[U_vals_heap[j]]);
808 unorm[row_i] += scalar_mag(cur_row[U_vals_heap[j]]);
811 unorm[row_i] /= (orig_U_len + U_vals_heaplen);
813 #ifdef IFPACK2_WRITE_ILUT_FACTORS
814 for (
int ii = 0; ii < U_tmp_idx[row_i].size(); ++ii) {
815 ofsU << row_i <<
" " << U_tmp_idx[row_i][ii] <<
" "
816 << U_tmpv[row_i][ii] << std::endl;
827 Array<size_t> nnzPerRow(myNumRows);
833 L_solver_->setMatrix(Teuchos::null);
834 U_solver_->setMatrix(Teuchos::null);
837 nnzPerRow[row_i] = L_tmp_idx[row_i].size();
844 L_->insertLocalValues(row_i, L_tmp_idx[row_i](), L_tmpv[row_i]());
850 nnzPerRow[row_i] = U_tmp_idx[row_i].size();
857 U_->insertLocalValues(row_i, U_tmp_idx[row_i](), U_tmpv[row_i]());
862 L_solver_->setMatrix(L_);
863 L_solver_->initialize();
864 L_solver_->compute();
866 U_solver_->setMatrix(U_);
867 U_solver_->initialize();
868 U_solver_->compute();
874 if (this->isComputed()) {
875 Kokkos::resize(L_rowmap_, L_rowmap_orig_.size());
876 Kokkos::resize(U_rowmap_, U_rowmap_orig_.size());
877 Kokkos::deep_copy(L_rowmap_, L_rowmap_orig_);
878 Kokkos::deep_copy(U_rowmap_, U_rowmap_orig_);
881 RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_local_);
883 if (A_local_crs.is_null()) {
885 Array<size_t> entriesPerRow(numRows);
887 entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
889 RCP<crs_matrix_type> A_local_crs_nc =
891 A_local_->getColMap(),
894 nonconst_local_inds_host_view_type indices(
"indices", A_local_->getLocalMaxNumRowEntries());
895 nonconst_values_host_view_type values(
"values", A_local_->getLocalMaxNumRowEntries());
897 size_t numEntries = 0;
898 A_local_->getLocalRowCopy(i, indices, values, numEntries);
899 A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
901 A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
904 auto lclMtx = A_local_crs->getLocalMatrixDevice();
905 A_local_rowmap_ = lclMtx.graph.row_map;
906 A_local_entries_ = lclMtx.graph.entries;
907 A_local_values_ = lclMtx.values;
911 auto par_ilut_handle = KernelHandle_->get_par_ilut_handle();
912 auto nnzL = par_ilut_handle->get_nnzL();
913 static_graph_entries_t L_entries_ = static_graph_entries_t(
"L_entries", nnzL);
914 local_matrix_values_t L_values_ = local_matrix_values_t(
"L_values", nnzL);
916 auto nnzU = par_ilut_handle->get_nnzU();
917 static_graph_entries_t U_entries_ = static_graph_entries_t(
"U_entries", nnzU);
918 local_matrix_values_t U_values_ = local_matrix_values_t(
"U_values", nnzU);
920 KokkosSparse::Experimental::par_ilut_numeric(KernelHandle_.getRawPtr(),
921 A_local_rowmap_, A_local_entries_, A_local_values_,
922 L_rowmap_, L_entries_, L_values_, U_rowmap_, U_entries_, U_values_);
924 auto L_kokkosCrsGraph = local_graph_device_type(L_entries_, L_rowmap_);
925 auto U_kokkosCrsGraph = local_graph_device_type(U_entries_, U_rowmap_);
927 local_matrix_device_type L_localCrsMatrix_device;
928 L_localCrsMatrix_device = local_matrix_device_type(
"L_Factor_localmatrix",
929 A_local_->getLocalNumRows(),
934 A_local_crs->getRowMap(),
935 A_local_crs->getColMap(),
936 A_local_crs->getDomainMap(),
937 A_local_crs->getRangeMap(),
938 A_local_crs->getGraph()->getImporter(),
939 A_local_crs->getGraph()->getExporter()));
941 local_matrix_device_type U_localCrsMatrix_device;
942 U_localCrsMatrix_device = local_matrix_device_type(
"U_Factor_localmatrix",
943 A_local_->getLocalNumRows(),
948 A_local_crs->getRowMap(),
949 A_local_crs->getColMap(),
950 A_local_crs->getDomainMap(),
951 A_local_crs->getRangeMap(),
952 A_local_crs->getGraph()->getImporter(),
953 A_local_crs->getGraph()->getExporter()));
955 L_solver_->setMatrix(L_);
956 L_solver_->compute();
957 U_solver_->setMatrix(U_);
958 U_solver_->compute();
962 ComputeTime_ += (timer.
wallTime() - startTime);
967 template <
class MatrixType>
969 apply(
const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
970 Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
976 using Teuchos::rcpFromRef;
979 !isComputed(), std::runtime_error,
980 "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
981 "factorization, before calling apply().");
984 X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
985 "Ifpack2::ILUT::apply: X and Y must have the same number of columns. "
987 << X.getNumVectors() <<
" columns, but Y has "
988 << Y.getNumVectors() <<
" columns.");
994 double startTime = timer.
wallTime();
998 if (alpha == one && beta == zero) {
1001 L_solver_->apply(X, Y, mode);
1004 U_solver_->apply(Y, Y, mode);
1008 U_solver_->apply(X, Y, mode);
1011 L_solver_->apply(Y, Y, mode);
1014 if (alpha == zero) {
1024 MV Y_tmp(Y.getMap(), Y.getNumVectors());
1025 apply(X, Y_tmp, mode);
1026 Y.update(alpha, Y_tmp, beta);
1032 ApplyTime_ += (timer.
wallTime() - startTime);
1035 template <
class MatrixType>
1037 std::ostringstream os;
1042 os <<
"\"Ifpack2::ILUT\": {";
1043 os <<
"Initialized: " << (isInitialized() ?
"true" :
"false") <<
", "
1044 <<
"Computed: " << (isComputed() ?
"true" :
"false") <<
", ";
1046 os <<
"Level-of-fill: " << getLevelOfFill() <<
", "
1047 <<
"absolute threshold: " << getAbsoluteThreshold() <<
", "
1048 <<
"relative threshold: " << getRelativeThreshold() <<
", "
1049 <<
"relaxation value: " << getRelaxValue() <<
", ";
1052 os <<
"Matrix: null";
1054 os <<
"Global matrix dimensions: ["
1055 << A_->getGlobalNumRows() <<
", " << A_->getGlobalNumCols() <<
"]"
1056 <<
", Global nnz: " << A_->getGlobalNumEntries();
1063 template <
class MatrixType>
1084 out <<
"\"Ifpack2::ILUT\":" << endl;
1086 out <<
"MatrixType: " << TypeNameTraits<MatrixType>::name() << endl;
1087 if (this->getObjectLabel() !=
"") {
1088 out <<
"Label: \"" << this->getObjectLabel() <<
"\"" << endl;
1090 out <<
"Initialized: " << (isInitialized() ?
"true" :
"false")
1092 <<
"Computed: " << (isComputed() ?
"true" :
"false")
1094 <<
"Level of fill: " << getLevelOfFill() << endl
1095 <<
"Absolute threshold: " << getAbsoluteThreshold() << endl
1096 <<
"Relative threshold: " << getRelativeThreshold() << endl
1097 <<
"Relax value: " << getRelaxValue() << endl;
1100 const double fillFraction =
1101 (double)getGlobalNumEntries() / (double)A_->getGlobalNumEntries();
1102 const double nnzToRows =
1103 (double)getGlobalNumEntries() / (double)U_->getGlobalNumRows();
1105 out <<
"Dimensions of L: [" << L_->getGlobalNumRows() <<
", "
1106 << L_->getGlobalNumRows() <<
"]" << endl
1107 <<
"Dimensions of U: [" << U_->getGlobalNumRows() <<
", "
1108 << U_->getGlobalNumRows() <<
"]" << endl
1109 <<
"Number of nonzeros in factors: " << getGlobalNumEntries() << endl
1110 <<
"Fill fraction of factors over A: " << fillFraction << endl
1111 <<
"Ratio of nonzeros to rows: " << nnzToRows << endl;
1114 out <<
"Number of initialize calls: " << getNumInitialize() << endl
1115 <<
"Number of compute calls: " << getNumCompute() << endl
1116 <<
"Number of apply calls: " << getNumApply() << endl
1117 <<
"Total time in seconds for initialize: " << getInitializeTime() << endl
1118 <<
"Total time in seconds for compute: " << getComputeTime() << endl
1119 <<
"Total time in seconds for apply: " << getApplyTime() << endl;
1121 out <<
"Local matrix:" << endl;
1122 A_local_->describe(out, vl);
1133 #define IFPACK2_ILUT_INSTANT(S, LO, GO, N) \
1134 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:97
bool hasTransposeApply() const
Whether this object's apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_ILUT_def.hpp:302
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:348
void initialize()
Clear any previously computed factors, and potentially compute sparsity patterns of factors...
Definition: Ifpack2_ILUT_def.hpp:429
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:1065
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_ILUT_def.hpp:1036
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:60
void rm_heap_root(Teuchos::Array< Ordinal > &heap, SizeType &heap_len)
Definition: Ifpack2_Heap.hpp:59
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:292
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_ILUT_decl.hpp:73
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:525
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:358
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:307
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_ILUT_def.hpp:332
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the input matrix's communicator.
Definition: Ifpack2_ILUT_def.hpp:264
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:102
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:281
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:35
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_ILUT_def.hpp:337
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_ILUT_def.hpp:322
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_ILUT_def.hpp:317
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:969
size_t getLocalNumEntries() const
Returns the number of nonzero entries in the local graph.
Definition: Ifpack2_ILUT_def.hpp:353
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_ILUT_def.hpp:327
void setParameters(const Teuchos::ParameterList ¶ms)
Set preconditioner parameters.
Definition: Ifpack2_ILUT_def.hpp:128
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_ILUT_def.hpp:312
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:76
Teuchos::RCP< const row_matrix_type > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:275