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