Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_ILUT_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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 // disable clang warnings
47 #if defined (__clang__) && !defined (__INTEL_COMPILER)
48 #pragma clang system_header
49 #endif
50 
51 #include "Ifpack2_Heap.hpp"
52 #include "Ifpack2_LocalFilter.hpp"
53 #include "Ifpack2_LocalSparseTriangularSolver_decl.hpp"
54 #include "Ifpack2_Parameters.hpp"
55 #include "Tpetra_CrsMatrix.hpp"
56 #include "Teuchos_Time.hpp"
57 #include "Teuchos_TypeNameTraits.hpp"
58 
59 namespace Ifpack2 {
60 
61  namespace {
62 
87  template<class ScalarType>
88  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
89  ilutDefaultDropTolerance () {
90  using Teuchos::as;
91  typedef Teuchos::ScalarTraits<ScalarType> STS;
92  typedef typename STS::magnitudeType magnitude_type;
93  typedef Teuchos::ScalarTraits<magnitude_type> STM;
94 
95  // 1/2. Hopefully this can be represented in magnitude_type.
96  const magnitude_type oneHalf = STM::one() / (STM::one() + STM::one());
97 
98  // The min ensures that in case magnitude_type has very low
99  // precision, we'll at least get some value strictly less than
100  // one.
101  return std::min (as<magnitude_type> (1000) * STS::magnitude (STS::eps ()), oneHalf);
102  }
103 
104  // Full specialization for ScalarType = double.
105  // This specialization preserves ILUT's previous default behavior.
106  template<>
107  Teuchos::ScalarTraits<double>::magnitudeType
108  ilutDefaultDropTolerance<double> () {
109  return 1e-12;
110  }
111 
112  } // namespace (anonymous)
113 
114 
115 template <class MatrixType>
116 ILUT<MatrixType>::ILUT (const Teuchos::RCP<const row_matrix_type>& A) :
117  A_ (A),
118  Athresh_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
119  Rthresh_ (Teuchos::ScalarTraits<magnitude_type>::one ()),
120  RelaxValue_ (Teuchos::ScalarTraits<magnitude_type>::zero ()),
121  LevelOfFill_ (1),
122  DropTolerance_ (ilutDefaultDropTolerance<scalar_type> ()),
123  InitializeTime_ (0.0),
124  ComputeTime_ (0.0),
125  ApplyTime_ (0.0),
126  NumInitialize_ (0),
127  NumCompute_ (0),
128  NumApply_ (0),
129  IsInitialized_ (false),
130  IsComputed_ (false)
131 {}
132 
133 template <class MatrixType>
135 {}
136 
137 template <class MatrixType>
138 void ILUT<MatrixType>::setParameters (const Teuchos::ParameterList& params)
139 {
140  using Teuchos::as;
141  using Teuchos::Exceptions::InvalidParameterName;
142  using Teuchos::Exceptions::InvalidParameterType;
143 
144  // Default values of the various parameters.
145  int fillLevel = 1;
146  magnitude_type absThresh = STM::zero ();
147  magnitude_type relThresh = STM::one ();
148  magnitude_type relaxValue = STM::zero ();
149  magnitude_type dropTol = ilutDefaultDropTolerance<scalar_type> ();
150 
151  bool gotFillLevel = false;
152  try {
153  // Try getting the fill level as an int.
154  fillLevel = params.get<int> ("fact: ilut level-of-fill");
155  gotFillLevel = true;
156  }
157  catch (InvalidParameterName&) {
158  // We didn't really get it, but this tells us to stop looking.
159  gotFillLevel = true;
160  }
161  catch (InvalidParameterType&) {
162  // The name is right, but the type is wrong; try different types.
163  // We don't have to check InvalidParameterName again, since we
164  // checked it above, and it has nothing to do with the type.
165  }
166 
167  if (! gotFillLevel) {
168  // Try magnitude_type, for backwards compatibility.
169  try {
170  fillLevel = as<int> (params.get<magnitude_type> ("fact: ilut level-of-fill"));
171  }
172  catch (InvalidParameterType&) {}
173  }
174  if (! gotFillLevel) {
175  // Try double, for backwards compatibility.
176  try {
177  fillLevel = as<int> (params.get<double> ("fact: ilut level-of-fill"));
178  }
179  catch (InvalidParameterType&) {}
180  }
181  // If none of the above attempts succeed, accept the default value.
182 
183  TEUCHOS_TEST_FOR_EXCEPTION(
184  fillLevel <= 0, std::runtime_error,
185  "Ifpack2::ILUT: The \"fact: ilut level-of-fill\" parameter must be "
186  "strictly greater than zero, but you specified a value of " << fillLevel
187  << ". Remember that for ILUT, the fill level p means something different "
188  "than it does for ILU(k). ILU(0) produces factors with the same sparsity "
189  "structure as the input matrix A; ILUT with p = 0 always produces a "
190  "diagonal matrix, and is thus probably not what you want.");
191 
192  try {
193  absThresh = params.get<magnitude_type> ("fact: absolute threshold");
194  }
195  catch (InvalidParameterType&) {
196  // Try double, for backwards compatibility.
197  // The cast from double to magnitude_type must succeed.
198  absThresh = as<magnitude_type> (params.get<double> ("fact: absolute threshold"));
199  }
200  catch (InvalidParameterName&) {
201  // Accept the default value.
202  }
203 
204  try {
205  relThresh = params.get<magnitude_type> ("fact: relative threshold");
206  }
207  catch (InvalidParameterType&) {
208  // Try double, for backwards compatibility.
209  // The cast from double to magnitude_type must succeed.
210  relThresh = as<magnitude_type> (params.get<double> ("fact: relative threshold"));
211  }
212  catch (InvalidParameterName&) {
213  // Accept the default value.
214  }
215 
216  try {
217  relaxValue = params.get<magnitude_type> ("fact: relax value");
218  }
219  catch (InvalidParameterType&) {
220  // Try double, for backwards compatibility.
221  // The cast from double to magnitude_type must succeed.
222  relaxValue = as<magnitude_type> (params.get<double> ("fact: relax value"));
223  }
224  catch (InvalidParameterName&) {
225  // Accept the default value.
226  }
227 
228  try {
229  dropTol = params.get<magnitude_type> ("fact: drop tolerance");
230  }
231  catch (InvalidParameterType&) {
232  // Try double, for backwards compatibility.
233  // The cast from double to magnitude_type must succeed.
234  dropTol = as<magnitude_type> (params.get<double> ("fact: drop tolerance"));
235  }
236  catch (InvalidParameterName&) {
237  // Accept the default value.
238  }
239 
240  // "Commit" the values only after validating all of them. This
241  // ensures that there are no side effects if this routine throws an
242  // exception.
243 
244  // mfh 28 Nov 2012: The previous code would not assign Athresh_,
245  // Rthresh_, RelaxValue_, or DropTolerance_ if the read-in value was
246  // -1. I don't know if keeping this behavior is correct, but I'll
247  // keep it just so as not to change previous behavior.
248 
249  LevelOfFill_ = fillLevel;
250  if (absThresh != -STM::one ()) {
251  Athresh_ = absThresh;
252  }
253  if (relThresh != -STM::one ()) {
254  Rthresh_ = relThresh;
255  }
256  if (relaxValue != -STM::one ()) {
257  RelaxValue_ = relaxValue;
258  }
259  if (dropTol != -STM::one ()) {
260  DropTolerance_ = dropTol;
261  }
262 }
263 
264 
265 template <class MatrixType>
266 Teuchos::RCP<const Teuchos::Comm<int> >
268  TEUCHOS_TEST_FOR_EXCEPTION(
269  A_.is_null (), std::runtime_error, "Ifpack2::ILUT::getComm: "
270  "The matrix is null. Please call setMatrix() with a nonnull input "
271  "before calling this method.");
272  return A_->getComm ();
273 }
274 
275 
276 template <class MatrixType>
277 Teuchos::RCP<const typename ILUT<MatrixType>::row_matrix_type>
279  return A_;
280 }
281 
282 
283 template <class MatrixType>
284 Teuchos::RCP<const typename ILUT<MatrixType>::map_type>
286 {
287  TEUCHOS_TEST_FOR_EXCEPTION(
288  A_.is_null (), std::runtime_error, "Ifpack2::ILUT::getDomainMap: "
289  "The matrix is null. Please call setMatrix() with a nonnull input "
290  "before calling this method.");
291  return A_->getDomainMap ();
292 }
293 
294 
295 template <class MatrixType>
296 Teuchos::RCP<const typename ILUT<MatrixType>::map_type>
298 {
299  TEUCHOS_TEST_FOR_EXCEPTION(
300  A_.is_null (), std::runtime_error, "Ifpack2::ILUT::getRangeMap: "
301  "The matrix is null. Please call setMatrix() with a nonnull input "
302  "before calling this method.");
303  return A_->getRangeMap ();
304 }
305 
306 
307 template <class MatrixType>
309  return true;
310 }
311 
312 
313 template <class MatrixType>
315  return NumInitialize_;
316 }
317 
318 
319 template <class MatrixType>
321  return NumCompute_;
322 }
323 
324 
325 template <class MatrixType>
327  return NumApply_;
328 }
329 
330 
331 template <class MatrixType>
333  return InitializeTime_;
334 }
335 
336 
337 template<class MatrixType>
339  return ComputeTime_;
340 }
341 
342 
343 template<class MatrixType>
345  return ApplyTime_;
346 }
347 
348 
349 template<class MatrixType>
351  return L_->getGlobalNumEntries () + U_->getGlobalNumEntries ();
352 }
353 
354 
355 template<class MatrixType>
357  return L_->getNodeNumEntries () + U_->getNodeNumEntries ();
358 }
359 
360 
361 template<class MatrixType>
362 void ILUT<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
363 {
364  if (A.getRawPtr () != A_.getRawPtr ()) {
365  // Check in serial or one-process mode if the matrix is square.
366  TEUCHOS_TEST_FOR_EXCEPTION(
367  ! A.is_null () && A->getComm ()->getSize () == 1 &&
368  A->getNodeNumRows () != A->getNodeNumCols (),
369  std::runtime_error, "Ifpack2::ILUT::setMatrix: If A's communicator only "
370  "contains one process, then A must be square. Instead, you provided a "
371  "matrix A with " << A->getNodeNumRows () << " rows and "
372  << A->getNodeNumCols () << " columns.");
373 
374  // It's legal for A to be null; in that case, you may not call
375  // initialize() until calling setMatrix() with a nonnull input.
376  // Regardless, setting the matrix invalidates any previous
377  // factorization.
378  IsInitialized_ = false;
379  IsComputed_ = false;
380  A_local_ = Teuchos::null;
381 
382  // The sparse triangular solvers get a triangular factor as their
383  // input matrix. The triangular factors L_ and U_ are getting
384  // reset, so we reset the solvers' matrices to null. Do that
385  // before setting L_ and U_ to null, so that latter step actually
386  // frees the factors.
387  if (! L_solver_.is_null ()) {
388  L_solver_->setMatrix (Teuchos::null);
389  }
390  if (! U_solver_.is_null ()) {
391  U_solver_->setMatrix (Teuchos::null);
392  }
393 
394  L_ = Teuchos::null;
395  U_ = Teuchos::null;
396  A_ = A;
397  }
398 }
399 
400 
401 template<class MatrixType>
403 {
404  Teuchos::Time timer ("ILUT::initialize");
405  {
406  Teuchos::TimeMonitor timeMon (timer);
407 
408  // Check that the matrix is nonnull.
409  TEUCHOS_TEST_FOR_EXCEPTION(
410  A_.is_null (), std::runtime_error, "Ifpack2::ILUT::initialize: "
411  "The matrix to precondition is null. Please call setMatrix() with a "
412  "nonnull input before calling this method.");
413 
414  // Clear any previous computations.
415  IsInitialized_ = false;
416  IsComputed_ = false;
417  A_local_ = Teuchos::null;
418  L_ = Teuchos::null;
419  U_ = Teuchos::null;
420 
421  A_local_ = makeLocalFilter (A_); // Compute the local filter.
422 
423  IsInitialized_ = true;
424  ++NumInitialize_;
425  }
426  InitializeTime_ += timer.totalElapsedTime ();
427 }
428 
429 
430 template<typename ScalarType>
431 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
432 scalar_mag (const ScalarType& s)
433 {
434  return Teuchos::ScalarTraits<ScalarType>::magnitude(s);
435 }
436 
437 
438 template<class MatrixType>
440 {
441  using Teuchos::Array;
442  using Teuchos::ArrayRCP;
443  using Teuchos::ArrayView;
444  using Teuchos::as;
445  using Teuchos::rcp;
446  using Teuchos::reduceAll;
447 
448  //--------------------------------------------------------------------------
449  // Ifpack2::ILUT is a translation of the Aztec ILUT implementation. The Aztec
450  // ILUT implementation was written by Ray Tuminaro.
451  //
452  // This isn't an exact translation of the Aztec ILUT algorithm, for the
453  // following reasons:
454  // 1. Minor differences result from the fact that Aztec factors a MSR format
455  // matrix in place, while the code below factors an input CrsMatrix which
456  // remains untouched and stores the resulting factors in separate L and U
457  // CrsMatrix objects.
458  // Also, the Aztec code begins by shifting the matrix pointers back
459  // by one, and the pointer contents back by one, and then using 1-based
460  // Fortran-style indexing in the algorithm. This Ifpack2 code uses C-style
461  // 0-based indexing throughout.
462  // 2. Aztec stores the inverse of the diagonal of U. This Ifpack2 code
463  // stores the non-inverted diagonal in U.
464  // The triangular solves (in Ifpack2::ILUT::apply()) are performed by
465  // calling the Tpetra::CrsMatrix::solve method on the L and U objects, and
466  // this requires U to contain the non-inverted diagonal.
467  //
468  // ABW.
469  //--------------------------------------------------------------------------
470 
471  // Don't count initialization in the compute() time.
472  if (! isInitialized ()) {
473  initialize ();
474  }
475 
476  Teuchos::Time timer ("ILUT::compute");
477  { // Timer scope for timing compute()
478  Teuchos::TimeMonitor timeMon (timer, true);
479  const scalar_type zero = STS::zero ();
480  const scalar_type one = STS::one ();
481 
482  const local_ordinal_type myNumRows = A_local_->getNodeNumRows ();
483  L_ = rcp (new crs_matrix_type (A_local_->getRowMap (), A_local_->getColMap (), 0));
484  U_ = rcp (new crs_matrix_type (A_local_->getRowMap (), A_local_->getColMap (), 0));
485 
486  // CGB: note, this caching approach may not be necessary anymore
487  // We will store ArrayView objects that are views of the rows of U, so that
488  // we don't have to repeatedly retrieve the view for each row. These will
489  // be populated row by row as the factorization proceeds.
490  Array<ArrayView<const local_ordinal_type> > Uindices (myNumRows);
491  Array<ArrayView<const scalar_type> > Ucoefs (myNumRows);
492 
493  // If this macro is defined, files containing the L and U factors
494  // will be written. DON'T CHECK IN THE CODE WITH THIS MACRO ENABLED!!!
495  // #define IFPACK2_WRITE_FACTORS
496 #ifdef IFPACK2_WRITE_FACTORS
497  std::ofstream ofsL("L.tif.mtx", std::ios::out);
498  std::ofstream ofsU("U.tif.mtx", std::ios::out);
499 #endif
500 
501  // The code here uses double for fill calculations, even though
502  // the fill level is actually an integer. The point is to avoid
503  // rounding and overflow for integer calculations. If int is <=
504  // 32 bits, it can never overflow double, so this cast is always
505  // OK as long as int has <= 32 bits.
506 
507  // Calculate how much fill will be allowed in addition to the
508  // space that corresponds to the input matrix entries.
509  double local_nnz = static_cast<double> (A_local_->getNodeNumEntries ());
510  double fill;
511  {
512  const double fillLevel = as<double> (getLevelOfFill ());
513  fill = ((fillLevel - 1) * local_nnz) / (2 * myNumRows);
514  }
515 
516  // std::ceil gives the smallest integer larger than the argument.
517  // this may give a slightly different result than Aztec's fill value in
518  // some cases.
519  double fill_ceil=std::ceil(fill);
520 
521  // Similarly to Aztec, we will allow the same amount of fill for each
522  // row, half in L and half in U.
523  size_type fillL = static_cast<size_type>(fill_ceil);
524  size_type fillU = static_cast<size_type>(fill_ceil);
525 
526  Array<scalar_type> InvDiagU (myNumRows, zero);
527 
528  Array<local_ordinal_type> tmp_idx;
529  Array<scalar_type> tmpv;
530 
531  enum { UNUSED, ORIG, FILL };
532  local_ordinal_type max_col = myNumRows;
533 
534  Array<int> pattern(max_col, UNUSED);
535  Array<scalar_type> cur_row(max_col, zero);
536  Array<magnitude_type> unorm(max_col);
537  magnitude_type rownorm;
538  Array<local_ordinal_type> L_cols_heap;
539  Array<local_ordinal_type> U_cols;
540  Array<local_ordinal_type> L_vals_heap;
541  Array<local_ordinal_type> U_vals_heap;
542 
543  // A comparison object which will be used to create 'heaps' of indices
544  // that are ordered according to the corresponding values in the
545  // 'cur_row' array.
546  greater_indirect<scalar_type,local_ordinal_type> vals_comp(cur_row);
547 
548  // =================== //
549  // start factorization //
550  // =================== //
551 
552  ArrayRCP<local_ordinal_type> ColIndicesARCP;
553  ArrayRCP<scalar_type> ColValuesARCP;
554  if (! A_local_->supportsRowViews ()) {
555  const size_t maxnz = A_local_->getNodeMaxNumRowEntries ();
556  ColIndicesARCP.resize (maxnz);
557  ColValuesARCP.resize (maxnz);
558  }
559 
560  for (local_ordinal_type row_i = 0 ; row_i < myNumRows ; ++row_i) {
561  ArrayView<const local_ordinal_type> ColIndicesA;
562  ArrayView<const scalar_type> ColValuesA;
563  size_t RowNnz;
564 
565  if (A_local_->supportsRowViews ()) {
566  A_local_->getLocalRowView (row_i, ColIndicesA, ColValuesA);
567  RowNnz = ColIndicesA.size ();
568  }
569  else {
570  A_local_->getLocalRowCopy (row_i, ColIndicesARCP (), ColValuesARCP (), RowNnz);
571  ColIndicesA = ColIndicesARCP (0, RowNnz);
572  ColValuesA = ColValuesARCP (0, RowNnz);
573  }
574 
575  // Always include the diagonal in the U factor. The value should get
576  // set in the next loop below.
577  U_cols.push_back(row_i);
578  cur_row[row_i] = zero;
579  pattern[row_i] = ORIG;
580 
581  size_type L_cols_heaplen = 0;
582  rownorm = STM::zero ();
583  for (size_t i = 0; i < RowNnz; ++i) {
584  if (ColIndicesA[i] < myNumRows) {
585  if (ColIndicesA[i] < row_i) {
586  add_to_heap(ColIndicesA[i], L_cols_heap, L_cols_heaplen);
587  }
588  else if (ColIndicesA[i] > row_i) {
589  U_cols.push_back(ColIndicesA[i]);
590  }
591 
592  cur_row[ColIndicesA[i]] = ColValuesA[i];
593  pattern[ColIndicesA[i]] = ORIG;
594  rownorm += scalar_mag(ColValuesA[i]);
595  }
596  }
597 
598  // Alter the diagonal according to the absolute-threshold and
599  // relative-threshold values. If not set, those values default
600  // to zero and one respectively.
601  const magnitude_type rthresh = getRelativeThreshold();
602  const scalar_type& v = cur_row[row_i];
603  cur_row[row_i] = as<scalar_type> (getAbsoluteThreshold() * IFPACK2_SGN(v)) + rthresh*v;
604 
605  size_type orig_U_len = U_cols.size();
606  RowNnz = L_cols_heap.size() + orig_U_len;
607  rownorm = getDropTolerance() * rownorm/RowNnz;
608 
609  // The following while loop corresponds to the 'L30' goto's in Aztec.
610  size_type L_vals_heaplen = 0;
611  while (L_cols_heaplen > 0) {
612  local_ordinal_type row_k = L_cols_heap.front();
613 
614  scalar_type multiplier = cur_row[row_k] * InvDiagU[row_k];
615  cur_row[row_k] = multiplier;
616  magnitude_type mag_mult = scalar_mag(multiplier);
617  if (mag_mult*unorm[row_k] < rownorm) {
618  pattern[row_k] = UNUSED;
619  rm_heap_root(L_cols_heap, L_cols_heaplen);
620  continue;
621  }
622  if (pattern[row_k] != ORIG) {
623  if (L_vals_heaplen < fillL) {
624  add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
625  }
626  else if (L_vals_heaplen==0 ||
627  mag_mult < scalar_mag(cur_row[L_vals_heap.front()])) {
628  pattern[row_k] = UNUSED;
629  rm_heap_root(L_cols_heap, L_cols_heaplen);
630  continue;
631  }
632  else {
633  pattern[L_vals_heap.front()] = UNUSED;
634  rm_heap_root(L_vals_heap, L_vals_heaplen, vals_comp);
635  add_to_heap(row_k, L_vals_heap, L_vals_heaplen, vals_comp);
636  }
637  }
638 
639  /* Reduce current row */
640 
641  ArrayView<const local_ordinal_type>& ColIndicesU = Uindices[row_k];
642  ArrayView<const scalar_type>& ColValuesU = Ucoefs[row_k];
643  size_type ColNnzU = ColIndicesU.size();
644 
645  for(size_type j=0; j<ColNnzU; ++j) {
646  if (ColIndicesU[j] > row_k) {
647  scalar_type tmp = multiplier * ColValuesU[j];
648  local_ordinal_type col_j = ColIndicesU[j];
649  if (pattern[col_j] != UNUSED) {
650  cur_row[col_j] -= tmp;
651  }
652  else if (scalar_mag(tmp) > rownorm) {
653  cur_row[col_j] = -tmp;
654  pattern[col_j] = FILL;
655  if (col_j > row_i) {
656  U_cols.push_back(col_j);
657  }
658  else {
659  add_to_heap(col_j, L_cols_heap, L_cols_heaplen);
660  }
661  }
662  }
663  }
664 
665  rm_heap_root(L_cols_heap, L_cols_heaplen);
666  }//end of while(L_cols_heaplen) loop
667 
668 
669  // Put indices and values for L into arrays and then into the L_ matrix.
670 
671  // first, the original entries from the L section of A:
672  for (size_type i = 0; i < ColIndicesA.size (); ++i) {
673  if (ColIndicesA[i] < row_i) {
674  tmp_idx.push_back(ColIndicesA[i]);
675  tmpv.push_back(cur_row[ColIndicesA[i]]);
676  pattern[ColIndicesA[i]] = UNUSED;
677  }
678  }
679 
680  // next, the L entries resulting from fill:
681  for (size_type j = 0; j < L_vals_heaplen; ++j) {
682  tmp_idx.push_back(L_vals_heap[j]);
683  tmpv.push_back(cur_row[L_vals_heap[j]]);
684  pattern[L_vals_heap[j]] = UNUSED;
685  }
686 
687  // L has a one on the diagonal, but we don't explicitly store
688  // it. If we don't store it, then the kernel which performs the
689  // triangular solve can assume a unit diagonal, take a short-cut
690  // and perform faster.
691 
692  L_->insertLocalValues (row_i, tmp_idx (), tmpv ());
693 #ifdef IFPACK2_WRITE_FACTORS
694  for (size_type ii = 0; ii < tmp_idx.size (); ++ii) {
695  ofsL << row_i << " " << tmp_idx[ii] << " " << tmpv[ii] << std::endl;
696  }
697 #endif
698 
699  tmp_idx.clear();
700  tmpv.clear();
701 
702  // Pick out the diagonal element, store its reciprocal.
703  if (cur_row[row_i] == zero) {
704  std::cerr << "Ifpack2::ILUT::Compute: zero pivot encountered! Replacing with rownorm and continuing...(You may need to set the parameter 'fact: absolute threshold'.)" << std::endl;
705  cur_row[row_i] = rownorm;
706  }
707  InvDiagU[row_i] = one / cur_row[row_i];
708 
709  // Non-inverted diagonal is stored for U:
710  tmp_idx.push_back(row_i);
711  tmpv.push_back(cur_row[row_i]);
712  unorm[row_i] = scalar_mag(cur_row[row_i]);
713  pattern[row_i] = UNUSED;
714 
715  // Now put indices and values for U into arrays and then into the U_ matrix.
716  // The first entry in U_cols is the diagonal, which we just handled, so we'll
717  // start our loop at j=1.
718 
719  size_type U_vals_heaplen = 0;
720  for(size_type j=1; j<U_cols.size(); ++j) {
721  local_ordinal_type col = U_cols[j];
722  if (pattern[col] != ORIG) {
723  if (U_vals_heaplen < fillU) {
724  add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
725  }
726  else if (U_vals_heaplen!=0 && scalar_mag(cur_row[col]) >
727  scalar_mag(cur_row[U_vals_heap.front()])) {
728  rm_heap_root(U_vals_heap, U_vals_heaplen, vals_comp);
729  add_to_heap(col, U_vals_heap, U_vals_heaplen, vals_comp);
730  }
731  }
732  else {
733  tmp_idx.push_back(col);
734  tmpv.push_back(cur_row[col]);
735  unorm[row_i] += scalar_mag(cur_row[col]);
736  }
737  pattern[col] = UNUSED;
738  }
739 
740  for(size_type j=0; j<U_vals_heaplen; ++j) {
741  tmp_idx.push_back(U_vals_heap[j]);
742  tmpv.push_back(cur_row[U_vals_heap[j]]);
743  unorm[row_i] += scalar_mag(cur_row[U_vals_heap[j]]);
744  }
745 
746  unorm[row_i] /= (orig_U_len + U_vals_heaplen);
747 
748  U_->insertLocalValues(row_i, tmp_idx(), tmpv() );
749 #ifdef IFPACK2_WRITE_FACTORS
750  for(int ii=0; ii<tmp_idx.size(); ++ii) {
751  ofsU <<row_i<< " " <<tmp_idx[ii]<< " " <<tmpv[ii]<< std::endl;
752  }
753 #endif
754  tmp_idx.clear();
755  tmpv.clear();
756 
757  U_->getLocalRowView(row_i, Uindices[row_i], Ucoefs[row_i] );
758 
759  L_cols_heap.clear();
760  U_cols.clear();
761  L_vals_heap.clear();
762  U_vals_heap.clear();
763  } // end of for(row_i) loop
764 
765  // FIXME (mfh 03 Apr 2013) Do we need to supply a domain and range Map?
766  L_->fillComplete();
767  U_->fillComplete();
768 
769  L_solver_ = Teuchos::rcp (new LocalSparseTriangularSolver<row_matrix_type> (L_));
770  L_solver_->initialize ();
771  L_solver_->compute ();
772 
773  U_solver_ = Teuchos::rcp (new LocalSparseTriangularSolver<row_matrix_type> (U_));
774  U_solver_->initialize ();
775  U_solver_->compute ();
776  }
777  ComputeTime_ += timer.totalElapsedTime ();
778  IsComputed_ = true;
779  ++NumCompute_;
780 }
781 
782 
783 template <class MatrixType>
785 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
786  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
787  Teuchos::ETransp mode,
788  scalar_type alpha,
789  scalar_type beta) const
790 {
791  using Teuchos::RCP;
792  using Teuchos::rcp;
793  using Teuchos::rcpFromRef;
794  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> MV;
795 
796  Teuchos::Time timer ("ILUT::apply");
797  { // Timer scope for timing apply()
798  Teuchos::TimeMonitor timeMon (timer, true);
799 
800  TEUCHOS_TEST_FOR_EXCEPTION(
801  ! isComputed (), std::runtime_error,
802  "Ifpack2::ILUT::apply: You must call compute() to compute the incomplete "
803  "factorization, before calling apply().");
804 
805  TEUCHOS_TEST_FOR_EXCEPTION(
806  X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
807  "Ifpack2::ILUT::apply: X and Y must have the same number of columns. "
808  "X has " << X.getNumVectors () << " columns, but Y has "
809  << Y.getNumVectors () << " columns.");
810 
811  if (alpha == Teuchos::ScalarTraits<scalar_type>::zero ()) {
812  // alpha == 0, so we don't need to apply the operator.
813  //
814  // The special case for beta == 0 ensures that if Y contains Inf
815  // or NaN values, we replace them with 0 (following BLAS
816  // convention), rather than multiplying them by 0 to get NaN.
817  if (beta == Teuchos::ScalarTraits<scalar_type>::zero ()) {
818  Y.putScalar (beta);
819  } else {
820  Y.scale (beta);
821  }
822  return;
823  }
824 
825  // If beta != 0, create a temporary multivector Y_temp to hold the
826  // contents of alpha*M^{-1}*X. Otherwise, alias Y_temp to Y.
827  RCP<MV> Y_temp;
828  if (beta == Teuchos::ScalarTraits<scalar_type>::zero ()) {
829  Y_temp = rcpFromRef (Y);
830  } else {
831  Y_temp = rcp (new MV (Y.getMap (), Y.getNumVectors ()));
832  }
833 
834  // If X and Y are pointing to the same memory location, create an
835  // auxiliary vector, X_temp, so that we don't clobber the input
836  // when computing the output. Otherwise, alias X_temp to X.
837  RCP<const MV> X_temp;
838  {
839  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
840  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
841  if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
842  X_temp = rcp (new MV (X, Teuchos::Copy));
843  } else {
844  X_temp = rcpFromRef (X);
845  }
846  }
847 
848  // Create a temporary multivector Y_mid to hold the intermediate
849  // between the L and U (or U and L, for the transpose or conjugate
850  // transpose case) solves.
851  RCP<MV> Y_mid = rcp (new MV (Y.getMap (), Y.getNumVectors ()));
852 
853  if (mode == Teuchos::NO_TRANS) { // Solve L U Y = X
854  L_solver_->apply (*X_temp, *Y_mid, mode);
855 
856  // FIXME (mfh 20 Aug 2013) Is it OK to use Y_temp for both the
857  // input and the output?
858 
859  U_solver_->apply (*Y_mid, *Y_temp, mode);
860  }
861  else { // Solve U^* L^* Y = X
862  U_solver_->apply (*X_temp, *Y_mid, mode);
863 
864  // FIXME (mfh 20 Aug 2013) Is it OK to use Y_temp for both the
865  // input and the output?
866 
867  L_solver_->apply (*Y_mid, *Y_temp, mode);
868  }
869 
870  if (beta == Teuchos::ScalarTraits<scalar_type>::zero ()) {
871  Y.scale (alpha);
872  } else { // beta != 0
873  Y.update (alpha, *Y_temp, beta);
874  }
875  }
876  ++NumApply_;
877  ApplyTime_ += timer.totalElapsedTime ();
878 }
879 
880 
881 template <class MatrixType>
882 std::string ILUT<MatrixType>::description () const
883 {
884  std::ostringstream os;
885 
886  // Output is a valid YAML dictionary in flow style. If you don't
887  // like everything on a single line, you should call describe()
888  // instead.
889  os << "\"Ifpack2::ILUT\": {";
890  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
891  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
892 
893  os << "Level-of-fill: " << getLevelOfFill() << ", "
894  << "absolute threshold: " << getAbsoluteThreshold() << ", "
895  << "relative threshold: " << getRelativeThreshold() << ", "
896  << "relaxation value: " << getRelaxValue() << ", ";
897 
898  if (A_.is_null ()) {
899  os << "Matrix: null";
900  }
901  else {
902  os << "Global matrix dimensions: ["
903  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]"
904  << ", Global nnz: " << A_->getGlobalNumEntries();
905  }
906 
907  os << "}";
908  return os.str ();
909 }
910 
911 
912 template <class MatrixType>
913 void
915 describe (Teuchos::FancyOStream& out,
916  const Teuchos::EVerbosityLevel verbLevel) const
917 {
918  using Teuchos::Comm;
919  using Teuchos::OSTab;
920  using Teuchos::RCP;
921  using Teuchos::TypeNameTraits;
922  using std::endl;
923  using Teuchos::VERB_DEFAULT;
924  using Teuchos::VERB_NONE;
925  using Teuchos::VERB_LOW;
926  using Teuchos::VERB_MEDIUM;
927  using Teuchos::VERB_HIGH;
928  using Teuchos::VERB_EXTREME;
929 
930  const Teuchos::EVerbosityLevel vl =
931  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
932  OSTab tab0 (out);
933 
934  if (vl > VERB_NONE) {
935  out << "\"Ifpack2::ILUT\":" << endl;
936  OSTab tab1 (out);
937  out << "MatrixType: " << TypeNameTraits<MatrixType>::name () << endl;
938  if (this->getObjectLabel () != "") {
939  out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
940  }
941  out << "Initialized: " << (isInitialized () ? "true" : "false")
942  << endl
943  << "Computed: " << (isComputed () ? "true" : "false")
944  << endl
945  << "Level of fill: " << getLevelOfFill () << endl
946  << "Absolute threshold: " << getAbsoluteThreshold () << endl
947  << "Relative threshold: " << getRelativeThreshold () << endl
948  << "Relax value: " << getRelaxValue () << endl;
949 
950  if (isComputed () && vl >= VERB_HIGH) {
951  const double fillFraction =
952  (double) getGlobalNumEntries () / (double) A_->getGlobalNumEntries ();
953  const double nnzToRows =
954  (double) getGlobalNumEntries () / (double) U_->getGlobalNumRows ();
955 
956  out << "Dimensions of L: [" << L_->getGlobalNumRows () << ", "
957  << L_->getGlobalNumRows () << "]" << endl
958  << "Dimensions of U: [" << U_->getGlobalNumRows () << ", "
959  << U_->getGlobalNumRows () << "]" << endl
960  << "Number of nonzeros in factors: " << getGlobalNumEntries () << endl
961  << "Fill fraction of factors over A: " << fillFraction << endl
962  << "Ratio of nonzeros to rows: " << nnzToRows << endl;
963  }
964 
965  out << "Number of initialize calls: " << getNumInitialize () << endl
966  << "Number of compute calls: " << getNumCompute () << endl
967  << "Number of apply calls: " << getNumApply () << endl
968  << "Total time in seconds for initialize: " << getInitializeTime () << endl
969  << "Total time in seconds for compute: " << getComputeTime () << endl
970  << "Total time in seconds for apply: " << getApplyTime () << endl;
971 
972  out << "Local matrix:" << endl;
973  A_local_->describe (out, vl);
974  }
975 }
976 
977 template <class MatrixType>
978 Teuchos::RCP<const typename ILUT<MatrixType>::row_matrix_type>
979 ILUT<MatrixType>::makeLocalFilter (const Teuchos::RCP<const row_matrix_type>& A)
980 {
981  if (A->getComm ()->getSize () > 1) {
982  return Teuchos::rcp (new LocalFilter<row_matrix_type> (A));
983  } else {
984  return A;
985  }
986 }
987 
988 }//namespace Ifpack2
989 
990 
991 // FIXME (mfh 16 Sep 2014) We should really only use RowMatrix here!
992 // There's no need to instantiate for CrsMatrix too. All Ifpack2
993 // preconditioners can and should do dynamic casts if they need a type
994 // more specific than RowMatrix.
995 
996 #define IFPACK2_ILUT_INSTANT(S,LO,GO,N) \
997  template class Ifpack2::ILUT< Tpetra::RowMatrix<S, LO, GO, N> >;
998 
999 #endif /* IFPACK2_ILUT_DEF_HPP */
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_ILUT_def.hpp:320
ILUT(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_ILUT_def.hpp:116
virtual ~ILUT()
Destructor.
Definition: Ifpack2_ILUT_def.hpp:134
global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_ILUT_def.hpp:350
bool hasTransposeApply() const
Whether this object&#39;s apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_ILUT_def.hpp:308
size_t getNodeNumEntries() const
Returns the number of nonzero entries in the local graph.
Definition: Ifpack2_ILUT_def.hpp:356
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_ILUT_decl.hpp:118
void initialize()
Clear any previously computed factors.
Definition: Ifpack2_ILUT_def.hpp:402
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_ILUT_def.hpp:326
int getNumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack2_ILUT_def.hpp:314
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:91
Teuchos::RCP< const map_type > getRangeMap() const
Tpetra::Map representing the range of this operator.
Definition: Ifpack2_ILUT_def.hpp:297
void rm_heap_root(Teuchos::Array< Ordinal > &heap, SizeType &heap_len)
Definition: Ifpack2_Heap.hpp:92
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_ILUT_def.hpp:332
"Preconditioner" that solves local sparse triangular systems.
Definition: Ifpack2_LocalSparseTriangularSolver_decl.hpp:83
Teuchos::RCP< const map_type > getDomainMap() const
Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_ILUT_def.hpp:285
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_ILUT_decl.hpp:106
void compute()
Compute factors L and U using the specified diagonal perturbation thresholds and relaxation parameter...
Definition: Ifpack2_ILUT_def.hpp:439
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:362
Teuchos::RCP< const row_matrix_type > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_ILUT_def.hpp:278
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_ILUT_def.hpp:882
Definition: Ifpack2_Container.hpp:761
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:915
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_ILUT_def.hpp:338
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:126
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the input matrix&#39;s communicator.
Definition: Ifpack2_ILUT_def.hpp:267
void add_to_heap(const Ordinal &idx, Teuchos::Array< Ordinal > &heap, SizeType &heap_len)
Definition: Ifpack2_Heap.hpp:70
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
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:785
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_ILUT_def.hpp:344
void setParameters(const Teuchos::ParameterList &params)
Set preconditioner parameters.
Definition: Ifpack2_ILUT_def.hpp:138
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_ILUT_decl.hpp:109