Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_RILUK_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 // ***********************************************************************
38 //@HEADER
39 */
40 
41 #ifndef IFPACK2_CRSRILUK_DEF_HPP
42 #define IFPACK2_CRSRILUK_DEF_HPP
43 
44 #include "Ifpack2_LocalFilter.hpp"
45 #include "Tpetra_CrsMatrix.hpp"
46 #include "Teuchos_StandardParameterEntryValidators.hpp"
47 #include "Ifpack2_LocalSparseTriangularSolver.hpp"
48 #include "Ifpack2_Details_getParamTryingTypes.hpp"
49 #include "Kokkos_Sort.hpp"
50 #include "KokkosKernels_SparseUtils.hpp"
51 #include "KokkosKernels_Sorting.hpp"
52 
53 namespace Ifpack2 {
54 
55 namespace Details {
56 struct IlukImplType {
57  enum Enum {
58  Serial,
59  KSPILUK
60  };
61 
62  static void loadPLTypeOption (Teuchos::Array<std::string>& type_strs, Teuchos::Array<Enum>& type_enums) {
63  type_strs.resize(2);
64  type_strs[0] = "Serial";
65  type_strs[1] = "KSPILUK";
66  type_enums.resize(2);
67  type_enums[0] = Serial;
68  type_enums[1] = KSPILUK;
69  }
70 };
71 }
72 
73 template<class MatrixType>
74 RILUK<MatrixType>::RILUK (const Teuchos::RCP<const row_matrix_type>& Matrix_in)
75  : A_ (Matrix_in),
76  LevelOfFill_ (0),
77  Overalloc_ (2.),
78  isAllocated_ (false),
79  isInitialized_ (false),
80  isComputed_ (false),
81  numInitialize_ (0),
82  numCompute_ (0),
83  numApply_ (0),
84  initializeTime_ (0.0),
85  computeTime_ (0.0),
86  applyTime_ (0.0),
87  RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
88  Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
89  Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
90  isKokkosKernelsSpiluk_(false)
91 {
92  allocateSolvers();
93 }
94 
95 
96 template<class MatrixType>
97 RILUK<MatrixType>::RILUK (const Teuchos::RCP<const crs_matrix_type>& Matrix_in)
98  : A_ (Matrix_in),
99  LevelOfFill_ (0),
100  Overalloc_ (2.),
101  isAllocated_ (false),
102  isInitialized_ (false),
103  isComputed_ (false),
104  numInitialize_ (0),
105  numCompute_ (0),
106  numApply_ (0),
107  initializeTime_ (0.0),
108  computeTime_ (0.0),
109  applyTime_ (0.0),
110  RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
111  Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
112  Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
113  isKokkosKernelsSpiluk_(false)
114 {
115  allocateSolvers();
116 }
117 
118 
119 template<class MatrixType>
121 {
122  if (Teuchos::nonnull (KernelHandle_))
123  {
124  KernelHandle_->destroy_spiluk_handle();
125  }
126 }
127 
128 template<class MatrixType>
130 {
131  L_solver_ = Teuchos::rcp (new LocalSparseTriangularSolver<row_matrix_type> ());
132  U_solver_ = Teuchos::rcp (new LocalSparseTriangularSolver<row_matrix_type> ());
133 }
134 
135 template<class MatrixType>
136 void
137 RILUK<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
138 {
139  // It's legal for A to be null; in that case, you may not call
140  // initialize() until calling setMatrix() with a nonnull input.
141  // Regardless, setting the matrix invalidates any previous
142  // factorization.
143  if (A.getRawPtr () != A_.getRawPtr ()) {
144  isAllocated_ = false;
145  isInitialized_ = false;
146  isComputed_ = false;
147  A_local_ = Teuchos::null;
148  Graph_ = Teuchos::null;
149 
150  // The sparse triangular solvers get a triangular factor as their
151  // input matrix. The triangular factors L_ and U_ are getting
152  // reset, so we reset the solvers' matrices to null. Do that
153  // before setting L_ and U_ to null, so that latter step actually
154  // frees the factors.
155  if (! L_solver_.is_null ()) {
156  L_solver_->setMatrix (Teuchos::null);
157  }
158  if (! U_solver_.is_null ()) {
159  U_solver_->setMatrix (Teuchos::null);
160  }
161 
162  L_ = Teuchos::null;
163  U_ = Teuchos::null;
164  D_ = Teuchos::null;
165  A_ = A;
166  }
167 }
168 
169 
170 
171 template<class MatrixType>
174 {
175  TEUCHOS_TEST_FOR_EXCEPTION(
176  L_.is_null (), std::runtime_error, "Ifpack2::RILUK::getL: The L factor "
177  "is null. Please call initialize() (and preferably also compute()) "
178  "before calling this method. If the input matrix has not yet been set, "
179  "you must first call setMatrix() with a nonnull input matrix before you "
180  "may call initialize() or compute().");
181  return *L_;
182 }
183 
184 
185 template<class MatrixType>
186 const Tpetra::Vector<typename RILUK<MatrixType>::scalar_type,
191 {
192  TEUCHOS_TEST_FOR_EXCEPTION(
193  D_.is_null (), std::runtime_error, "Ifpack2::RILUK::getD: The D factor "
194  "(of diagonal entries) is null. Please call initialize() (and "
195  "preferably also compute()) before calling this method. If the input "
196  "matrix has not yet been set, you must first call setMatrix() with a "
197  "nonnull input matrix before you may call initialize() or compute().");
198  return *D_;
199 }
200 
201 
202 template<class MatrixType>
205 {
206  TEUCHOS_TEST_FOR_EXCEPTION(
207  U_.is_null (), std::runtime_error, "Ifpack2::RILUK::getU: The U factor "
208  "is null. Please call initialize() (and preferably also compute()) "
209  "before calling this method. If the input matrix has not yet been set, "
210  "you must first call setMatrix() with a nonnull input matrix before you "
211  "may call initialize() or compute().");
212  return *U_;
213 }
214 
215 
216 template<class MatrixType>
218  TEUCHOS_TEST_FOR_EXCEPTION(
219  A_.is_null (), std::runtime_error, "Ifpack2::RILUK::getNodeSmootherComplexity: "
220  "The input matrix A is null. Please call setMatrix() with a nonnull "
221  "input matrix, then call compute(), before calling this method.");
222  // RILUK methods cost roughly one apply + the nnz in the upper+lower triangles
223  if(!L_.is_null() && !U_.is_null())
224  return A_->getNodeNumEntries() + L_->getNodeNumEntries() + U_->getNodeNumEntries();
225  else
226  return 0;
227 }
228 
229 
230 template<class MatrixType>
231 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
233  typename RILUK<MatrixType>::node_type> >
235  TEUCHOS_TEST_FOR_EXCEPTION(
236  A_.is_null (), std::runtime_error, "Ifpack2::RILUK::getDomainMap: "
237  "The matrix is null. Please call setMatrix() with a nonnull input "
238  "before calling this method.");
239 
240  // FIXME (mfh 25 Jan 2014) Shouldn't this just come from A_?
241  TEUCHOS_TEST_FOR_EXCEPTION(
242  Graph_.is_null (), std::runtime_error, "Ifpack2::RILUK::getDomainMap: "
243  "The computed graph is null. Please call initialize() before calling "
244  "this method.");
245  return Graph_->getL_Graph ()->getDomainMap ();
246 }
247 
248 
249 template<class MatrixType>
250 Teuchos::RCP<const Tpetra::Map<typename RILUK<MatrixType>::local_ordinal_type,
252  typename RILUK<MatrixType>::node_type> >
254  TEUCHOS_TEST_FOR_EXCEPTION(
255  A_.is_null (), std::runtime_error, "Ifpack2::RILUK::getRangeMap: "
256  "The matrix is null. Please call setMatrix() with a nonnull input "
257  "before calling this method.");
258 
259  // FIXME (mfh 25 Jan 2014) Shouldn't this just come from A_?
260  TEUCHOS_TEST_FOR_EXCEPTION(
261  Graph_.is_null (), std::runtime_error, "Ifpack2::RILUK::getRangeMap: "
262  "The computed graph is null. Please call initialize() before calling "
263  "this method.");
264  return Graph_->getL_Graph ()->getRangeMap ();
265 }
266 
267 
268 template<class MatrixType>
270 {
271  using Teuchos::null;
272  using Teuchos::rcp;
273 
274  if (! isAllocated_) {
275  // Deallocate any existing storage. This avoids storing 2x
276  // memory, since RCP op= does not deallocate until after the
277  // assignment.
278  L_ = null;
279  U_ = null;
280  D_ = null;
281 
282  // Allocate Matrix using ILUK graphs
283  L_ = rcp (new crs_matrix_type (Graph_->getL_Graph ()));
284  U_ = rcp (new crs_matrix_type (Graph_->getU_Graph ()));
285  L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
286  U_->setAllToScalar (STS::zero ());
287 
288  // FIXME (mfh 24 Jan 2014) This assumes domain == range Map for L and U.
289  L_->fillComplete ();
290  U_->fillComplete ();
291  D_ = rcp (new vec_type (Graph_->getL_Graph ()->getRowMap ()));
292  }
293  isAllocated_ = true;
294 }
295 
296 
297 template<class MatrixType>
298 void
299 RILUK<MatrixType>::
300 setParameters (const Teuchos::ParameterList& params)
301 {
302  using Teuchos::RCP;
303  using Teuchos::ParameterList;
304  using Teuchos::Array;
305  using Details::getParamTryingTypes;
306  const char prefix[] = "Ifpack2::RILUK: ";
307 
308  // Default values of the various parameters.
309  int fillLevel = 0;
310  magnitude_type absThresh = STM::zero ();
311  magnitude_type relThresh = STM::one ();
312  magnitude_type relaxValue = STM::zero ();
313  double overalloc = 2.;
314 
315  // "fact: iluk level-of-fill" parsing is more complicated, because
316  // we want to allow as many types as make sense. int is the native
317  // type, but we also want to accept double (for backwards
318  // compatibilty with ILUT). You can't cast arbitrary magnitude_type
319  // (e.g., Sacado::MP::Vector) to int, so we use float instead, to
320  // get coverage of the most common magnitude_type cases. Weirdly,
321  // there's an Ifpack2 test that sets the fill level as a
322  // global_ordinal_type.
323  {
324  const std::string paramName ("fact: iluk level-of-fill");
325  getParamTryingTypes<int, int, global_ordinal_type, double, float>
326  (fillLevel, params, paramName, prefix);
327  }
328  // For the other parameters, we prefer magnitude_type, but allow
329  // double for backwards compatibility.
330  {
331  const std::string paramName ("fact: absolute threshold");
332  getParamTryingTypes<magnitude_type, magnitude_type, double>
333  (absThresh, params, paramName, prefix);
334  }
335  {
336  const std::string paramName ("fact: relative threshold");
337  getParamTryingTypes<magnitude_type, magnitude_type, double>
338  (relThresh, params, paramName, prefix);
339  }
340  {
341  const std::string paramName ("fact: relax value");
342  getParamTryingTypes<magnitude_type, magnitude_type, double>
343  (relaxValue, params, paramName, prefix);
344  }
345  {
346  const std::string paramName ("fact: iluk overalloc");
347  getParamTryingTypes<double, double>
348  (overalloc, params, paramName, prefix);
349  }
350 
351  // Parsing implementation type
352  Details::IlukImplType::Enum ilukimplType = Details::IlukImplType::Serial;
353  do {
354  static const char typeName[] = "fact: type";
355 
356  if ( ! params.isType<std::string>(typeName)) break;
357 
358  // Map std::string <-> IlukImplType::Enum.
359  Array<std::string> ilukimplTypeStrs;
360  Array<Details::IlukImplType::Enum> ilukimplTypeEnums;
361  Details::IlukImplType::loadPLTypeOption (ilukimplTypeStrs, ilukimplTypeEnums);
362  Teuchos::StringToIntegralParameterEntryValidator<Details::IlukImplType::Enum>
363  s2i(ilukimplTypeStrs (), ilukimplTypeEnums (), typeName, false);
364 
365  ilukimplType = s2i.getIntegralValue(params.get<std::string>(typeName));
366  } while (0);
367 
368  if (ilukimplType == Details::IlukImplType::KSPILUK) {
369  this->isKokkosKernelsSpiluk_ = true;
370  }
371  else {
372  this->isKokkosKernelsSpiluk_ = false;
373  }
374 
375  // Forward to trisolvers.
376  L_solver_->setParameters(params);
377  U_solver_->setParameters(params);
378 
379  // "Commit" the values only after validating all of them. This
380  // ensures that there are no side effects if this routine throws an
381  // exception.
382 
383  LevelOfFill_ = fillLevel;
384  Overalloc_ = overalloc;
385 
386  // mfh 28 Nov 2012: The previous code would not assign Athresh_,
387  // Rthresh_, or RelaxValue_, if the read-in value was -1. I don't
388  // know if keeping this behavior is correct, but I'll keep it just
389  // so as not to change previous behavior.
390 
391  if (absThresh != -STM::one ()) {
392  Athresh_ = absThresh;
393  }
394  if (relThresh != -STM::one ()) {
395  Rthresh_ = relThresh;
396  }
397  if (relaxValue != -STM::one ()) {
398  RelaxValue_ = relaxValue;
399  }
400 }
401 
402 
403 template<class MatrixType>
404 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
406  return Teuchos::rcp_implicit_cast<const row_matrix_type> (A_);
407 }
408 
409 
410 template<class MatrixType>
411 Teuchos::RCP<const typename RILUK<MatrixType>::crs_matrix_type>
413  return Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_, true);
414 }
415 
416 
417 template<class MatrixType>
418 Teuchos::RCP<const typename RILUK<MatrixType>::row_matrix_type>
419 RILUK<MatrixType>::makeLocalFilter (const Teuchos::RCP<const row_matrix_type>& A)
420 {
421  using Teuchos::RCP;
422  using Teuchos::rcp;
423  using Teuchos::rcp_dynamic_cast;
424  using Teuchos::rcp_implicit_cast;
425 
426  // If A_'s communicator only has one process, or if its column and
427  // row Maps are the same, then it is already local, so use it
428  // directly.
429  if (A->getRowMap ()->getComm ()->getSize () == 1 ||
430  A->getRowMap ()->isSameAs (* (A->getColMap ()))) {
431  return A;
432  }
433 
434  // If A_ is already a LocalFilter, then use it directly. This
435  // should be the case if RILUK is being used through
436  // AdditiveSchwarz, for example.
437  RCP<const LocalFilter<row_matrix_type> > A_lf_r =
438  rcp_dynamic_cast<const LocalFilter<row_matrix_type> > (A);
439  if (! A_lf_r.is_null ()) {
440  return rcp_implicit_cast<const row_matrix_type> (A_lf_r);
441  }
442  else {
443  // A_'s communicator has more than one process, its row Map and
444  // its column Map differ, and A_ is not a LocalFilter. Thus, we
445  // have to wrap it in a LocalFilter.
446  return rcp (new LocalFilter<row_matrix_type> (A));
447  }
448 }
449 
450 
451 template<class MatrixType>
453 {
454  using Teuchos::RCP;
455  using Teuchos::rcp;
456  using Teuchos::rcp_const_cast;
457  using Teuchos::rcp_dynamic_cast;
458  using Teuchos::rcp_implicit_cast;
459  using Teuchos::Array;
460  using Teuchos::ArrayView;
461  typedef Tpetra::CrsGraph<local_ordinal_type,
463  node_type> crs_graph_type;
464  const char prefix[] = "Ifpack2::RILUK::initialize: ";
465 
466  TEUCHOS_TEST_FOR_EXCEPTION
467  (A_.is_null (), std::runtime_error, prefix << "The matrix is null. Please "
468  "call setMatrix() with a nonnull input before calling this method.");
469  TEUCHOS_TEST_FOR_EXCEPTION
470  (! A_->isFillComplete (), std::runtime_error, prefix << "The matrix is not "
471  "fill complete. You may not invoke initialize() or compute() with this "
472  "matrix until the matrix is fill complete. If your matrix is a "
473  "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
474  "range Maps, if appropriate) before calling this method.");
475 
476  Teuchos::Time timer ("RILUK::initialize");
477  double startTime = timer.wallTime();
478  { // Start timing
479  Teuchos::TimeMonitor timeMon (timer);
480 
481  // Calling initialize() means that the user asserts that the graph
482  // of the sparse matrix may have changed. We must not just reuse
483  // the previous graph in that case.
484  //
485  // Regarding setting isAllocated_ to false: Eventually, we may want
486  // some kind of clever memory reuse strategy, but it's always
487  // correct just to blow everything away and start over.
488  isInitialized_ = false;
489  isAllocated_ = false;
490  isComputed_ = false;
491  Graph_ = Teuchos::null;
492 
493  A_local_ = makeLocalFilter (A_);
494  TEUCHOS_TEST_FOR_EXCEPTION(
495  A_local_.is_null (), std::logic_error, "Ifpack2::RILUK::initialize: "
496  "makeLocalFilter returned null; it failed to compute A_local. "
497  "Please report this bug to the Ifpack2 developers.");
498 
499  // FIXME (mfh 24 Jan 2014, 26 Mar 2014) It would be more efficient
500  // to rewrite RILUK so that it works with any RowMatrix input, not
501  // just CrsMatrix. (That would require rewriting IlukGraph to
502  // handle a Tpetra::RowGraph.) However, to make it work for now,
503  // we just copy the input matrix if it's not a CrsMatrix.
504 
505  if (this->isKokkosKernelsSpiluk_) {
506  this->KernelHandle_ = Teuchos::rcp (new kk_handle_type ());
507  KernelHandle_->create_spiluk_handle( KokkosSparse::Experimental::SPILUKAlgorithm::SEQLVLSCHD_TP1,
508  A_local_->getNodeNumRows(),
509  2*A_local_->getNodeNumEntries()*(LevelOfFill_+1),
510  2*A_local_->getNodeNumEntries()*(LevelOfFill_+1) );
511  }
512 
513  {
514  RCP<const crs_matrix_type> A_local_crs =
515  rcp_dynamic_cast<const crs_matrix_type> (A_local_);
516  if (A_local_crs.is_null ()) {
517  local_ordinal_type numRows = A_local_->getNodeNumRows();
518  Array<size_t> entriesPerRow(numRows);
519  for(local_ordinal_type i = 0; i < numRows; i++) {
520  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
521  }
522  RCP<crs_matrix_type> A_local_crs_nc =
523  rcp (new crs_matrix_type (A_local_->getRowMap (),
524  A_local_->getColMap (),
525  entriesPerRow()));
526  // copy entries into A_local_crs
527  nonconst_local_inds_host_view_type indices("indices",A_local_->getNodeMaxNumRowEntries());
528  nonconst_values_host_view_type values("values",A_local_->getNodeMaxNumRowEntries());
529  for(local_ordinal_type i = 0; i < numRows; i++) {
530  size_t numEntries = 0;
531  A_local_->getLocalRowCopy(i, indices, values, numEntries);
532  A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()), indices.data());
533  }
534  A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
535  A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
536  }
537  Graph_ = rcp (new Ifpack2::IlukGraph<crs_graph_type,kk_handle_type> (A_local_crs->getCrsGraph (),
538  LevelOfFill_, 0, Overalloc_));
539  }
540 
541  if (this->isKokkosKernelsSpiluk_) Graph_->initialize (KernelHandle_);
542  else Graph_->initialize ();
543 
544  allocate_L_and_U ();
545  checkOrderingConsistency (*A_local_);
546  L_solver_->setMatrix (L_);
547  L_solver_->initialize ();
548  U_solver_->setMatrix (U_);
549  U_solver_->initialize ();
550 
551  // Do not call initAllValues. compute() always calls initAllValues to
552  // fill L and U with possibly new numbers. initialize() is concerned
553  // only with the nonzero pattern.
554  } // Stop timing
555 
556  isInitialized_ = true;
557  ++numInitialize_;
558  initializeTime_ += (timer.wallTime() - startTime);
559 }
560 
561 template<class MatrixType>
562 void
564 checkOrderingConsistency (const row_matrix_type& A)
565 {
566  // First check that the local row map ordering is the same as the local portion of the column map.
567  // The extraction of the strictly lower/upper parts of A, as well as the factorization,
568  // implicitly assume that this is the case.
569  Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getNodeElementList();
570  Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getNodeElementList();
571  bool gidsAreConsistentlyOrdered=true;
572  global_ordinal_type indexOfInconsistentGID=0;
573  for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
574  if (rowGIDs[i] != colGIDs[i]) {
575  gidsAreConsistentlyOrdered=false;
576  indexOfInconsistentGID=i;
577  break;
578  }
579  }
580  TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==false, std::runtime_error,
581  "The ordering of the local GIDs in the row and column maps is not the same"
582  << std::endl << "at index " << indexOfInconsistentGID
583  << ". Consistency is required, as all calculations are done with"
584  << std::endl << "local indexing.");
585 }
586 
587 template<class MatrixType>
588 void
589 RILUK<MatrixType>::
590 initAllValues (const row_matrix_type& A)
591 {
592  using Teuchos::ArrayRCP;
593  using Teuchos::Comm;
594  using Teuchos::ptr;
595  using Teuchos::RCP;
596  using Teuchos::REDUCE_SUM;
597  using Teuchos::reduceAll;
598  typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
599 
600  size_t NumIn = 0, NumL = 0, NumU = 0;
601  bool DiagFound = false;
602  size_t NumNonzeroDiags = 0;
603  size_t MaxNumEntries = A.getGlobalMaxNumRowEntries();
604 
605  // Allocate temporary space for extracting the strictly
606  // lower and upper parts of the matrix A.
607  nonconst_local_inds_host_view_type InI("InI",MaxNumEntries);
608  Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
609  Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
610  nonconst_values_host_view_type InV("InV",MaxNumEntries);
611  Teuchos::Array<scalar_type> LV(MaxNumEntries);
612  Teuchos::Array<scalar_type> UV(MaxNumEntries);
613 
614  // Check if values should be inserted or replaced
615  const bool ReplaceValues = L_->isStaticGraph () || L_->isLocallyIndexed ();
616 
617  L_->resumeFill ();
618  U_->resumeFill ();
619  if (ReplaceValues) {
620  L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
621  U_->setAllToScalar (STS::zero ());
622  }
623 
624  D_->putScalar (STS::zero ()); // Set diagonal values to zero
625  auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
626 
627  RCP<const map_type> rowMap = L_->getRowMap ();
628 
629  // First we copy the user's matrix into L and U, regardless of fill level.
630  // It is important to note that L and U are populated using local indices.
631  // This means that if the row map GIDs are not monotonically increasing
632  // (i.e., permuted or gappy), then the strictly lower (upper) part of the
633  // matrix is not the one that you would get if you based L (U) on GIDs.
634  // This is ok, as the *order* of the GIDs in the rowmap is a better
635  // expression of the user's intent than the GIDs themselves.
636 
637  Teuchos::ArrayView<const global_ordinal_type> nodeGIDs = rowMap->getNodeElementList();
638  for (size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
639  local_ordinal_type local_row = myRow;
640 
641  //TODO JJH 4April2014 An optimization is to use getLocalRowView. Not all matrices support this,
642  // we'd need to check via the Tpetra::RowMatrix method supportsRowViews().
643  A.getLocalRowCopy (local_row, InI, InV, NumIn); // Get Values and Indices
644 
645  // Split into L and U (we don't assume that indices are ordered).
646 
647  NumL = 0;
648  NumU = 0;
649  DiagFound = false;
650 
651  for (size_t j = 0; j < NumIn; ++j) {
652  const local_ordinal_type k = InI[j];
653 
654  if (k == local_row) {
655  DiagFound = true;
656  // Store perturbed diagonal in Tpetra::Vector D_
657  DV(local_row) += Rthresh_ * InV[j] + IFPACK2_SGN(InV[j]) * Athresh_;
658  }
659  else if (k < 0) { // Out of range
660  TEUCHOS_TEST_FOR_EXCEPTION(
661  true, std::runtime_error, "Ifpack2::RILUK::initAllValues: current "
662  "GID k = " << k << " < 0. I'm not sure why this is an error; it is "
663  "probably an artifact of the undocumented assumptions of the "
664  "original implementation (likely copied and pasted from Ifpack). "
665  "Nevertheless, the code I found here insisted on this being an error "
666  "state, so I will throw an exception here.");
667  }
668  else if (k < local_row) {
669  LI[NumL] = k;
670  LV[NumL] = InV[j];
671  NumL++;
672  }
673  else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
674  UI[NumU] = k;
675  UV[NumU] = InV[j];
676  NumU++;
677  }
678  }
679 
680  // Check in things for this row of L and U
681 
682  if (DiagFound) {
683  ++NumNonzeroDiags;
684  } else {
685  DV(local_row) = Athresh_;
686  }
687 
688  if (NumL) {
689  if (ReplaceValues) {
690  L_->replaceLocalValues(local_row, LI(0, NumL), LV(0,NumL));
691  } else {
692  //FIXME JJH 24April2014 Is this correct? I believe this case is when there aren't already values
693  //FIXME in this row in the column locations corresponding to UI.
694  L_->insertLocalValues(local_row, LI(0,NumL), LV(0,NumL));
695  }
696  }
697 
698  if (NumU) {
699  if (ReplaceValues) {
700  U_->replaceLocalValues(local_row, UI(0,NumU), UV(0,NumU));
701  } else {
702  //FIXME JJH 24April2014 Is this correct? I believe this case is when there aren't already values
703  //FIXME in this row in the column locations corresponding to UI.
704  U_->insertLocalValues(local_row, UI(0,NumU), UV(0,NumU));
705  }
706  }
707  }
708 
709  // At this point L and U have the values of A in the structure of L
710  // and U, and diagonal vector D
711 
712  isInitialized_ = true;
713 }
714 
715 
716 template<class MatrixType>
718 {
719  using Teuchos::RCP;
720  using Teuchos::rcp;
721  using Teuchos::rcp_const_cast;
722  using Teuchos::rcp_dynamic_cast;
723  using Teuchos::Array;
724  using Teuchos::ArrayView;
725  const char prefix[] = "Ifpack2::RILUK::compute: ";
726 
727  // initialize() checks this too, but it's easier for users if the
728  // error shows them the name of the method that they actually
729  // called, rather than the name of some internally called method.
730  TEUCHOS_TEST_FOR_EXCEPTION
731  (A_.is_null (), std::runtime_error, prefix << "The matrix is null. Please "
732  "call setMatrix() with a nonnull input before calling this method.");
733  TEUCHOS_TEST_FOR_EXCEPTION
734  (! A_->isFillComplete (), std::runtime_error, prefix << "The matrix is not "
735  "fill complete. You may not invoke initialize() or compute() with this "
736  "matrix until the matrix is fill complete. If your matrix is a "
737  "Tpetra::CrsMatrix, please call fillComplete on it (with the domain and "
738  "range Maps, if appropriate) before calling this method.");
739 
740  if (! isInitialized ()) {
741  initialize (); // Don't count this in the compute() time
742  }
743 
744  Teuchos::Time timer ("RILUK::compute");
745 
746  // Start timing
747  Teuchos::TimeMonitor timeMon (timer);
748  double startTime = timer.wallTime();
749 
750  isComputed_ = false;
751 
752  if (!this->isKokkosKernelsSpiluk_) {
753  // Fill L and U with numbers. This supports nonzero pattern reuse by calling
754  // initialize() once and then compute() multiple times.
755  initAllValues (*A_local_);
756 
757  // MinMachNum should be officially defined, for now pick something a little
758  // bigger than IEEE underflow value
759 
760  const scalar_type MinDiagonalValue = STS::rmin ();
761  const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
762 
763  size_t NumIn, NumL, NumU;
764 
765  // Get Maximum Row length
766  const size_t MaxNumEntries =
767  L_->getNodeMaxNumRowEntries () + U_->getNodeMaxNumRowEntries () + 1;
768 
769  Teuchos::Array<local_ordinal_type> InI(MaxNumEntries); // Allocate temp space
770  Teuchos::Array<scalar_type> InV(MaxNumEntries);
771  size_t num_cols = U_->getColMap()->getNodeNumElements();
772  Teuchos::Array<int> colflag(num_cols);
773 
774  auto DV = Kokkos::subview(D_->getLocalViewHost(Tpetra::Access::ReadWrite), Kokkos::ALL(), 0);
775 
776  // Now start the factorization.
777 
778  for (size_t j = 0; j < num_cols; ++j) {
779  colflag[j] = -1;
780  }
781  using IST = typename row_matrix_type::impl_scalar_type;
782  for (size_t i = 0; i < L_->getNodeNumRows (); ++i) {
783  local_ordinal_type local_row = i;
784  // Need some integer workspace and pointers
785  size_t NumUU;
786  local_inds_host_view_type UUI;
787  values_host_view_type UUV;
788 
789  // Fill InV, InI with current row of L, D and U combined
790 
791  NumIn = MaxNumEntries;
792  nonconst_local_inds_host_view_type InI_v(InI.data(),MaxNumEntries);
793  nonconst_values_host_view_type InV_v(reinterpret_cast<IST*>(InV.data()),MaxNumEntries);
794 
795  L_->getLocalRowCopy (local_row, InI_v , InV_v, NumL);
796 
797  InV[NumL] = DV(i); // Put in diagonal
798  InI[NumL] = local_row;
799 
800  nonconst_local_inds_host_view_type InI_sub(InI.data()+NumL+1,MaxNumEntries-NumL-1);
801  nonconst_values_host_view_type InV_sub(reinterpret_cast<IST*>(InV.data())+NumL+1,MaxNumEntries-NumL-1);
802 
803  U_->getLocalRowCopy (local_row, InI_sub,InV_sub, NumU);
804  NumIn = NumL+NumU+1;
805 
806  // Set column flags
807  for (size_t j = 0; j < NumIn; ++j) {
808  colflag[InI[j]] = j;
809  }
810 
811  scalar_type diagmod = STS::zero (); // Off-diagonal accumulator
812 
813  for (size_t jj = 0; jj < NumL; ++jj) {
814  local_ordinal_type j = InI[jj];
815  IST multiplier = InV[jj]; // current_mults++;
816 
817  InV[jj] *= static_cast<scalar_type>(DV(j));
818 
819  U_->getLocalRowView(j, UUI, UUV); // View of row above
820  NumUU = UUI.size();
821 
822  if (RelaxValue_ == STM::zero ()) {
823  for (size_t k = 0; k < NumUU; ++k) {
824  const int kk = colflag[UUI[k]];
825  // FIXME (mfh 23 Dec 2013) Wait a second, we just set
826  // colflag above using size_t (which is generally unsigned),
827  // but now we're querying it using int (which is signed).
828  if (kk > -1) {
829  InV[kk] -= static_cast<scalar_type>(multiplier * UUV[k]);
830  }
831  }
832 
833  }
834  else {
835  for (size_t k = 0; k < NumUU; ++k) {
836  // FIXME (mfh 23 Dec 2013) Wait a second, we just set
837  // colflag above using size_t (which is generally unsigned),
838  // but now we're querying it using int (which is signed).
839  const int kk = colflag[UUI[k]];
840  if (kk > -1) {
841  InV[kk] -= static_cast<scalar_type>(multiplier*UUV[k]);
842  }
843  else {
844  diagmod -= static_cast<scalar_type>(multiplier*UUV[k]);
845  }
846  }
847  }
848  }
849 
850  if (NumL) {
851  // Replace current row of L
852  L_->replaceLocalValues (local_row, InI (0, NumL), InV (0, NumL));
853  }
854 
855  DV(i) = InV[NumL]; // Extract Diagonal value
856 
857  if (RelaxValue_ != STM::zero ()) {
858  DV(i) += RelaxValue_*diagmod; // Add off diagonal modifications
859  }
860 
861  if (STS::magnitude (DV(i)) > STS::magnitude (MaxDiagonalValue)) {
862  if (STS::real (DV(i)) < STM::zero ()) {
863  DV(i) = -MinDiagonalValue;
864  }
865  else {
866  DV(i) = MinDiagonalValue;
867  }
868  }
869  else {
870  DV(i) = static_cast<impl_scalar_type>(STS::one ()) / DV(i); // Invert diagonal value
871  }
872 
873  for (size_t j = 0; j < NumU; ++j) {
874  InV[NumL+1+j] *= static_cast<scalar_type>(DV(i)); // Scale U by inverse of diagonal
875  }
876 
877  if (NumU) {
878  // Replace current row of L and U
879  U_->replaceLocalValues (local_row, InI (NumL+1, NumU), InV (NumL+1, NumU));
880  }
881 
882  // Reset column flags
883  for (size_t j = 0; j < NumIn; ++j) {
884  colflag[InI[j]] = -1;
885  }
886  }
887 
888  // The domain of L and the range of U are exactly their own row maps
889  // (there is no communication). The domain of U and the range of L
890  // must be the same as those of the original matrix, However if the
891  // original matrix is a VbrMatrix, these two latter maps are
892  // translation from a block map to a point map.
893  // FIXME (mfh 23 Dec 2013) Do we know that the column Map of L_ is
894  // always one-to-one?
895  L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
896  U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
897 
898  // If L_solver_ or U_solver store modified factors internally, we need to reset those
899  L_solver_->setMatrix (L_);
900  L_solver_->compute ();
901  U_solver_->setMatrix (U_);
902  U_solver_->compute ();
903  }
904  else {
905  {//Make sure values in A is picked up even in case of pattern reuse
906  RCP<const crs_matrix_type> A_local_crs =
907  rcp_dynamic_cast<const crs_matrix_type> (A_local_);
908  if (A_local_crs.is_null ()) {
909  local_ordinal_type numRows = A_local_->getNodeNumRows();
910  Array<size_t> entriesPerRow(numRows);
911  for(local_ordinal_type i = 0; i < numRows; i++) {
912  entriesPerRow[i] = A_local_->getNumEntriesInLocalRow(i);
913  }
914  RCP<crs_matrix_type> A_local_crs_nc =
915  rcp (new crs_matrix_type (A_local_->getRowMap (),
916  A_local_->getColMap (),
917  entriesPerRow()));
918  // copy entries into A_local_crs
919  nonconst_local_inds_host_view_type indices("indices",A_local_->getNodeMaxNumRowEntries());
920  nonconst_values_host_view_type values("values",A_local_->getNodeMaxNumRowEntries());
921  for(local_ordinal_type i = 0; i < numRows; i++) {
922  size_t numEntries = 0;
923  A_local_->getLocalRowCopy(i, indices, values, numEntries);
924  A_local_crs_nc->insertLocalValues(i, numEntries, reinterpret_cast<scalar_type*>(values.data()),indices.data());
925  }
926  A_local_crs_nc->fillComplete (A_local_->getDomainMap (), A_local_->getRangeMap ());
927  A_local_crs = rcp_const_cast<const crs_matrix_type> (A_local_crs_nc);
928  }
929  A_local_rowmap_ = A_local_crs->getLocalMatrixDevice().graph.row_map;
930  A_local_entries_ = A_local_crs->getLocalMatrixDevice().graph.entries;
931  A_local_values_ = A_local_crs->getLocalValuesView();
932  }
933 
934  L_->resumeFill ();
935  U_->resumeFill ();
936 
937  if (L_->isStaticGraph () || L_->isLocallyIndexed ()) {
938  L_->setAllToScalar (STS::zero ()); // Zero out L and U matrices
939  U_->setAllToScalar (STS::zero ());
940  }
941 
942  using row_map_type = typename crs_matrix_type::local_matrix_device_type::row_map_type;
943 
944  row_map_type L_rowmap = L_->getLocalMatrixDevice().graph.row_map;
945  auto L_entries = L_->getLocalMatrixDevice().graph.entries;
946  auto L_values = L_->getLocalValuesView();
947  row_map_type U_rowmap = U_->getLocalMatrixDevice().graph.row_map;
948  auto U_entries = U_->getLocalMatrixDevice().graph.entries;
949  auto U_values = U_->getLocalValuesView();
950 
951  KokkosSparse::Experimental::spiluk_numeric( KernelHandle_.getRawPtr(), LevelOfFill_,
952  A_local_rowmap_, A_local_entries_, A_local_values_,
953  L_rowmap, L_entries, L_values, U_rowmap, U_entries, U_values );
954 
955  L_->fillComplete (L_->getColMap (), A_local_->getRangeMap ());
956  U_->fillComplete (A_local_->getDomainMap (), U_->getRowMap ());
957 
958  L_solver_->setMatrix (L_);
959  L_solver_->compute ();
960  U_solver_->setMatrix (U_);
961  U_solver_->compute ();
962  }
963 
964  isComputed_ = true;
965  ++numCompute_;
966  computeTime_ += (timer.wallTime() - startTime);
967 }
968 
969 
970 template<class MatrixType>
971 void
973 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
974  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
975  Teuchos::ETransp mode,
976  scalar_type alpha,
977  scalar_type beta) const
978 {
979  using Teuchos::RCP;
980  using Teuchos::rcpFromRef;
981 
982  TEUCHOS_TEST_FOR_EXCEPTION(
983  A_.is_null (), std::runtime_error, "Ifpack2::RILUK::apply: The matrix is "
984  "null. Please call setMatrix() with a nonnull input, then initialize() "
985  "and compute(), before calling this method.");
986  TEUCHOS_TEST_FOR_EXCEPTION(
987  ! isComputed (), std::runtime_error,
988  "Ifpack2::RILUK::apply: If you have not yet called compute(), "
989  "you must call compute() before calling this method.");
990  TEUCHOS_TEST_FOR_EXCEPTION(
991  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
992  "Ifpack2::RILUK::apply: X and Y do not have the same number of columns. "
993  "X.getNumVectors() = " << X.getNumVectors ()
994  << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
995  TEUCHOS_TEST_FOR_EXCEPTION(
996  STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
997  "Ifpack2::RILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
998  "complex Scalar type. Please talk to the Ifpack2 developers to get this "
999  "fixed. There is a FIXME in this file about this very issue.");
1000 #ifdef HAVE_IFPACK2_DEBUG
1001  {
1002  const magnitude_type D_nrm1 = D_->norm1 ();
1003  TEUCHOS_TEST_FOR_EXCEPTION( STM::isnaninf (D_nrm1), std::runtime_error, "Ifpack2::RILUK::apply: The 1-norm of the stored diagonal is NaN or Inf.");
1004  Teuchos::Array<magnitude_type> norms (X.getNumVectors ());
1005  X.norm1 (norms ());
1006  bool good = true;
1007  for (size_t j = 0; j < X.getNumVectors (); ++j) {
1008  if (STM::isnaninf (norms[j])) {
1009  good = false;
1010  break;
1011  }
1012  }
1013  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::RILUK::apply: The 1-norm of the input X is NaN or Inf.");
1014  }
1015 #endif // HAVE_IFPACK2_DEBUG
1016 
1017  const scalar_type one = STS::one ();
1018  const scalar_type zero = STS::zero ();
1019 
1020  Teuchos::Time timer ("RILUK::apply");
1021  double startTime = timer.wallTime();
1022  { // Start timing
1023  Teuchos::TimeMonitor timeMon (timer);
1024  if (alpha == one && beta == zero) {
1025  if (mode == Teuchos::NO_TRANS) { // Solve L (D (U Y)) = X for Y.
1026  // Start by solving L Y = X for Y.
1027  L_solver_->apply (X, Y, mode);
1028 
1029  if (!this->isKokkosKernelsSpiluk_) {
1030  // Solve D Y = Y. The operation lets us do this in place in Y, so we can
1031  // write "solve D Y = Y for Y."
1032  Y.elementWiseMultiply (one, *D_, Y, zero);
1033  }
1034 
1035  U_solver_->apply (Y, Y, mode); // Solve U Y = Y.
1036  }
1037  else { // Solve U^P (D^P (L^P Y)) = X for Y (where P is * or T).
1038  // Start by solving U^P Y = X for Y.
1039  U_solver_->apply (X, Y, mode);
1040 
1041  if (!this->isKokkosKernelsSpiluk_) {
1042  // Solve D^P Y = Y.
1043  //
1044  // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we
1045  // need to do an elementwise multiply with the conjugate of
1046  // D_, not just with D_ itself.
1047  Y.elementWiseMultiply (one, *D_, Y, zero);
1048  }
1049 
1050  L_solver_->apply (Y, Y, mode); // Solve L^P Y = Y.
1051  }
1052  }
1053  else { // alpha != 1 or beta != 0
1054  if (alpha == zero) {
1055  // The special case for beta == 0 ensures that if Y contains Inf
1056  // or NaN values, we replace them with 0 (following BLAS
1057  // convention), rather than multiplying them by 0 to get NaN.
1058  if (beta == zero) {
1059  Y.putScalar (zero);
1060  } else {
1061  Y.scale (beta);
1062  }
1063  } else { // alpha != zero
1064  MV Y_tmp (Y.getMap (), Y.getNumVectors ());
1065  apply (X, Y_tmp, mode);
1066  Y.update (alpha, Y_tmp, beta);
1067  }
1068  }
1069  }//end timing
1070 
1071 #ifdef HAVE_IFPACK2_DEBUG
1072  {
1073  Teuchos::Array<magnitude_type> norms (Y.getNumVectors ());
1074  Y.norm1 (norms ());
1075  bool good = true;
1076  for (size_t j = 0; j < Y.getNumVectors (); ++j) {
1077  if (STM::isnaninf (norms[j])) {
1078  good = false;
1079  break;
1080  }
1081  }
1082  TEUCHOS_TEST_FOR_EXCEPTION( ! good, std::runtime_error, "Ifpack2::RILUK::apply: The 1-norm of the output Y is NaN or Inf.");
1083  }
1084 #endif // HAVE_IFPACK2_DEBUG
1085 
1086  ++numApply_;
1087  applyTime_ += (timer.wallTime() - startTime);
1088 }
1089 
1090 
1091 //VINH comment out since multiply() is not needed anywhere
1092 //template<class MatrixType>
1093 //void RILUK<MatrixType>::
1094 //multiply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
1095 // Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
1096 // const Teuchos::ETransp mode) const
1097 //{
1098 // const scalar_type zero = STS::zero ();
1099 // const scalar_type one = STS::one ();
1100 //
1101 // if (mode != Teuchos::NO_TRANS) {
1102 // U_->apply (X, Y, mode); //
1103 // Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
1104 //
1105 // // FIXME (mfh 24 Jan 2014) If mode = Teuchos::CONJ_TRANS, we need
1106 // // to do an elementwise multiply with the conjugate of D_, not
1107 // // just with D_ itself.
1108 // Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
1109 //
1110 // MV Y_tmp (Y, Teuchos::Copy); // Need a temp copy of Y
1111 // L_->apply (Y_tmp, Y, mode);
1112 // Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
1113 // }
1114 // else {
1115 // L_->apply (X, Y, mode);
1116 // Y.update (one, X, one); // Y = Y + X (account for implicit unit diagonal)
1117 // Y.elementWiseMultiply (one, *D_, Y, zero); // y = D*y (D_ has inverse of diagonal)
1118 // MV Y_tmp (Y, Teuchos::Copy); // Need a temp copy of Y1
1119 // U_->apply (Y_tmp, Y, mode);
1120 // Y.update (one, Y_tmp, one); // (account for implicit unit diagonal)
1121 // }
1122 //}
1123 
1124 template<class MatrixType>
1126 {
1127  std::ostringstream os;
1128 
1129  // Output is a valid YAML dictionary in flow style. If you don't
1130  // like everything on a single line, you should call describe()
1131  // instead.
1132  os << "\"Ifpack2::RILUK\": {";
1133  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
1134  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
1135 
1136  os << "Level-of-fill: " << getLevelOfFill() << ", ";
1137 
1138  if (A_.is_null ()) {
1139  os << "Matrix: null";
1140  }
1141  else {
1142  os << "Global matrix dimensions: ["
1143  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]"
1144  << ", Global nnz: " << A_->getGlobalNumEntries();
1145  }
1146 
1147  if (! L_solver_.is_null ()) os << ", " << L_solver_->description ();
1148  if (! U_solver_.is_null ()) os << ", " << U_solver_->description ();
1149 
1150  os << "}";
1151  return os.str ();
1152 }
1153 
1154 } // namespace Ifpack2
1155 
1156 #define IFPACK2_RILUK_INSTANT(S,LO,GO,N) \
1157  template class Ifpack2::RILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
1158 
1159 #endif
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:263
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_RILUK_decl.hpp:281
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_RILUK_def.hpp:137
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:245
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Tpetra::CrsMatrix specialization used by this class for representing L and U.
Definition: Ifpack2_RILUK_decl.hpp:284
Ifpack2 implementation details.
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:266
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:260
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Returns the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_RILUK_def.hpp:234
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:79
crs_matrix_type::impl_scalar_type impl_scalar_type
Scalar type stored in Kokkos::Views (CrsMatrix and MultiVector)
Definition: Ifpack2_RILUK_decl.hpp:293
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:100
Tpetra::global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_RILUK_decl.hpp:546
Definition: Ifpack2_Container_decl.hpp:576
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Returns the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_RILUK_def.hpp:253
Tpetra::RowMatrix< MatrixType::scalar_type, MatrixType::local_ordinal_type, MatrixType::global_ordinal_type, MatrixType::node_type > ::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:257