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