Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_ILUT_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4 //
5 // Copyright 2009 NTESS and the Ifpack2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef IFPACK2_ILUT_DEF_HPP
11 #define IFPACK2_ILUT_DEF_HPP
12 
13 #include <type_traits>
15 #include "Teuchos_StandardParameterEntryValidators.hpp"
16 #include "Teuchos_Time.hpp"
17 #include "Tpetra_CrsMatrix.hpp"
18 #include "KokkosSparse_par_ilut.hpp"
19 
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"
25 
26 namespace Ifpack2 {
27 
28 namespace {
29 
30 struct IlutImplType {
31  enum Enum {
32  Serial,
33  PAR_ILUT
34  };
35 
36  static void loadPLTypeOption(Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
37  type_strs.resize(2);
38  type_strs[0] = "serial";
39  type_strs[1] = "par_ilut";
40  type_enums.resize(2);
41  type_enums[0] = Serial;
42  type_enums[1] = PAR_ILUT;
43  }
44 };
45 
70 template <class ScalarType>
72 ilutDefaultDropTolerance() {
74  typedef typename STS::magnitudeType magnitude_type;
76 
77  // 1/2. Hopefully this can be represented in magnitude_type.
78  const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
79 
80  // The min ensures that in case magnitude_type has very low
81  // precision, we'll at least get some value strictly less than
82  // one.
83  return std::min(static_cast<magnitude_type>(1000) * STS::magnitude(STS::eps()), oneHalf);
84 }
85 
86 // Full specialization for ScalarType = double.
87 // This specialization preserves ILUT's previous default behavior.
88 template <>
90 ilutDefaultDropTolerance<double>() {
91  return 1e-12;
92 }
93 
94 } // namespace
95 
96 template <class MatrixType>
98  : A_(A)
99  , Athresh_(Teuchos::ScalarTraits<magnitude_type>::zero())
100  , Rthresh_(Teuchos::ScalarTraits<magnitude_type>::one())
101  , RelaxValue_(Teuchos::ScalarTraits<magnitude_type>::zero())
102  , LevelOfFill_(1.0)
103  , DropTolerance_(ilutDefaultDropTolerance<scalar_type>())
104  , par_ilut_options_{1, 0., -1, -1, 0.75, false}
105  , InitializeTime_(0.0)
106  , ComputeTime_(0.0)
107  , ApplyTime_(0.0)
108  , NumInitialize_(0)
109  , NumCompute_(0)
110  , NumApply_(0)
111  , IsInitialized_(false)
112  , IsComputed_(false)
113  , useKokkosKernelsParILUT_(false)
114 
115 {
116  allocateSolvers();
117 }
118 
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");
125 }
126 
127 template <class MatrixType>
129  using Ifpack2::Details::getParamTryingTypes;
130  const char prefix[] = "Ifpack2::ILUT: ";
131 
132  // Don't actually change the instance variables until we've checked
133  // all parameters. This ensures that setParameters satisfies the
134  // strong exception guarantee (i.e., is transactional).
135 
136  // Parsing implementation type
137  IlutImplType::Enum ilutimplType = IlutImplType::Serial;
138  do {
139  static const char typeName[] = "fact: type";
140 
141  if (!params.isType<std::string>(typeName)) break;
142 
143  // Map std::string <-> IlutImplType::Enum.
144  Teuchos::Array<std::string> ilutimplTypeStrs;
145  Teuchos::Array<IlutImplType::Enum> ilutimplTypeEnums;
146  IlutImplType::loadPLTypeOption(ilutimplTypeStrs, ilutimplTypeEnums);
148  s2i(ilutimplTypeStrs(), ilutimplTypeEnums(), typeName, false);
149 
150  ilutimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
151  } while (0);
152 
153  if (ilutimplType == IlutImplType::PAR_ILUT) {
154  this->useKokkosKernelsParILUT_ = true;
155  } else {
156  this->useKokkosKernelsParILUT_ = false;
157  }
158 
159  // Fill level in ILUT is a double, not a magnitude_type, because it
160  // depends on LO and GO, not on Scalar. Also, you can't cast
161  // arbitrary magnitude_type (e.g., Sacado::MP::Vector) to double.
162  double fillLevel = LevelOfFill_;
163  {
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);
169  TEUCHOS_TEST_FOR_EXCEPTION(fillLevel < 1.0, std::runtime_error,
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 "
177  "fill-in.");
178  }
179 
180  magnitude_type absThresh = Athresh_;
181  {
182  const std::string paramName("fact: absolute threshold");
183  getParamTryingTypes<magnitude_type, magnitude_type, double>(absThresh, params, paramName, prefix);
184  }
185 
186  magnitude_type relThresh = Rthresh_;
187  {
188  const std::string paramName("fact: relative threshold");
189  getParamTryingTypes<magnitude_type, magnitude_type, double>(relThresh, params, paramName, prefix);
190  }
191 
192  magnitude_type relaxValue = RelaxValue_;
193  {
194  const std::string paramName("fact: relax value");
195  getParamTryingTypes<magnitude_type, magnitude_type, double>(relaxValue, params, paramName, prefix);
196  }
197 
198  magnitude_type dropTol = DropTolerance_;
199  {
200  const std::string paramName("fact: drop tolerance");
201  getParamTryingTypes<magnitude_type, magnitude_type, double>(dropTol, params, paramName, prefix);
202  }
203 
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;
217 
218  std::string par_ilut_plist_name("parallel ILUT options");
219  if (params.isSublist(par_ilut_plist_name)) {
220  Teuchos::ParameterList const& par_ilut_plist = params.sublist(par_ilut_plist_name);
221 
222  std::string paramName("maximum iterations");
223  getParamTryingTypes<int, int>(par_ilut_max_iter, par_ilut_plist, paramName, prefix);
224 
225  paramName = "residual norm delta stop";
226  getParamTryingTypes<magnitude_type, magnitude_type, double>(par_ilut_residual_norm_delta_stop, par_ilut_plist, paramName, prefix);
227 
228  paramName = "team size";
229  getParamTryingTypes<int, int>(par_ilut_team_size, par_ilut_plist, paramName, prefix);
230 
231  paramName = "vector size";
232  getParamTryingTypes<int, int>(par_ilut_vector_size, par_ilut_plist, paramName, prefix);
233 
234  paramName = "fill in limit";
235  getParamTryingTypes<float, float, double>(par_ilut_fill_in_limit, par_ilut_plist, paramName, prefix);
236 
237  paramName = "verbose";
238  getParamTryingTypes<bool, bool>(par_ilut_verbose, par_ilut_plist, paramName, prefix);
239 
240  } // if (params.isSublist(par_ilut_plist_name))
241 
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;
248 
249  } // if (this->useKokkosKernelsParILUT_)
250 
251  // Forward to trisolvers.
252  L_solver_->setParameters(params);
253  U_solver_->setParameters(params);
254 
255  LevelOfFill_ = fillLevel;
256  Athresh_ = absThresh;
257  Rthresh_ = relThresh;
258  RelaxValue_ = relaxValue;
259  DropTolerance_ = dropTol;
260 }
261 
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();
271 }
272 
273 template <class MatrixType>
276  return A_;
277 }
278 
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();
288 }
289 
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();
299 }
300 
301 template <class MatrixType>
303  return true;
304 }
305 
306 template <class MatrixType>
308  return NumInitialize_;
309 }
310 
311 template <class MatrixType>
313  return NumCompute_;
314 }
315 
316 template <class MatrixType>
318  return NumApply_;
319 }
320 
321 template <class MatrixType>
323  return InitializeTime_;
324 }
325 
326 template <class MatrixType>
328  return ComputeTime_;
329 }
330 
331 template <class MatrixType>
333  return ApplyTime_;
334 }
335 
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.");
343  // ILUT methods cost roughly one apply + the nnz in the upper+lower triangles
344  return A_->getLocalNumEntries() + getLocalNumEntries();
345 }
346 
347 template <class MatrixType>
349  return L_->getGlobalNumEntries() + U_->getGlobalNumEntries();
350 }
351 
352 template <class MatrixType>
354  return L_->getLocalNumEntries() + U_->getLocalNumEntries();
355 }
356 
357 template <class MatrixType>
359  if (A.getRawPtr() != A_.getRawPtr()) {
360  // Check in serial or one-process mode if the matrix is square.
362  !A.is_null() && A->getComm()->getSize() == 1 &&
363  A->getLocalNumRows() != A->getLocalNumCols(),
364  std::runtime_error,
365  "Ifpack2::ILUT::setMatrix: If A's communicator only "
366  "contains one process, then A must be square. Instead, you provided a "
367  "matrix A with "
368  << A->getLocalNumRows() << " rows and "
369  << A->getLocalNumCols() << " columns.");
370 
371  // It's legal for A to be null; in that case, you may not call
372  // initialize() until calling setMatrix() with a nonnull input.
373  // Regardless, setting the matrix invalidates any previous
374  // factorization.
375  IsInitialized_ = false;
376  IsComputed_ = false;
377  A_local_ = Teuchos::null;
378 
379  // The sparse triangular solvers get a triangular factor as their
380  // input matrix. The triangular factors L_ and U_ are getting
381  // reset, so we reset the solvers' matrices to null. Do that
382  // before setting L_ and U_ to null, so that latter step actually
383  // frees the factors.
384  if (!L_solver_.is_null()) {
385  L_solver_->setMatrix(Teuchos::null);
386  }
387  if (!U_solver_.is_null()) {
388  U_solver_->setMatrix(Teuchos::null);
389  }
390 
391  L_ = Teuchos::null;
392  U_ = Teuchos::null;
393  A_ = A;
394  }
395 }
396 
397 template <class MatrixType>
400  using Teuchos::RCP;
401  using Teuchos::rcp;
402  using Teuchos::rcp_dynamic_cast;
403  using Teuchos::rcp_implicit_cast;
404 
405  // If A_'s communicator only has one process, or if its column and
406  // row Maps are the same, then it is already local, so use it
407  // directly.
408  if (A->getRowMap()->getComm()->getSize() == 1 ||
409  A->getRowMap()->isSameAs(*(A->getColMap()))) {
410  return A;
411  }
412 
413  // If A_ is already a LocalFilter, then use it directly. This
414  // should be the case if RILUT is being used through
415  // AdditiveSchwarz, for example.
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);
420  } else {
421  // A_'s communicator has more than one process, its row Map and
422  // its column Map differ, and A_ is not a LocalFilter. Thus, we
423  // have to wrap it in a LocalFilter.
424  return rcp(new LocalFilter<row_matrix_type>(A));
425  }
426 }
427 
428 template <class MatrixType>
430  using Teuchos::Array;
431  using Teuchos::RCP;
432  using Teuchos::rcp_const_cast;
433  Teuchos::Time timer("ILUT::initialize");
434  double startTime = timer.wallTime();
435  {
436  Teuchos::TimeMonitor timeMon(timer);
437 
438  // Check that the matrix is nonnull.
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.");
444 
445  // Clear any previous computations.
446  IsInitialized_ = false;
447  IsComputed_ = false;
448  A_local_ = Teuchos::null;
449  L_ = Teuchos::null;
450  U_ = Teuchos::null;
451 
452  A_local_ = makeLocalFilter(A_); // Compute the local filter.
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.");
458 
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);
469 
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()) {
472  // the result is a host-based matrix, which is the same as what happens in RILUK
473  local_ordinal_type numRows = A_local_->getLocalNumRows();
474  Array<size_t> entriesPerRow(numRows);
475  for (local_ordinal_type i = 0; i < numRows; i++) {
476  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
477  }
478  RCP<crs_matrix_type> A_local_crs_nc =
479  rcp(new crs_matrix_type(A_local_->getRowMap(),
480  A_local_->getColMap(),
481  entriesPerRow()));
482  // copy entries into A_local_crs
483  nonconst_local_inds_host_view_type indices("indices", A_local_->getLocalMaxNumRowEntries());
484  nonconst_values_host_view_type values("values", A_local_->getLocalMaxNumRowEntries());
485  for (local_ordinal_type i = 0; i < numRows; i++) {
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());
489  }
490  A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
491  A_local_crs = rcp_const_cast<const crs_matrix_type>(A_local_crs_nc);
492  }
493  auto A_local_crs_device = A_local_crs->getLocalMatrixDevice();
494 
495  // KokkosKernels requires unsigned
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);
502 
503  KokkosSparse::Experimental::par_ilut_symbolic(KernelHandle_.getRawPtr(),
504  A_local_crs_device.graph.row_map, A_local_crs_device.graph.entries,
505  L_rowmap_,
506  U_rowmap_);
507 
508  Kokkos::deep_copy(L_rowmap_orig_, L_rowmap_);
509  Kokkos::deep_copy(U_rowmap_orig_, U_rowmap_);
510  }
511 
512  IsInitialized_ = true;
513  ++NumInitialize_;
514  } // timer scope
515  InitializeTime_ += (timer.wallTime() - startTime);
516 }
517 
518 template <typename ScalarType>
520 scalar_mag(const ScalarType& s) {
522 }
523 
524 template <class MatrixType>
526  using Teuchos::Array;
527  using Teuchos::ArrayRCP;
528  using Teuchos::ArrayView;
529  using Teuchos::as;
530  using Teuchos::rcp;
531  using Teuchos::RCP;
532  using Teuchos::rcp_const_cast;
533  using Teuchos::reduceAll;
534 
535  // Don't count initialization in the compute() time.
536  if (!isInitialized()) {
537  initialize();
538  }
539 
540  Teuchos::Time timer("ILUT::compute");
541  double startTime = timer.wallTime();
542  { // Timer scope for timing compute()
543  Teuchos::TimeMonitor timeMon(timer, true);
544 
545  if (!this->useKokkosKernelsParILUT_) {
546  //--------------------------------------------------------------------------
547  // Ifpack2::ILUT's serial version is a translation of the Aztec ILUT
548  // implementation. The Aztec ILUT implementation was written by Ray Tuminaro.
549  //
550  // This isn't an exact translation of the Aztec ILUT algorithm, for the
551  // following reasons:
552  // 1. Minor differences result from the fact that Aztec factors a MSR format
553  // matrix in place, while the code below factors an input CrsMatrix which
554  // remains untouched and stores the resulting factors in separate L and U
555  // CrsMatrix objects.
556  // Also, the Aztec code begins by shifting the matrix pointers back
557  // by one, and the pointer contents back by one, and then using 1-based
558  // Fortran-style indexing in the algorithm. This Ifpack2 code uses C-style
559  // 0-based indexing throughout.
560  // 2. Aztec stores the inverse of the diagonal of U. This Ifpack2 code
561  // stores the non-inverted diagonal in U.
562  // The triangular solves (in Ifpack2::ILUT::apply()) are performed by
563  // calling the Tpetra::CrsMatrix::solve method on the L and U objects, and
564  // this requires U to contain the non-inverted diagonal.
565  //
566  // ABW.
567  //--------------------------------------------------------------------------
568 
569  const scalar_type zero = STS::zero();
570  const scalar_type one = STS::one();
571 
572  const local_ordinal_type myNumRows = A_local_->getLocalNumRows();
573 
574  // If this macro is defined, files containing the L and U factors
575  // will be written. DON'T CHECK IN THE CODE WITH THIS MACRO ENABLED!!!
576  // #define IFPACK2_WRITE_ILUT_FACTORS
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);
580 #endif
581 
582  // Calculate how much fill will be allowed in addition to the
583  // space that corresponds to the input matrix entries.
584  double local_nnz = static_cast<double>(A_local_->getLocalNumEntries());
585  double fill = ((getLevelOfFill() - 1.0) * local_nnz) / (2 * myNumRows);
586 
587  // std::ceil gives the smallest integer larger than the argument.
588  // this may give a slightly different result than Aztec's fill value in
589  // some cases.
590  double fill_ceil = std::ceil(fill);
591 
592  // Similarly to Aztec, we will allow the same amount of fill for each
593  // row, half in L and half in U.
594  size_type fillL = static_cast<size_type>(fill_ceil);
595  size_type fillU = static_cast<size_type>(fill_ceil);
596 
597  Array<scalar_type> InvDiagU(myNumRows, zero);
598 
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);
603 
604  enum { UNUSED,
605  ORIG,
606  FILL };
607  local_ordinal_type max_col = myNumRows;
608 
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;
617 
618  // A comparison object which will be used to create 'heaps' of indices
619  // that are ordered according to the corresponding values in the
620  // 'cur_row' array.
621  greater_indirect<scalar_type, local_ordinal_type> vals_comp(cur_row);
622 
623  // =================== //
624  // start factorization //
625  // =================== //
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);
632  }
633 
634  for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
635  local_inds_host_view_type ColIndicesA;
636  values_host_view_type ColValuesA;
637  size_t RowNnz;
638 
639  if (A_local_->supportsRowViews()) {
640  A_local_->getLocalRowView(row_i, ColIndicesA, ColValuesA);
641  RowNnz = ColIndicesA.size();
642  } else {
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));
646  }
647 
648  // Always include the diagonal in the U factor. The value should get
649  // set in the next loop below.
650  U_cols.push_back(row_i);
651  cur_row[row_i] = zero;
652  pattern[row_i] = ORIG;
653 
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]);
662  }
663 
664  cur_row[ColIndicesA[i]] = ColValuesA[i];
665  pattern[ColIndicesA[i]] = ORIG;
666  rownorm += scalar_mag(ColValuesA[i]);
667  }
668  }
669 
670  // Alter the diagonal according to the absolute-threshold and
671  // relative-threshold values. If not set, those values default
672  // to zero and one respectively.
673  const magnitude_type rthresh = getRelativeThreshold();
674  const scalar_type& v = cur_row[row_i];
675  cur_row[row_i] = as<scalar_type>(getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh * v;
676 
677  size_type orig_U_len = U_cols.size();
678  RowNnz = L_cols_heap.size() + orig_U_len;
679  rownorm = getDropTolerance() * rownorm / RowNnz;
680 
681  // The following while loop corresponds to the 'L30' goto's in Aztec.
682  size_type L_vals_heaplen = 0;
683  while (L_cols_heaplen > 0) {
684  local_ordinal_type row_k = L_cols_heap.front();
685 
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;
691  rm_heap_root(L_cols_heap, L_cols_heaplen);
692  continue;
693  }
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;
700  rm_heap_root(L_cols_heap, L_cols_heaplen);
701  continue;
702  } else {
703  pattern[L_vals_heap.front()] = UNUSED;
704  rm_heap_root(L_vals_heap, L_vals_heaplen, vals_comp);
705  add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
706  }
707  }
708 
709  /* Reduce current row */
710 
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();
714 
715  for (size_type j = 0; j < ColNnzU; ++j) {
716  if (ColIndicesU[j] > row_k) {
717  scalar_type tmp = multiplier * ColValuesU[j];
718  local_ordinal_type col_j = ColIndicesU[j];
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;
724  if (col_j > row_i) {
725  U_cols.push_back(col_j);
726  } else {
727  add_to_heap(col_j, L_cols_heap, L_cols_heaplen);
728  }
729  }
730  }
731  }
732 
733  rm_heap_root(L_cols_heap, L_cols_heaplen);
734  } // end of while(L_cols_heaplen) loop
735 
736  // Put indices and values for L into arrays and then into the L_ matrix.
737 
738  // first, the original entries from the L section of A:
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;
744  }
745  }
746 
747  // next, the L entries resulting from fill:
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;
752  }
753 
754  // L has a one on the diagonal, but we don't explicitly store
755  // it. If we don't store it, then the kernel which performs the
756  // triangular solve can assume a unit diagonal, take a short-cut
757  // and perform faster.
758 
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;
763  }
764 #endif
765 
766  // Pick out the diagonal element, store its reciprocal.
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;
773  }
774  InvDiagU[row_i] = one / cur_row[row_i];
775 
776  // Non-inverted diagonal is stored for U:
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;
781 
782  // Now put indices and values for U into arrays and then into the U_ matrix.
783  // The first entry in U_cols is the diagonal, which we just handled, so we'll
784  // start our loop at j=1.
785 
786  size_type U_vals_heaplen = 0;
787  for (size_type j = 1; j < U_cols.size(); ++j) {
788  local_ordinal_type col = U_cols[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()])) {
794  rm_heap_root(U_vals_heap, U_vals_heaplen, vals_comp);
795  add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
796  }
797  } else {
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]);
801  }
802  pattern[col] = UNUSED;
803  }
804 
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]]);
809  }
810 
811  unorm[row_i] /= (orig_U_len + U_vals_heaplen);
812 
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;
817  }
818 #endif
819 
820  L_cols_heap.clear();
821  U_cols.clear();
822  L_vals_heap.clear();
823  U_vals_heap.clear();
824  } // end of for(row_i) loop
825 
826  // Now allocate and fill the matrices
827  Array<size_t> nnzPerRow(myNumRows);
828 
829  // Make sure to release the old memory for L & U prior to recomputing to
830  // avoid bloating the high-water mark.
831  L_ = Teuchos::null;
832  U_ = Teuchos::null;
833  L_solver_->setMatrix(Teuchos::null);
834  U_solver_->setMatrix(Teuchos::null);
835 
836  for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
837  nnzPerRow[row_i] = L_tmp_idx[row_i].size();
838  }
839 
840  L_ = rcp(new crs_matrix_type(A_local_->getRowMap(), A_local_->getColMap(),
841  nnzPerRow()));
842 
843  for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
844  L_->insertLocalValues(row_i, L_tmp_idx[row_i](), L_tmpv[row_i]());
845  }
846 
847  L_->fillComplete();
848 
849  for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
850  nnzPerRow[row_i] = U_tmp_idx[row_i].size();
851  }
852 
853  U_ = rcp(new crs_matrix_type(A_local_->getRowMap(), A_local_->getColMap(),
854  nnzPerRow()));
855 
856  for (local_ordinal_type row_i = 0; row_i < myNumRows; ++row_i) {
857  U_->insertLocalValues(row_i, U_tmp_idx[row_i](), U_tmpv[row_i]());
858  }
859 
860  U_->fillComplete();
861 
862  L_solver_->setMatrix(L_);
863  L_solver_->initialize();
864  L_solver_->compute();
865 
866  U_solver_->setMatrix(U_);
867  U_solver_->initialize();
868  U_solver_->compute();
869 
870  } // if (!this->useKokkosKernelsParILUT_)
871  else {
872  // Set L, U rowmaps back to original state. Par_ilut can change them, which invalidates them
873  // if compute is called again.
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_);
879  }
880 
881  RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_local_);
882  { // Make sure values in A is picked up even in case of pattern reuse
883  if (A_local_crs.is_null()) {
884  local_ordinal_type numRows = A_local_->getLocalNumRows();
885  Array<size_t> entriesPerRow(numRows);
886  for (local_ordinal_type i = 0; i < numRows; i++) {
887  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
888  }
889  RCP<crs_matrix_type> A_local_crs_nc =
890  rcp(new crs_matrix_type(A_local_->getRowMap(),
891  A_local_->getColMap(),
892  entriesPerRow()));
893  // copy entries into A_local_crs
894  nonconst_local_inds_host_view_type indices("indices", A_local_->getLocalMaxNumRowEntries());
895  nonconst_values_host_view_type values("values", A_local_->getLocalMaxNumRowEntries());
896  for (local_ordinal_type i = 0; i < numRows; i++) {
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());
900  }
901  A_local_crs_nc->fillComplete(A_local_->getDomainMap(), A_local_->getRangeMap());
902  A_local_crs = rcp_const_cast<const crs_matrix_type>(A_local_crs_nc);
903  }
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;
908  }
909 
910  // JHU TODO Should allocation of L & U's column (aka entry) and value arrays occur here or in init()?
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);
915 
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);
919 
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_);
923 
924  auto L_kokkosCrsGraph = local_graph_device_type(L_entries_, L_rowmap_);
925  auto U_kokkosCrsGraph = local_graph_device_type(U_entries_, U_rowmap_);
926 
927  local_matrix_device_type L_localCrsMatrix_device;
928  L_localCrsMatrix_device = local_matrix_device_type("L_Factor_localmatrix",
929  A_local_->getLocalNumRows(),
930  L_values_,
931  L_kokkosCrsGraph);
932 
933  L_ = rcp(new crs_matrix_type(L_localCrsMatrix_device,
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()));
940 
941  local_matrix_device_type U_localCrsMatrix_device;
942  U_localCrsMatrix_device = local_matrix_device_type("U_Factor_localmatrix",
943  A_local_->getLocalNumRows(),
944  U_values_,
945  U_kokkosCrsGraph);
946 
947  U_ = rcp(new crs_matrix_type(U_localCrsMatrix_device,
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()));
954 
955  L_solver_->setMatrix(L_);
956  L_solver_->compute(); // NOTE: Only do compute if the pointer changed. Otherwise, do nothing
957  U_solver_->setMatrix(U_);
958  U_solver_->compute(); // NOTE: Only do compute if the pointer changed. Otherwise, do nothing
959  } // if (!this->useKokkosKernelsParILUT_) ... else ...
960 
961  } // Timer scope for timing compute()
962  ComputeTime_ += (timer.wallTime() - startTime);
963  IsComputed_ = true;
964  ++NumCompute_;
965 } // compute()
966 
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,
971  Teuchos::ETransp mode,
972  scalar_type alpha,
973  scalar_type beta) const {
974  using Teuchos::RCP;
975  using Teuchos::rcp;
976  using Teuchos::rcpFromRef;
977 
979  !isComputed(), std::runtime_error,
980  "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
981  "factorization, before calling apply().");
982 
984  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
985  "Ifpack2::ILUT::apply: X and Y must have the same number of columns. "
986  "X has "
987  << X.getNumVectors() << " columns, but Y has "
988  << Y.getNumVectors() << " columns.");
989 
990  const scalar_type one = STS::one();
991  const scalar_type zero = STS::zero();
992 
993  Teuchos::Time timer("ILUT::apply");
994  double startTime = timer.wallTime();
995  { // Start timing
996  Teuchos::TimeMonitor timeMon(timer, true);
997 
998  if (alpha == one && beta == zero) {
999  if (mode == Teuchos::NO_TRANS) { // Solve L (U Y) = X for Y.
1000  // Start by solving L Y = X for Y.
1001  L_solver_->apply(X, Y, mode);
1002 
1003  // Solve U Y = Y.
1004  U_solver_->apply(Y, Y, mode);
1005  } else { // Solve U^P (L^P Y)) = X for Y (where P is * or T).
1006 
1007  // Start by solving U^P Y = X for Y.
1008  U_solver_->apply(X, Y, mode);
1009 
1010  // Solve L^P Y = Y.
1011  L_solver_->apply(Y, Y, mode);
1012  }
1013  } else { // alpha != 1 or beta != 0
1014  if (alpha == zero) {
1015  // The special case for beta == 0 ensures that if Y contains Inf
1016  // or NaN values, we replace them with 0 (following BLAS
1017  // convention), rather than multiplying them by 0 to get NaN.
1018  if (beta == zero) {
1019  Y.putScalar(zero);
1020  } else {
1021  Y.scale(beta);
1022  }
1023  } else { // alpha != zero
1024  MV Y_tmp(Y.getMap(), Y.getNumVectors());
1025  apply(X, Y_tmp, mode);
1026  Y.update(alpha, Y_tmp, beta);
1027  }
1028  }
1029  } // end timing
1030 
1031  ++NumApply_;
1032  ApplyTime_ += (timer.wallTime() - startTime);
1033 } // apply()
1034 
1035 template <class MatrixType>
1036 std::string ILUT<MatrixType>::description() const {
1037  std::ostringstream os;
1038 
1039  // Output is a valid YAML dictionary in flow style. If you don't
1040  // like everything on a single line, you should call describe()
1041  // instead.
1042  os << "\"Ifpack2::ILUT\": {";
1043  os << "Initialized: " << (isInitialized() ? "true" : "false") << ", "
1044  << "Computed: " << (isComputed() ? "true" : "false") << ", ";
1045 
1046  os << "Level-of-fill: " << getLevelOfFill() << ", "
1047  << "absolute threshold: " << getAbsoluteThreshold() << ", "
1048  << "relative threshold: " << getRelativeThreshold() << ", "
1049  << "relaxation value: " << getRelaxValue() << ", ";
1050 
1051  if (A_.is_null()) {
1052  os << "Matrix: null";
1053  } else {
1054  os << "Global matrix dimensions: ["
1055  << A_->getGlobalNumRows() << ", " << A_->getGlobalNumCols() << "]"
1056  << ", Global nnz: " << A_->getGlobalNumEntries();
1057  }
1058 
1059  os << "}";
1060  return os.str();
1061 }
1062 
1063 template <class MatrixType>
1064 void ILUT<MatrixType>::
1066  const Teuchos::EVerbosityLevel verbLevel) const {
1067  using std::endl;
1068  using Teuchos::Comm;
1069  using Teuchos::OSTab;
1070  using Teuchos::RCP;
1072  using Teuchos::VERB_DEFAULT;
1073  using Teuchos::VERB_EXTREME;
1074  using Teuchos::VERB_HIGH;
1075  using Teuchos::VERB_LOW;
1076  using Teuchos::VERB_MEDIUM;
1077  using Teuchos::VERB_NONE;
1078 
1079  const Teuchos::EVerbosityLevel vl =
1080  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1081  OSTab tab0(out);
1082 
1083  if (vl > VERB_NONE) {
1084  out << "\"Ifpack2::ILUT\":" << endl;
1085  OSTab tab1(out);
1086  out << "MatrixType: " << TypeNameTraits<MatrixType>::name() << endl;
1087  if (this->getObjectLabel() != "") {
1088  out << "Label: \"" << this->getObjectLabel() << "\"" << endl;
1089  }
1090  out << "Initialized: " << (isInitialized() ? "true" : "false")
1091  << endl
1092  << "Computed: " << (isComputed() ? "true" : "false")
1093  << endl
1094  << "Level of fill: " << getLevelOfFill() << endl
1095  << "Absolute threshold: " << getAbsoluteThreshold() << endl
1096  << "Relative threshold: " << getRelativeThreshold() << endl
1097  << "Relax value: " << getRelaxValue() << endl;
1098 
1099  if (isComputed() && vl >= VERB_HIGH) {
1100  const double fillFraction =
1101  (double)getGlobalNumEntries() / (double)A_->getGlobalNumEntries();
1102  const double nnzToRows =
1103  (double)getGlobalNumEntries() / (double)U_->getGlobalNumRows();
1104 
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;
1112  }
1113 
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;
1120 
1121  out << "Local matrix:" << endl;
1122  A_local_->describe(out, vl);
1123  }
1124 }
1125 
1126 } // namespace Ifpack2
1127 
1128 // FIXME (mfh 16 Sep 2014) We should really only use RowMatrix here!
1129 // There's no need to instantiate for CrsMatrix too. All Ifpack2
1130 // preconditioners can and should do dynamic casts if they need a type
1131 // more specific than RowMatrix.
1132 
1133 #define IFPACK2_ILUT_INSTANT(S, LO, GO, N) \
1134  template class Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> >;
1135 
1136 #endif /* IFPACK2_ILUT_DEF_HPP */
ILUT(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_ILUT_def.hpp:97
bool hasTransposeApply() const
Whether this object&#39;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
T * getRawPtr() 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 &paramName="", 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&#39;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
static double wallTime()
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 &params)
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
bool is_null() const
Teuchos::RCP< const row_matrix_type > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:275