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 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_ILUT_DEF_HPP
44 #define IFPACK2_ILUT_DEF_HPP
45 
46 #include <type_traits>
47 #include "Kokkos_StaticCrsGraph.hpp"
49 #include "Teuchos_StandardParameterEntryValidators.hpp"
50 #include "Teuchos_Time.hpp"
51 #include "Tpetra_CrsMatrix.hpp"
52 #include "KokkosSparse_par_ilut.hpp"
53 
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"
59 
60 namespace Ifpack2 {
61 
62  namespace {
63 
64  struct IlutImplType {
65  enum Enum {
66  Serial,
67  PAR_ILUT
68  };
69 
70  static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
71  type_strs.resize(2);
72  type_strs[0] = "serial";
73  type_strs[1] = "par_ilut";
74  type_enums.resize(2);
75  type_enums[0] = Serial;
76  type_enums[1] = PAR_ILUT;
77  }
78  };
79 
80 
105  template<class ScalarType>
107  ilutDefaultDropTolerance () {
109  typedef typename STS::magnitudeType magnitude_type;
111 
112  // 1/2. Hopefully this can be represented in magnitude_type.
113  const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
114 
115  // The min ensures that in case magnitude_type has very low
116  // precision, we'll at least get some value strictly less than
117  // one.
118  return std::min (static_cast<magnitude_type> (1000) * STS::magnitude (STS::eps ()), oneHalf);
119  }
120 
121  // Full specialization for ScalarType = double.
122  // This specialization preserves ILUT's previous default behavior.
123  template<>
125  ilutDefaultDropTolerance<double> () {
126  return 1e-12;
127  }
128 
129  } // namespace (anonymous)
130 
131 
132 template <class MatrixType>
134  A_ (A),
135  Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
136  Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
137  RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
138  LevelOfFill_ (1.0),
139  DropTolerance_ (ilutDefaultDropTolerance<scalar_type> ()),
140  par_ilut_options_{1, 0., -1, -1, 0.75, false},
141  InitializeTime_ (0.0),
142  ComputeTime_ (0.0),
143  ApplyTime_ (0.0),
144  NumInitialize_ (0),
145  NumCompute_ (0),
146  NumApply_ (0),
147  IsInitialized_ (false),
148  IsComputed_ (false),
149  useKokkosKernelsParILUT_(false)
150 
151 {
152  allocateSolvers();
153 }
154 
155 template<class MatrixType>
156 void ILUT<MatrixType>::allocateSolvers ()
157 {
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");
162 }
163 
164 template <class MatrixType>
166 {
167  using Ifpack2::Details::getParamTryingTypes;
168  const char prefix[] = "Ifpack2::ILUT: ";
169 
170  // Don't actually change the instance variables until we've checked
171  // all parameters. This ensures that setParameters satisfies the
172  // strong exception guarantee (i.e., is transactional).
173 
174  // Parsing implementation type
175  IlutImplType::Enum ilutimplType = IlutImplType::Serial;
176  do {
177  static const char typeName[] = "fact: type";
178 
179  if ( ! params.isType<std::string>(typeName)) break;
180 
181  // Map std::string <-> IlutImplType::Enum.
182  Teuchos::Array<std::string> ilutimplTypeStrs;
183  Teuchos::Array<IlutImplType::Enum> ilutimplTypeEnums;
184  IlutImplType::loadPLTypeOption (ilutimplTypeStrs, ilutimplTypeEnums);
186  s2i(ilutimplTypeStrs (), ilutimplTypeEnums (), typeName, false);
187 
188  ilutimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
189  } while (0);
190 
191  if (ilutimplType == IlutImplType::PAR_ILUT) {
192  this->useKokkosKernelsParILUT_ = true;
193  }
194  else {
195  this->useKokkosKernelsParILUT_ = false;
196  }
197 
198  // Fill level in ILUT is a double, not a magnitude_type, because it
199  // depends on LO and GO, not on Scalar. Also, you can't cast
200  // arbitrary magnitude_type (e.g., Sacado::MP::Vector) to double.
201  double fillLevel = LevelOfFill_;
202  {
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 "
217  "fill-in.");
218  }
219 
220  magnitude_type absThresh = Athresh_;
221  {
222  const std::string paramName ("fact: absolute threshold");
223  getParamTryingTypes<magnitude_type, magnitude_type, double>
224  (absThresh, params, paramName, prefix);
225  }
226 
227  magnitude_type relThresh = Rthresh_;
228  {
229  const std::string paramName ("fact: relative threshold");
230  getParamTryingTypes<magnitude_type, magnitude_type, double>
231  (relThresh, params, paramName, prefix);
232  }
233 
234  magnitude_type relaxValue = RelaxValue_;
235  {
236  const std::string paramName ("fact: relax value");
237  getParamTryingTypes<magnitude_type, magnitude_type, double>
238  (relaxValue, params, paramName, prefix);
239  }
240 
241  magnitude_type dropTol = DropTolerance_;
242  {
243  const std::string paramName ("fact: drop tolerance");
244  getParamTryingTypes<magnitude_type, magnitude_type, double>
245  (dropTol, params, paramName, prefix);
246  }
247 
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;
261 
262  std::string par_ilut_plist_name("parallel ILUT options");
263  if (params.isSublist(par_ilut_plist_name)) {
264  Teuchos::ParameterList const &par_ilut_plist = params.sublist(par_ilut_plist_name);
265 
266  std::string paramName("maximum iterations");
267  getParamTryingTypes<int, int>(par_ilut_max_iter, par_ilut_plist, paramName, prefix);
268 
269  paramName = "residual norm delta stop";
270  getParamTryingTypes<magnitude_type, magnitude_type, double>(par_ilut_residual_norm_delta_stop, par_ilut_plist, paramName, prefix);
271 
272  paramName = "team size";
273  getParamTryingTypes<int, int>(par_ilut_team_size, par_ilut_plist, paramName, prefix);
274 
275  paramName = "vector size";
276  getParamTryingTypes<int, int>(par_ilut_vector_size, par_ilut_plist, paramName, prefix);
277 
278  paramName = "fill in limit";
279  getParamTryingTypes<float, float, double>(par_ilut_fill_in_limit, par_ilut_plist, paramName, prefix);
280 
281  paramName = "verbose";
282  getParamTryingTypes<bool, bool>(par_ilut_verbose, par_ilut_plist, paramName, prefix);
283 
284  } // if (params.isSublist(par_ilut_plist_name))
285 
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;
292 
293  } //if (this->useKokkosKernelsParILUT_)
294 
295  // Forward to trisolvers.
296  L_solver_->setParameters(params);
297  U_solver_->setParameters(params);
298 
299  LevelOfFill_ = fillLevel;
300  Athresh_ = absThresh;
301  Rthresh_ = relThresh;
302  RelaxValue_ = relaxValue;
303  DropTolerance_ = dropTol;
304 }
305 
306 
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 ();
315 }
316 
317 
318 template <class MatrixType>
321  return A_;
322 }
323 
324 
325 template <class MatrixType>
328 {
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 ();
334 }
335 
336 
337 template <class MatrixType>
340 {
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 ();
346 }
347 
348 
349 template <class MatrixType>
351  return true;
352 }
353 
354 
355 template <class MatrixType>
357  return NumInitialize_;
358 }
359 
360 
361 template <class MatrixType>
363  return NumCompute_;
364 }
365 
366 
367 template <class MatrixType>
369  return NumApply_;
370 }
371 
372 
373 template <class MatrixType>
375  return InitializeTime_;
376 }
377 
378 
379 template<class MatrixType>
381  return ComputeTime_;
382 }
383 
384 
385 template<class MatrixType>
387  return ApplyTime_;
388 }
389 
390 
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.");
397  // ILUT methods cost roughly one apply + the nnz in the upper+lower triangles
398  return A_->getLocalNumEntries() + getLocalNumEntries();
399 }
400 
401 
402 template<class MatrixType>
404  return L_->getGlobalNumEntries () + U_->getGlobalNumEntries ();
405 }
406 
407 
408 template<class MatrixType>
410  return L_->getLocalNumEntries () + U_->getLocalNumEntries ();
411 }
412 
413 
414 template<class MatrixType>
416 {
417  if (A.getRawPtr () != A_.getRawPtr ()) {
418  // Check in serial or one-process mode if the matrix is square.
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.");
426 
427  // It's legal for A to be null; in that case, you may not call
428  // initialize() until calling setMatrix() with a nonnull input.
429  // Regardless, setting the matrix invalidates any previous
430  // factorization.
431  IsInitialized_ = false;
432  IsComputed_ = false;
433  A_local_ = Teuchos::null;
434 
435  // The sparse triangular solvers get a triangular factor as their
436  // input matrix. The triangular factors L_ and U_ are getting
437  // reset, so we reset the solvers' matrices to null. Do that
438  // before setting L_ and U_ to null, so that latter step actually
439  // frees the factors.
440  if (! L_solver_.is_null ()) {
441  L_solver_->setMatrix (Teuchos::null);
442  }
443  if (! U_solver_.is_null ()) {
444  U_solver_->setMatrix (Teuchos::null);
445  }
446 
447  L_ = Teuchos::null;
448  U_ = Teuchos::null;
449  A_ = A;
450  }
451 }
452 
453 template <class MatrixType>
456 {
457  using Teuchos::RCP;
458  using Teuchos::rcp;
459  using Teuchos::rcp_dynamic_cast;
460  using Teuchos::rcp_implicit_cast;
461 
462  // If A_'s communicator only has one process, or if its column and
463  // row Maps are the same, then it is already local, so use it
464  // directly.
465  if (A->getRowMap ()->getComm ()->getSize () == 1 ||
466  A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
467  return A;
468  }
469 
470  // If A_ is already a LocalFilter, then use it directly. This
471  // should be the case if RILUT is being used through
472  // AdditiveSchwarz, for example.
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);
477  }
478  else {
479  // A_'s communicator has more than one process, its row Map and
480  // its column Map differ, and A_ is not a LocalFilter. Thus, we
481  // have to wrap it in a LocalFilter.
482  return rcp (new LocalFilter<row_matrix_type> (A));
483  }
484 }
485 
486 
487 template<class MatrixType>
489 {
490  using Teuchos::RCP;
491  using Teuchos::Array;
492  using Teuchos::rcp_const_cast;
493  Teuchos::Time timer ("ILUT::initialize");
494  double startTime = timer.wallTime();
495  {
496  Teuchos::TimeMonitor timeMon (timer);
497 
498  // Check that the matrix is nonnull.
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.");
503 
504  // Clear any previous computations.
505  IsInitialized_ = false;
506  IsComputed_ = false;
507  A_local_ = Teuchos::null;
508  L_ = Teuchos::null;
509  U_ = Teuchos::null;
510 
511  A_local_ = makeLocalFilter(A_); // Compute the local filter.
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.");
516 
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);
527 
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()) {
530  // the result is a host-based matrix, which is the same as what happens in RILUK
531  local_ordinal_type numRows = A_local_->getLocalNumRows();
532  Array<size_t> entriesPerRow(numRows);
533  for(local_ordinal_type i = 0; i < numRows; i++) {
534  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
535  }
536  RCP<crs_matrix_type> A_local_crs_nc =
537  rcp (new crs_matrix_type (A_local_->getRowMap (),
538  A_local_->getColMap (),
539  entriesPerRow()));
540  // copy entries into A_local_crs
541  nonconst_local_inds_host_view_type indices("indices",A_local_->getLocalMaxNumRowEntries());
542  nonconst_values_host_view_type values("values",A_local_->getLocalMaxNumRowEntries());
543  for(local_ordinal_type i = 0; i < numRows; i++) {
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());
547  }
548  A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
549  A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
550  }
551  auto A_local_crs_device = A_local_crs->getLocalMatrixDevice();
552 
553  //KokkosKernels requires unsigned
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);
558 
559  KokkosSparse::Experimental::par_ilut_symbolic(KernelHandle_.getRawPtr(),
560  A_local_crs_device.graph.row_map, A_local_crs_device.graph.entries,
561  L_rowmap_,
562  U_rowmap_);
563  }
564 
565  IsInitialized_ = true;
566  ++NumInitialize_;
567  } //timer scope
568  InitializeTime_ += (timer.wallTime() - startTime);
569 }
570 
571 
572 template<typename ScalarType>
574 scalar_mag (const ScalarType& s)
575 {
577 }
578 
579 
580 template<class MatrixType>
582 {
583  using Teuchos::Array;
584  using Teuchos::ArrayRCP;
585  using Teuchos::ArrayView;
586  using Teuchos::as;
587  using Teuchos::rcp;
588  using Teuchos::reduceAll;
589  using Teuchos::RCP;
590  using Teuchos::rcp_const_cast;
591 
592  // Don't count initialization in the compute() time.
593  if (! isInitialized ()) {
594  initialize ();
595  }
596 
597  Teuchos::Time timer ("ILUT::compute");
598  double startTime = timer.wallTime();
599  { // Timer scope for timing compute()
600  Teuchos::TimeMonitor timeMon (timer, true);
601 
602  if (!this->useKokkosKernelsParILUT_)
603  {
604  //--------------------------------------------------------------------------
605  // Ifpack2::ILUT's serial version is a translation of the Aztec ILUT
606  // implementation. The Aztec ILUT implementation was written by Ray Tuminaro.
607  //
608  // This isn't an exact translation of the Aztec ILUT algorithm, for the
609  // following reasons:
610  // 1. Minor differences result from the fact that Aztec factors a MSR format
611  // matrix in place, while the code below factors an input CrsMatrix which
612  // remains untouched and stores the resulting factors in separate L and U
613  // CrsMatrix objects.
614  // Also, the Aztec code begins by shifting the matrix pointers back
615  // by one, and the pointer contents back by one, and then using 1-based
616  // Fortran-style indexing in the algorithm. This Ifpack2 code uses C-style
617  // 0-based indexing throughout.
618  // 2. Aztec stores the inverse of the diagonal of U. This Ifpack2 code
619  // stores the non-inverted diagonal in U.
620  // The triangular solves (in Ifpack2::ILUT::apply()) are performed by
621  // calling the Tpetra::CrsMatrix::solve method on the L and U objects, and
622  // this requires U to contain the non-inverted diagonal.
623  //
624  // ABW.
625  //--------------------------------------------------------------------------
626 
627  const scalar_type zero = STS::zero ();
628  const scalar_type one = STS::one ();
629 
630  const local_ordinal_type myNumRows = A_local_->getLocalNumRows ();
631 
632  // If this macro is defined, files containing the L and U factors
633  // will be written. DON'T CHECK IN THE CODE WITH THIS MACRO ENABLED!!!
634  // #define IFPACK2_WRITE_ILUT_FACTORS
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);
638 #endif
639 
640  // Calculate how much fill will be allowed in addition to the
641  // space that corresponds to the input matrix entries.
642  double local_nnz = static_cast<double> (A_local_->getLocalNumEntries ());
643  double fill = ((getLevelOfFill () - 1.0) * local_nnz) / (2 * myNumRows);
644 
645  // std::ceil gives the smallest integer larger than the argument.
646  // this may give a slightly different result than Aztec's fill value in
647  // some cases.
648  double fill_ceil=std::ceil(fill);
649 
650  // Similarly to Aztec, we will allow the same amount of fill for each
651  // row, half in L and half in U.
652  size_type fillL = static_cast<size_type>(fill_ceil);
653  size_type fillU = static_cast<size_type>(fill_ceil);
654 
655  Array<scalar_type> InvDiagU (myNumRows, zero);
656 
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);
661 
662  enum { UNUSED, ORIG, FILL };
663  local_ordinal_type max_col = myNumRows;
664 
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;
673 
674  // A comparison object which will be used to create 'heaps' of indices
675  // that are ordered according to the corresponding values in the
676  // 'cur_row' array.
677  greater_indirect<scalar_type,local_ordinal_type> vals_comp(cur_row);
678 
679  // =================== //
680  // start factorization //
681  // =================== //
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);
688  }
689 
690  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
691  local_inds_host_view_type ColIndicesA;
692  values_host_view_type ColValuesA;
693  size_t RowNnz;
694 
695  if (A_local_->supportsRowViews ()) {
696  A_local_->getLocalRowView (row_i, ColIndicesA, ColValuesA);
697  RowNnz = ColIndicesA.size ();
698  }
699  else {
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));
703  }
704 
705  // Always include the diagonal in the U factor. The value should get
706  // set in the next loop below.
707  U_cols.push_back(row_i);
708  cur_row[row_i] = zero;
709  pattern[row_i] = ORIG;
710 
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);
717  }
718  else if (ColIndicesA[i] > row_i) {
719  U_cols.push_back(ColIndicesA[i]);
720  }
721 
722  cur_row[ColIndicesA[i]] = ColValuesA[i];
723  pattern[ColIndicesA[i]] = ORIG;
724  rownorm += scalar_mag(ColValuesA[i]);
725  }
726  }
727 
728  // Alter the diagonal according to the absolute-threshold and
729  // relative-threshold values. If not set, those values default
730  // to zero and one respectively.
731  const magnitude_type rthresh = getRelativeThreshold();
732  const scalar_type& v = cur_row[row_i];
733  cur_row[row_i] = as<scalar_type> (getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh*v;
734 
735  size_type orig_U_len = U_cols.size();
736  RowNnz = L_cols_heap.size() + orig_U_len;
737  rownorm = getDropTolerance() * rownorm/RowNnz;
738 
739  // The following while loop corresponds to the 'L30' goto's in Aztec.
740  size_type L_vals_heaplen = 0;
741  while (L_cols_heaplen > 0) {
742  local_ordinal_type row_k = L_cols_heap.front();
743 
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;
749  rm_heap_root(L_cols_heap, L_cols_heaplen);
750  continue;
751  }
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);
755  }
756  else if (L_vals_heaplen==0 ||
757  mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
758  pattern[row_k] = UNUSED;
759  rm_heap_root(L_cols_heap, L_cols_heaplen);
760  continue;
761  }
762  else {
763  pattern[L_vals_heap.front()] = UNUSED;
764  rm_heap_root(L_vals_heap, L_vals_heaplen, vals_comp);
765  add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
766  }
767  }
768 
769  /* Reduce current row */
770 
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();
774 
775  for(size_type j=0; j<ColNnzU; ++j) {
776  if (ColIndicesU[j] > row_k) {
777  scalar_type tmp = multiplier * ColValuesU[j];
778  local_ordinal_type col_j = ColIndicesU[j];
779  if (pattern[col_j] != UNUSED) {
780  cur_row[col_j] -= tmp;
781  }
782  else if (scalar_mag(tmp) > rownorm) {
783  cur_row[col_j] = -tmp;
784  pattern[col_j] = FILL;
785  if (col_j > row_i) {
786  U_cols.push_back(col_j);
787  }
788  else {
789  add_to_heap(col_j, L_cols_heap, L_cols_heaplen);
790  }
791  }
792  }
793  }
794 
795  rm_heap_root(L_cols_heap, L_cols_heaplen);
796  }//end of while(L_cols_heaplen) loop
797 
798 
799  // Put indices and values for L into arrays and then into the L_ matrix.
800 
801  // first, the original entries from the L section of A:
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;
807  }
808  }
809 
810  // next, the L entries resulting from fill:
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;
815  }
816 
817  // L has a one on the diagonal, but we don't explicitly store
818  // it. If we don't store it, then the kernel which performs the
819  // triangular solve can assume a unit diagonal, take a short-cut
820  // and perform faster.
821 
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;
826  }
827 #endif
828 
829 
830  // Pick out the diagonal element, store its reciprocal.
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;
837  }
838  InvDiagU[row_i] = one / cur_row[row_i];
839 
840  // Non-inverted diagonal is stored for U:
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;
845 
846  // Now put indices and values for U into arrays and then into the U_ matrix.
847  // The first entry in U_cols is the diagonal, which we just handled, so we'll
848  // start our loop at j=1.
849 
850  size_type U_vals_heaplen = 0;
851  for(size_type j=1; j<U_cols.size(); ++j) {
852  local_ordinal_type col = U_cols[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);
856  }
857  else if (U_vals_heaplen!=0 && scalar_mag(cur_row[col]) >
858  scalar_mag(cur_row[U_vals_heap.front()])) {
859  rm_heap_root(U_vals_heap, U_vals_heaplen, vals_comp);
860  add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
861  }
862  }
863  else {
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]);
867  }
868  pattern[col] = UNUSED;
869  }
870 
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]]);
875  }
876 
877  unorm[row_i] /= (orig_U_len + U_vals_heaplen);
878 
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;
883  }
884 #endif
885 
886  L_cols_heap.clear();
887  U_cols.clear();
888  L_vals_heap.clear();
889  U_vals_heap.clear();
890  } // end of for(row_i) loop
891 
892  // Now allocate and fill the matrices
893  Array<size_t> nnzPerRow(myNumRows);
894 
895  // Make sure to release the old memory for L & U prior to recomputing to
896  // avoid bloating the high-water mark.
897  L_ = Teuchos::null;
898  U_ = Teuchos::null;
899  L_solver_->setMatrix(Teuchos::null);
900  U_solver_->setMatrix(Teuchos::null);
901 
902  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
903  nnzPerRow[row_i] = L_tmp_idx[row_i].size();
904  }
905 
906  L_ = rcp (new crs_matrix_type (A_local_->getRowMap(), A_local_->getColMap(),
907  nnzPerRow()));
908 
909  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
910  L_->insertLocalValues (row_i, L_tmp_idx[row_i](), L_tmpv[row_i]());
911  }
912 
913  L_->fillComplete();
914 
915  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
916  nnzPerRow[row_i] = U_tmp_idx[row_i].size();
917  }
918 
919  U_ = rcp (new crs_matrix_type (A_local_->getRowMap(), A_local_->getColMap(),
920  nnzPerRow()));
921 
922  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
923  U_->insertLocalValues (row_i, U_tmp_idx[row_i](), U_tmpv[row_i]());
924  }
925 
926  U_->fillComplete();
927 
928  L_solver_->setMatrix(L_);
929  L_solver_->initialize ();
930  L_solver_->compute ();
931 
932  U_solver_->setMatrix(U_);
933  U_solver_->initialize ();
934  U_solver_->compute ();
935 
936  } //if (!this->useKokkosKernelsParILUT_)
937  else {
938  RCP<const crs_matrix_type> A_local_crs = Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_local_);
939  {//Make sure values in A is picked up even in case of pattern reuse
940  if(A_local_crs.is_null()) {
941  local_ordinal_type numRows = A_local_->getLocalNumRows();
942  Array<size_t> entriesPerRow(numRows);
943  for(local_ordinal_type i = 0; i < numRows; i++) {
944  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
945  }
946  RCP<crs_matrix_type> A_local_crs_nc =
947  rcp (new crs_matrix_type (A_local_->getRowMap (),
948  A_local_->getColMap (),
949  entriesPerRow()));
950  // copy entries into A_local_crs
951  nonconst_local_inds_host_view_type indices("indices",A_local_->getLocalMaxNumRowEntries());
952  nonconst_values_host_view_type values("values",A_local_->getLocalMaxNumRowEntries());
953  for(local_ordinal_type i = 0; i < numRows; i++) {
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());
957  }
958  A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
959  A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
960  }
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;
965  }
966 
967  //JHU TODO Should allocation of L & U's column (aka entry) and value arrays occur here or in init()?
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);
972 
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);
976 
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_);
980 
981  auto L_kokkosCrsGraph = local_graph_device_type(L_entries_, L_rowmap_);
982  auto U_kokkosCrsGraph = local_graph_device_type(U_entries_, U_rowmap_);
983 
984  local_matrix_device_type L_localCrsMatrix_device;
985  L_localCrsMatrix_device = local_matrix_device_type("L_Factor_localmatrix",
986  A_local_->getLocalNumRows(),
987  L_values_,
988  L_kokkosCrsGraph);
989 
990  L_ = rcp (new crs_matrix_type (L_localCrsMatrix_device,
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()));
997 
998  local_matrix_device_type U_localCrsMatrix_device;
999  U_localCrsMatrix_device = local_matrix_device_type("U_Factor_localmatrix",
1000  A_local_->getLocalNumRows(),
1001  U_values_,
1002  U_kokkosCrsGraph);
1003 
1004  U_ = rcp (new crs_matrix_type (U_localCrsMatrix_device,
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()));
1011 
1012  L_solver_->setMatrix (L_);
1013  L_solver_->compute ();//NOTE: Only do compute if the pointer changed. Otherwise, do nothing
1014  U_solver_->setMatrix (U_);
1015  U_solver_->compute ();//NOTE: Only do compute if the pointer changed. Otherwise, do nothing
1016 
1017  } //if (!this->useKokkosKernelsParILUT_) ... else ...
1018 
1019  } // Timer scope for timing compute()
1020  ComputeTime_ += (timer.wallTime() - startTime);
1021  IsComputed_ = true;
1022  ++NumCompute_;
1023 } //compute()
1024 
1025 
1026 template <class MatrixType>
1027 void ILUT<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,
1030  Teuchos::ETransp mode,
1031  scalar_type alpha,
1032  scalar_type beta) const
1033 {
1034  using Teuchos::RCP;
1035  using Teuchos::rcp;
1036  using Teuchos::rcpFromRef;
1037 
1039  ! isComputed (), std::runtime_error,
1040  "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
1041  "factorization, before calling apply().");
1042 
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.");
1048 
1049  const scalar_type one = STS::one ();
1050  const scalar_type zero = STS::zero ();
1051 
1052  Teuchos::Time timer ("ILUT::apply");
1053  double startTime = timer.wallTime();
1054  { // Start timing
1055  Teuchos::TimeMonitor timeMon (timer, true);
1056 
1057  if (alpha == one && beta == zero) {
1058  if (mode == Teuchos::NO_TRANS) { // Solve L (U Y) = X for Y.
1059  // Start by solving L Y = X for Y.
1060  L_solver_->apply (X, Y, mode);
1061 
1062  // Solve U Y = Y.
1063  U_solver_->apply (Y, Y, mode);
1064  }
1065  else { // Solve U^P (L^P Y)) = X for Y (where P is * or T).
1066 
1067  // Start by solving U^P Y = X for Y.
1068  U_solver_->apply (X, Y, mode);
1069 
1070  // Solve L^P Y = Y.
1071  L_solver_->apply (Y, Y, mode);
1072  }
1073  }
1074  else { // alpha != 1 or beta != 0
1075  if (alpha == zero) {
1076  // The special case for beta == 0 ensures that if Y contains Inf
1077  // or NaN values, we replace them with 0 (following BLAS
1078  // convention), rather than multiplying them by 0 to get NaN.
1079  if (beta == zero) {
1080  Y.putScalar (zero);
1081  } else {
1082  Y.scale (beta);
1083  }
1084  } else { // alpha != zero
1085  MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1086  apply (X, Y_tmp, mode);
1087  Y.update (alpha, Y_tmp, beta);
1088  }
1089  }
1090  }//end timing
1091 
1092  ++NumApply_;
1093  ApplyTime_ += (timer.wallTime() - startTime);
1094 } //apply()
1095 
1096 
1097 template <class MatrixType>
1098 std::string ILUT<MatrixType>::description () const
1099 {
1100  std::ostringstream os;
1101 
1102  // Output is a valid YAML dictionary in flow style. If you don't
1103  // like everything on a single line, you should call describe()
1104  // instead.
1105  os << "\"Ifpack2::ILUT\": {";
1106  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
1107  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
1108 
1109  os << "Level-of-fill: " << getLevelOfFill() << ", "
1110  << "absolute threshold: " << getAbsoluteThreshold() << ", "
1111  << "relative threshold: " << getRelativeThreshold() << ", "
1112  << "relaxation value: " << getRelaxValue() << ", ";
1113 
1114  if (A_.is_null ()) {
1115  os << "Matrix: null";
1116  }
1117  else {
1118  os << "Global matrix dimensions: ["
1119  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]"
1120  << ", Global nnz: " << A_->getGlobalNumEntries();
1121  }
1122 
1123  os << "}";
1124  return os.str ();
1125 }
1126 
1127 
1128 template <class MatrixType>
1129 void
1132  const Teuchos::EVerbosityLevel verbLevel) const
1133 {
1134  using Teuchos::Comm;
1135  using Teuchos::OSTab;
1136  using Teuchos::RCP;
1138  using std::endl;
1139  using Teuchos::VERB_DEFAULT;
1140  using Teuchos::VERB_NONE;
1141  using Teuchos::VERB_LOW;
1142  using Teuchos::VERB_MEDIUM;
1143  using Teuchos::VERB_HIGH;
1144  using Teuchos::VERB_EXTREME;
1145 
1146  const Teuchos::EVerbosityLevel vl =
1147  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1148  OSTab tab0 (out);
1149 
1150  if (vl > VERB_NONE) {
1151  out << "\"Ifpack2::ILUT\":" << endl;
1152  OSTab tab1 (out);
1153  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
1154  if (this->getObjectLabel () != "") {
1155  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
1156  }
1157  out << "Initialized: " << (isInitialized () ? "true" : "false")
1158  << endl
1159  << "Computed: " << (isComputed () ? "true" : "false")
1160  << endl
1161  << "Level of fill: " << getLevelOfFill () << endl
1162  << "Absolute threshold: " << getAbsoluteThreshold () << endl
1163  << "Relative threshold: " << getRelativeThreshold () << endl
1164  << "Relax value: " << getRelaxValue () << endl;
1165 
1166  if (isComputed () && vl >= VERB_HIGH) {
1167  const double fillFraction =
1168  (double) getGlobalNumEntries () / (double) A_->getGlobalNumEntries ();
1169  const double nnzToRows =
1170  (double) getGlobalNumEntries () / (double) U_->getGlobalNumRows ();
1171 
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;
1179  }
1180 
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;
1187 
1188  out << "Local matrix:" << endl;
1189  A_local_->describe (out, vl);
1190  }
1191 }
1192 
1193 }//namespace Ifpack2
1194 
1195 
1196 // FIXME (mfh 16 Sep 2014) We should really only use RowMatrix here!
1197 // There's no need to instantiate for CrsMatrix too. All Ifpack2
1198 // preconditioners can and should do dynamic casts if they need a type
1199 // more specific than RowMatrix.
1200 
1201 #define IFPACK2_ILUT_INSTANT(S,LO,GO,N) \
1202  template class Ifpack2::ILUT< Tpetra::RowMatrix<S, LO, GO, N> >;
1203 
1204 #endif /* IFPACK2_ILUT_DEF_HPP */
ILUT(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_ILUT_def.hpp:133
bool hasTransposeApply() const
Whether this object&#39;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
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: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 &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: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&#39;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
static double wallTime()
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 &params)
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
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:320