Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Experimental_RBILUK_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_EXPERIMENTAL_CRSRBILUK_DEF_HPP
44 #define IFPACK2_EXPERIMENTAL_CRSRBILUK_DEF_HPP
45 
46 #include "Tpetra_Experimental_BlockMultiVector.hpp"
47 #include "Tpetra_Experimental_BlockView.hpp"
48 #include "Ifpack2_OverlappingRowMatrix.hpp"
49 #include "Ifpack2_LocalFilter.hpp"
50 #include "Ifpack2_Utilities.hpp"
51 #include "Ifpack2_RILUK.hpp"
52 
53 //#define IFPACK2_RBILUK_INITIAL
54 #define IFPACK2_RBILUK_INITIAL_NOKK
55 
56 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
57 #include "KokkosBatched_Gemm_Decl.hpp"
58 #include "KokkosBatched_Gemm_Serial_Impl.hpp"
59 #include "KokkosBatched_Util.hpp"
60 #endif
61 
62 namespace Ifpack2 {
63 
64 namespace Experimental {
65 
66 template<class MatrixType>
68  : RILUK<row_matrix_type>(Teuchos::rcp_dynamic_cast<const row_matrix_type>(Matrix_in) ),
69  A_(Matrix_in),
70  A_block_(Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(Matrix_in))
71 {}
72 
73 template<class MatrixType>
75  : RILUK<row_matrix_type>(Teuchos::rcp_dynamic_cast<const row_matrix_type>(Matrix_in) ),
76  A_block_(Matrix_in)
77 {}
78 
79 
80 template<class MatrixType>
82 
83 
84 template<class MatrixType>
85 void
87 {
88  // FIXME (mfh 04 Nov 2015) What about A_? When does that get (re)set?
89 
90  // It's legal for A to be null; in that case, you may not call
91  // initialize() until calling setMatrix() with a nonnull input.
92  // Regardless, setting the matrix invalidates any previous
93  // factorization.
94  if (A.getRawPtr () != A_block_.getRawPtr ())
95  {
96  this->isAllocated_ = false;
97  this->isInitialized_ = false;
98  this->isComputed_ = false;
99  this->Graph_ = Teuchos::null;
100  L_block_ = Teuchos::null;
101  U_block_ = Teuchos::null;
102  D_block_ = Teuchos::null;
103  A_block_ = A;
104  }
105 }
106 
107 
108 
109 template<class MatrixType>
110 const typename RBILUK<MatrixType>::block_crs_matrix_type&
112 {
114  L_block_.is_null (), std::runtime_error, "Ifpack2::RILUK::getL: The L factor "
115  "is null. Please call initialize() (and preferably also compute()) "
116  "before calling this method. If the input matrix has not yet been set, "
117  "you must first call setMatrix() with a nonnull input matrix before you "
118  "may call initialize() or compute().");
119  return *L_block_;
120 }
121 
122 
123 template<class MatrixType>
124 const typename RBILUK<MatrixType>::block_crs_matrix_type&
126 {
128  D_block_.is_null (), std::runtime_error, "Ifpack2::RILUK::getD: The D factor "
129  "(of diagonal entries) is null. Please call initialize() (and "
130  "preferably also compute()) before calling this method. If the input "
131  "matrix has not yet been set, you must first call setMatrix() with a "
132  "nonnull input matrix before you may call initialize() or compute().");
133  return *D_block_;
134 }
135 
136 
137 template<class MatrixType>
138 const typename RBILUK<MatrixType>::block_crs_matrix_type&
140 {
142  U_block_.is_null (), std::runtime_error, "Ifpack2::RILUK::getU: The U factor "
143  "is null. Please call initialize() (and preferably also compute()) "
144  "before calling this method. If the input matrix has not yet been set, "
145  "you must first call setMatrix() with a nonnull input matrix before you "
146  "may call initialize() or compute().");
147  return *U_block_;
148 }
149 
150 template<class MatrixType>
152 {
153  using Teuchos::null;
154  using Teuchos::rcp;
155 
156  if (! this->isAllocated_) {
157  // Deallocate any existing storage. This avoids storing 2x
158  // memory, since RCP op= does not deallocate until after the
159  // assignment.
160  this->L_ = null;
161  this->U_ = null;
162  this->D_ = null;
163  L_block_ = null;
164  U_block_ = null;
165  D_block_ = null;
166 
167  // Allocate Matrix using ILUK graphs
168  L_block_ = rcp(new block_crs_matrix_type(*this->Graph_->getL_Graph (), blockSize_) );
169  U_block_ = rcp(new block_crs_matrix_type(*this->Graph_->getU_Graph (), blockSize_) );
170  D_block_ = rcp(new block_crs_matrix_type(*(Ifpack2::Details::computeDiagonalGraph(*(this->Graph_->getOverlapGraph()))),
171  blockSize_) );
172  L_block_->setAllToScalar (STM::zero ()); // Zero out L and U matrices
173  U_block_->setAllToScalar (STM::zero ());
174  D_block_->setAllToScalar (STM::zero ());
175 
176  }
177  this->isAllocated_ = true;
178 }
179 
180 template<class MatrixType>
183  return A_block_;
184 }
185 
186 template<class MatrixType>
188 {
189  using Teuchos::RCP;
190  using Teuchos::rcp;
191  using Teuchos::rcp_dynamic_cast;
192  const char prefix[] = "Ifpack2::Experimental::RBILUK::initialize: ";
193 
194  // FIXME (mfh 04 Nov 2015) Apparently it's OK for A_ to be null.
195  // That probably means that this preconditioner was created with a
196  // BlockCrsMatrix directly, so it doesn't need the LocalFilter.
197 
198  // TEUCHOS_TEST_FOR_EXCEPTION
199  // (A_.is_null (), std::runtime_error, prefix << "The matrix (A_, the "
200  // "RowMatrix) is null. Please call setMatrix() with a nonnull input "
201  // "before calling this method.");
202 
203  if (A_block_.is_null ()) {
204  // FIXME (mfh 04 Nov 2015) Why does the input have to be a
205  // LocalFilter? Why can't we just take a regular matrix, and
206  // apply a LocalFilter only if necessary, like other "local"
207  // Ifpack2 preconditioners already do?
208  RCP<const LocalFilter<row_matrix_type> > filteredA =
209  rcp_dynamic_cast<const LocalFilter<row_matrix_type> >(A_);
211  (filteredA.is_null (), std::runtime_error, prefix <<
212  "Cannot cast to filtered matrix.");
213  RCP<const OverlappingRowMatrix<row_matrix_type> > overlappedA =
214  rcp_dynamic_cast<const OverlappingRowMatrix<row_matrix_type> > (filteredA->getUnderlyingMatrix ());
215  if (! overlappedA.is_null ()) {
216  A_block_ = rcp_dynamic_cast<const block_crs_matrix_type>(overlappedA->getUnderlyingMatrix());
217  } else {
218  //If there is no overlap, filteredA could be the block CRS matrix
219  A_block_ = rcp_dynamic_cast<const block_crs_matrix_type>(filteredA->getUnderlyingMatrix());
220  }
221  }
222 
224  (A_block_.is_null (), std::runtime_error, prefix << "The matrix (A_block_, "
225  "the BlockCrsMatrix) is null. Please call setMatrix() with a nonnull "
226  "input before calling this method.");
228  (! A_block_->isFillComplete (), std::runtime_error, prefix << "The matrix "
229  "(A_block_, the BlockCrsMatrix) is not fill complete. You may not invoke "
230  "initialize() or compute() with this matrix until the matrix is fill "
231  "complete. Note: BlockCrsMatrix is fill complete if and only if its "
232  "underlying graph is fill complete.");
233 
234  blockSize_ = A_block_->getBlockSize();
235 
236  Teuchos::Time timer ("RBILUK::initialize");
237  { // Start timing
238  Teuchos::TimeMonitor timeMon (timer);
239 
240  // Calling initialize() means that the user asserts that the graph
241  // of the sparse matrix may have changed. We must not just reuse
242  // the previous graph in that case.
243  //
244  // Regarding setting isAllocated_ to false: Eventually, we may want
245  // some kind of clever memory reuse strategy, but it's always
246  // correct just to blow everything away and start over.
247  this->isInitialized_ = false;
248  this->isAllocated_ = false;
249  this->isComputed_ = false;
250  this->Graph_ = Teuchos::null;
251 
252  typedef Tpetra::CrsGraph<local_ordinal_type,
254  node_type> crs_graph_type;
255 
256  RCP<const crs_graph_type> matrixCrsGraph = Teuchos::rcpFromRef(A_block_->getCrsGraph() );
257  this->Graph_ = rcp (new Ifpack2::IlukGraph<crs_graph_type> (matrixCrsGraph,
258  this->LevelOfFill_, 0));
259 
260  this->Graph_->initialize ();
261  allocate_L_and_U_blocks ();
262 #ifdef IFPACK2_RBILUK_INITIAL
263  initAllValues (*A_block_);
264 #endif
265  } // Stop timing
266 
267  this->isInitialized_ = true;
268  this->numInitialize_ += 1;
269  this->initializeTime_ += timer.totalElapsedTime ();
270 }
271 
272 
273 template<class MatrixType>
274 void
276 initAllValues (const block_crs_matrix_type& A)
277 {
278  using Teuchos::RCP;
279  typedef Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> map_type;
280 
281  local_ordinal_type NumIn = 0, NumL = 0, NumU = 0;
282  bool DiagFound = false;
283  size_t NumNonzeroDiags = 0;
284  size_t MaxNumEntries = A.getNodeMaxNumRowEntries();
285  local_ordinal_type blockMatSize = blockSize_*blockSize_;
286 
287  // First check that the local row map ordering is the same as the local portion of the column map.
288  // The extraction of the strictly lower/upper parts of A, as well as the factorization,
289  // implicitly assume that this is the case.
290  Teuchos::ArrayView<const global_ordinal_type> rowGIDs = A.getRowMap()->getNodeElementList();
291  Teuchos::ArrayView<const global_ordinal_type> colGIDs = A.getColMap()->getNodeElementList();
292  bool gidsAreConsistentlyOrdered=true;
293  global_ordinal_type indexOfInconsistentGID=0;
294  for (global_ordinal_type i=0; i<rowGIDs.size(); ++i) {
295  if (rowGIDs[i] != colGIDs[i]) {
296  gidsAreConsistentlyOrdered=false;
297  indexOfInconsistentGID=i;
298  break;
299  }
300  }
301  TEUCHOS_TEST_FOR_EXCEPTION(gidsAreConsistentlyOrdered==false, std::runtime_error,
302  "The ordering of the local GIDs in the row and column maps is not the same"
303  << std::endl << "at index " << indexOfInconsistentGID
304  << ". Consistency is required, as all calculations are done with"
305  << std::endl << "local indexing.");
306 
307  // Allocate temporary space for extracting the strictly
308  // lower and upper parts of the matrix A.
309  Teuchos::Array<local_ordinal_type> LI(MaxNumEntries);
310  Teuchos::Array<local_ordinal_type> UI(MaxNumEntries);
311  Teuchos::Array<scalar_type> LV(MaxNumEntries*blockMatSize);
312  Teuchos::Array<scalar_type> UV(MaxNumEntries*blockMatSize);
313 
314  Teuchos::Array<scalar_type> diagValues(blockMatSize);
315 
316  L_block_->setAllToScalar (STM::zero ()); // Zero out L and U matrices
317  U_block_->setAllToScalar (STM::zero ());
318  D_block_->setAllToScalar (STM::zero ()); // Set diagonal values to zero
319 
320  // NOTE (mfh 27 May 2016) The factorization below occurs entirely on
321  // host, so sync to host first. The const_cast is unfortunate but
322  // is our only option to make this correct.
323 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
324  const_cast<block_crs_matrix_type&> (A).template sync<Kokkos::HostSpace> ();
325  L_block_->template sync<Kokkos::HostSpace> ();
326  U_block_->template sync<Kokkos::HostSpace> ();
327  D_block_->template sync<Kokkos::HostSpace> ();
328  // NOTE (mfh 27 May 2016) We're modifying L, U, and D on host.
329  L_block_->template modify<Kokkos::HostSpace> ();
330  U_block_->template modify<Kokkos::HostSpace> ();
331  D_block_->template modify<Kokkos::HostSpace> ();
332 #else
333  const_cast<block_crs_matrix_type&> (A).sync_host ();
334  L_block_->sync_host ();
335  U_block_->sync_host ();
336  D_block_->sync_host ();
337  // NOTE (mfh 27 May 2016) We're modifying L, U, and D on host.
338  L_block_->modify_host ();
339  U_block_->modify_host ();
340  D_block_->modify_host ();
341 #endif
342 
343  RCP<const map_type> rowMap = L_block_->getRowMap ();
344 
345  // First we copy the user's matrix into L and U, regardless of fill level.
346  // It is important to note that L and U are populated using local indices.
347  // This means that if the row map GIDs are not monotonically increasing
348  // (i.e., permuted or gappy), then the strictly lower (upper) part of the
349  // matrix is not the one that you would get if you based L (U) on GIDs.
350  // This is ok, as the *order* of the GIDs in the rowmap is a better
351  // expression of the user's intent than the GIDs themselves.
352 
353  for (size_t myRow=0; myRow<A.getNodeNumRows(); ++myRow) {
354  local_ordinal_type local_row = myRow;
355 
356  //TODO JJH 4April2014 An optimization is to use getLocalRowView. Not all matrices support this,
357  // we'd need to check via the Tpetra::RowMatrix method supportsRowViews().
358  const local_ordinal_type * InI = 0;
359  scalar_type * InV = 0;
360  A.getLocalRowView(local_row, InI, InV, NumIn);
361 
362  // Split into L and U (we don't assume that indices are ordered).
363 
364  NumL = 0;
365  NumU = 0;
366  DiagFound = false;
367 
368  for (local_ordinal_type j = 0; j < NumIn; ++j) {
369  const local_ordinal_type k = InI[j];
370  const local_ordinal_type blockOffset = blockMatSize*j;
371 
372  if (k == local_row) {
373  DiagFound = true;
374  // Store perturbed diagonal in Tpetra::Vector D_
375  for (local_ordinal_type jj = 0; jj < blockMatSize; ++jj)
376  diagValues[jj] = this->Rthresh_ * InV[blockOffset+jj] + IFPACK2_SGN(InV[blockOffset+jj]) * this->Athresh_;
377  D_block_->replaceLocalValues(local_row, &InI[j], diagValues.getRawPtr(), 1);
378  }
379  else if (k < 0) { // Out of range
381  true, std::runtime_error, "Ifpack2::RILUK::initAllValues: current "
382  "GID k = " << k << " < 0. I'm not sure why this is an error; it is "
383  "probably an artifact of the undocumented assumptions of the "
384  "original implementation (likely copied and pasted from Ifpack). "
385  "Nevertheless, the code I found here insisted on this being an error "
386  "state, so I will throw an exception here.");
387  }
388  else if (k < local_row) {
389  LI[NumL] = k;
390  const local_ordinal_type LBlockOffset = NumL*blockMatSize;
391  for (local_ordinal_type jj = 0; jj < blockMatSize; ++jj)
392  LV[LBlockOffset+jj] = InV[blockOffset+jj];
393  NumL++;
394  }
395  else if (Teuchos::as<size_t>(k) <= rowMap->getNodeNumElements()) {
396  UI[NumU] = k;
397  const local_ordinal_type UBlockOffset = NumU*blockMatSize;
398  for (local_ordinal_type jj = 0; jj < blockMatSize; ++jj)
399  UV[UBlockOffset+jj] = InV[blockOffset+jj];
400  NumU++;
401  }
402  }
403 
404  // Check in things for this row of L and U
405 
406  if (DiagFound) {
407  ++NumNonzeroDiags;
408  } else
409  {
410  for (local_ordinal_type jj = 0; jj < blockSize_; ++jj)
411  diagValues[jj*(blockSize_+1)] = this->Athresh_;
412  D_block_->replaceLocalValues(local_row, &local_row, diagValues.getRawPtr(), 1);
413  }
414 
415  if (NumL) {
416  L_block_->replaceLocalValues(local_row, &LI[0], &LV[0], NumL);
417  }
418 
419  if (NumU) {
420  U_block_->replaceLocalValues(local_row, &UI[0], &UV[0], NumU);
421  }
422  }
423 
424  // NOTE (mfh 27 May 2016) Sync back to device, in case compute()
425  // ever gets a device implementation.
426  {
427  typedef typename block_crs_matrix_type::device_type device_type;
428  const_cast<block_crs_matrix_type&> (A).template sync<device_type> ();
429  L_block_->template sync<device_type> ();
430  U_block_->template sync<device_type> ();
431  D_block_->template sync<device_type> ();
432  }
433  this->isInitialized_ = true;
434 }
435 
436 namespace { // (anonymous)
437 
438 // For a given Kokkos::View type, possibly unmanaged, get the
439 // corresponding managed Kokkos::View type. This is handy for
440 // translating from little_block_type or little_vec_type (both
441 // possibly unmanaged) to their managed versions.
442 template<class LittleBlockType>
443 struct GetManagedView {
444  static_assert (Kokkos::Impl::is_view<LittleBlockType>::value,
445  "LittleBlockType must be a Kokkos::View.");
446  typedef Kokkos::View<typename LittleBlockType::non_const_data_type,
447  typename LittleBlockType::array_layout,
448  typename LittleBlockType::device_type> managed_non_const_type;
449  static_assert (static_cast<int> (managed_non_const_type::rank) ==
450  static_cast<int> (LittleBlockType::rank),
451  "managed_non_const_type::rank != LittleBlock::rank. "
452  "Please report this bug to the Ifpack2 developers.");
453 };
454 
455 } // namespace (anonymous)
456 
457 template<class MatrixType>
459 {
460  typedef impl_scalar_type IST;
461  const char prefix[] = "Ifpack2::Experimental::RBILUK::compute: ";
462 
463  // initialize() checks this too, but it's easier for users if the
464  // error shows them the name of the method that they actually
465  // called, rather than the name of some internally called method.
467  (A_block_.is_null (), std::runtime_error, prefix << "The matrix (A_block_, "
468  "the BlockCrsMatrix) is null. Please call setMatrix() with a nonnull "
469  "input before calling this method.");
471  (! A_block_->isFillComplete (), std::runtime_error, prefix << "The matrix "
472  "(A_block_, the BlockCrsMatrix) is not fill complete. You may not invoke "
473  "initialize() or compute() with this matrix until the matrix is fill "
474  "complete. Note: BlockCrsMatrix is fill complete if and only if its "
475  "underlying graph is fill complete.");
476 
477  if (! this->isInitialized ()) {
478  initialize (); // Don't count this in the compute() time
479  }
480 
481  // NOTE (mfh 27 May 2016) The factorization below occurs entirely on
482  // host, so sync to host first. The const_cast is unfortunate but
483  // is our only option to make this correct.
484  if (! A_block_.is_null ()) {
486  Teuchos::rcp_const_cast<block_crs_matrix_type> (A_block_);
487 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
488  A_nc->template sync<Kokkos::HostSpace> ();
489 #else
490  A_nc->sync_host ();
491 #endif
492  }
493 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
494  L_block_->template sync<Kokkos::HostSpace> ();
495  U_block_->template sync<Kokkos::HostSpace> ();
496  D_block_->template sync<Kokkos::HostSpace> ();
497  // NOTE (mfh 27 May 2016) We're modifying L, U, and D on host.
498  L_block_->template modify<Kokkos::HostSpace> ();
499  U_block_->template modify<Kokkos::HostSpace> ();
500  D_block_->template modify<Kokkos::HostSpace> ();
501 #else
502  L_block_->sync_host ();
503  U_block_->sync_host ();
504  D_block_->sync_host ();
505  // NOTE (mfh 27 May 2016) We're modifying L, U, and D on host.
506  L_block_->modify_host ();
507  U_block_->modify_host ();
508  D_block_->modify_host ();
509 #endif
510 
511  Teuchos::Time timer ("RBILUK::compute");
512  { // Start timing
513  Teuchos::TimeMonitor timeMon (timer);
514  this->isComputed_ = false;
515 
516  // MinMachNum should be officially defined, for now pick something a little
517  // bigger than IEEE underflow value
518 
519 // const scalar_type MinDiagonalValue = STS::rmin ();
520 // const scalar_type MaxDiagonalValue = STS::one () / MinDiagonalValue;
521  initAllValues (*A_block_);
522 
523  size_t NumIn;
524  local_ordinal_type NumL, NumU, NumURead;
525 
526  // Get Maximum Row length
527  const size_t MaxNumEntries =
528  L_block_->getNodeMaxNumRowEntries () + U_block_->getNodeMaxNumRowEntries () + 1;
529 
530  const local_ordinal_type blockMatSize = blockSize_*blockSize_;
531 
532  // FIXME (mfh 08 Nov 2015, 24 May 2016) We need to move away from
533  // expressing these strides explicitly, in order to finish #177
534  // (complete Kokkos-ization of BlockCrsMatrix) thoroughly.
535  const local_ordinal_type rowStride = blockSize_;
536 
537  Teuchos::Array<int> ipiv_teuchos(blockSize_);
538  Kokkos::View<int*, Kokkos::HostSpace,
539  Kokkos::MemoryUnmanaged> ipiv (ipiv_teuchos.getRawPtr (), blockSize_);
540  Teuchos::Array<IST> work_teuchos(blockSize_);
541  Kokkos::View<IST*, Kokkos::HostSpace,
542  Kokkos::MemoryUnmanaged> work (work_teuchos.getRawPtr (), blockSize_);
543 
544  size_t num_cols = U_block_->getColMap()->getNodeNumElements();
545  Teuchos::Array<int> colflag(num_cols);
546 
547  typename GetManagedView<little_block_type>::managed_non_const_type diagModBlock ("diagModBlock", blockSize_, blockSize_);
548  typename GetManagedView<little_block_type>::managed_non_const_type matTmp ("matTmp", blockSize_, blockSize_);
549  typename GetManagedView<little_block_type>::managed_non_const_type multiplier ("multiplier", blockSize_, blockSize_);
550 
551 // Teuchos::ArrayRCP<scalar_type> DV = D_->get1dViewNonConst(); // Get view of diagonal
552 
553  // Now start the factorization.
554 
555  // Need some integer workspace and pointers
556  local_ordinal_type NumUU;
557  for (size_t j = 0; j < num_cols; ++j) {
558  colflag[j] = -1;
559  }
560  Teuchos::Array<local_ordinal_type> InI(MaxNumEntries, 0);
561  Teuchos::Array<scalar_type> InV(MaxNumEntries*blockMatSize,STM::zero());
562 
563  const local_ordinal_type numLocalRows = L_block_->getNodeNumRows ();
564  for (local_ordinal_type local_row = 0; local_row < numLocalRows; ++local_row) {
565 
566  // Fill InV, InI with current row of L, D and U combined
567 
568  NumIn = MaxNumEntries;
569  const local_ordinal_type * colValsL;
570  scalar_type * valsL;
571 
572  L_block_->getLocalRowView(local_row, colValsL, valsL, NumL);
573  for (local_ordinal_type j = 0; j < NumL; ++j)
574  {
575  const local_ordinal_type matOffset = blockMatSize*j;
576  little_block_type lmat((typename little_block_type::value_type*) &valsL[matOffset],blockSize_,rowStride);
577  little_block_type lmatV((typename little_block_type::value_type*) &InV[matOffset],blockSize_,rowStride);
578  //lmatV.assign(lmat);
579  Tpetra::Experimental::COPY (lmat, lmatV);
580  InI[j] = colValsL[j];
581  }
582 
583  little_block_type dmat = D_block_->getLocalBlock(local_row, local_row);
584  little_block_type dmatV((typename little_block_type::value_type*) &InV[NumL*blockMatSize], blockSize_, rowStride);
585  //dmatV.assign(dmat);
586  Tpetra::Experimental::COPY (dmat, dmatV);
587  InI[NumL] = local_row;
588 
589  const local_ordinal_type * colValsU;
590  scalar_type * valsU;
591  U_block_->getLocalRowView(local_row, colValsU, valsU, NumURead);
592  NumU = 0;
593  for (local_ordinal_type j = 0; j < NumURead; ++j)
594  {
595  if (!(colValsU[j] < numLocalRows)) continue;
596  InI[NumL+1+j] = colValsU[j];
597  const local_ordinal_type matOffset = blockMatSize*(NumL+1+j);
598  little_block_type umat((typename little_block_type::value_type*) &valsU[blockMatSize*j], blockSize_, rowStride);
599  little_block_type umatV((typename little_block_type::value_type*) &InV[matOffset], blockSize_, rowStride);
600  //umatV.assign(umat);
601  Tpetra::Experimental::COPY (umat, umatV);
602  NumU += 1;
603  }
604  NumIn = NumL+NumU+1;
605 
606  // Set column flags
607  for (size_t j = 0; j < NumIn; ++j) {
608  colflag[InI[j]] = j;
609  }
610 
611 #ifndef IFPACK2_RBILUK_INITIAL
612  for (local_ordinal_type i = 0; i < blockSize_; ++i)
613  for (local_ordinal_type j = 0; j < blockSize_; ++j){
614  {
615  diagModBlock(i,j) = 0;
616  }
617  }
618 #else
619  scalar_type diagmod = STM::zero (); // Off-diagonal accumulator
620  Kokkos::deep_copy (diagModBlock, diagmod);
621 #endif
622 
623  for (local_ordinal_type jj = 0; jj < NumL; ++jj) {
624  local_ordinal_type j = InI[jj];
625  little_block_type currentVal((typename little_block_type::value_type*) &InV[jj*blockMatSize], blockSize_, rowStride); // current_mults++;
626  //multiplier.assign(currentVal);
627  Tpetra::Experimental::COPY (currentVal, multiplier);
628 
629  const little_block_type dmatInverse = D_block_->getLocalBlock(j,j);
630  // alpha = 1, beta = 0
631 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
632  KokkosBatched::Experimental::SerialGemm
633  <KokkosBatched::Experimental::Trans::NoTranspose,
634  KokkosBatched::Experimental::Trans::NoTranspose,
635  KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
636  (STS::one (), currentVal, dmatInverse, STS::zero (), matTmp);
637 #else
638  Tpetra::Experimental::GEMM ("N", "N", STS::one (), currentVal, dmatInverse,
639  STS::zero (), matTmp);
640 #endif
641  //blockMatOpts.square_matrix_matrix_multiply(reinterpret_cast<impl_scalar_type*> (currentVal.data ()), reinterpret_cast<impl_scalar_type*> (dmatInverse.data ()), reinterpret_cast<impl_scalar_type*> (matTmp.data ()), blockSize_);
642  //currentVal.assign(matTmp);
643  Tpetra::Experimental::COPY (matTmp, currentVal);
644 
645  const local_ordinal_type * UUI;
646  scalar_type * UUV;
647  U_block_->getLocalRowView(j, UUI, UUV, NumUU);
648 
649  if (this->RelaxValue_ == STM::zero ()) {
650  for (local_ordinal_type k = 0; k < NumUU; ++k) {
651  if (!(UUI[k] < numLocalRows)) continue;
652  const int kk = colflag[UUI[k]];
653  if (kk > -1) {
654  little_block_type kkval((typename little_block_type::value_type*) &InV[kk*blockMatSize], blockSize_, rowStride);
655  little_block_type uumat((typename little_block_type::value_type*) &UUV[k*blockMatSize], blockSize_, rowStride);
656 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
657  KokkosBatched::Experimental::SerialGemm
658  <KokkosBatched::Experimental::Trans::NoTranspose,
659  KokkosBatched::Experimental::Trans::NoTranspose,
660  KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
661  ( magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), kkval);
662 #else
663  Tpetra::Experimental::GEMM ("N", "N", magnitude_type(-STM::one ()), multiplier, uumat,
664  STM::one (), kkval);
665 #endif
666  //blockMatOpts.square_matrix_matrix_multiply(reinterpret_cast<impl_scalar_type*> (multiplier.data ()), reinterpret_cast<impl_scalar_type*> (uumat.data ()), reinterpret_cast<impl_scalar_type*> (kkval.data ()), blockSize_, -STM::one(), STM::one());
667  }
668  }
669  }
670  else {
671  for (local_ordinal_type k = 0; k < NumUU; ++k) {
672  if (!(UUI[k] < numLocalRows)) continue;
673  const int kk = colflag[UUI[k]];
674  little_block_type uumat((typename little_block_type::value_type*) &UUV[k*blockMatSize], blockSize_, rowStride);
675  if (kk > -1) {
676  little_block_type kkval((typename little_block_type::value_type*) &InV[kk*blockMatSize], blockSize_, rowStride);
677 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
678  KokkosBatched::Experimental::SerialGemm
679  <KokkosBatched::Experimental::Trans::NoTranspose,
680  KokkosBatched::Experimental::Trans::NoTranspose,
681  KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
682  (magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), kkval);
683 #else
684  Tpetra::Experimental::GEMM ("N", "N", magnitude_type(-STM::one ()), multiplier, uumat,
685  STM::one (), kkval);
686 #endif
687  //blockMatOpts.square_matrix_matrix_multiply(reinterpret_cast<impl_scalar_type*>(multiplier.data ()), reinterpret_cast<impl_scalar_type*>(uumat.data ()), reinterpret_cast<impl_scalar_type*>(kkval.data ()), blockSize_, -STM::one(), STM::one());
688  }
689  else {
690 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
691  KokkosBatched::Experimental::SerialGemm
692  <KokkosBatched::Experimental::Trans::NoTranspose,
693  KokkosBatched::Experimental::Trans::NoTranspose,
694  KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
695  (magnitude_type(-STM::one ()), multiplier, uumat, STM::one (), diagModBlock);
696 #else
697  Tpetra::Experimental::GEMM ("N", "N", magnitude_type(-STM::one ()), multiplier, uumat,
698  STM::one (), diagModBlock);
699 #endif
700  //blockMatOpts.square_matrix_matrix_multiply(reinterpret_cast<impl_scalar_type*>(multiplier.data ()), reinterpret_cast<impl_scalar_type*>(uumat.data ()), reinterpret_cast<impl_scalar_type*>(diagModBlock.data ()), blockSize_, -STM::one(), STM::one());
701  }
702  }
703  }
704  }
705  if (NumL) {
706  // Replace current row of L
707  L_block_->replaceLocalValues (local_row, InI.getRawPtr (), InV.getRawPtr (), NumL);
708  }
709 
710  // dmat.assign(dmatV);
711  Tpetra::Experimental::COPY (dmatV, dmat);
712 
713  if (this->RelaxValue_ != STM::zero ()) {
714  //dmat.update(this->RelaxValue_, diagModBlock);
715  Tpetra::Experimental::AXPY (this->RelaxValue_, diagModBlock, dmat);
716  }
717 
718 // if (STS::magnitude (DV[i]) > STS::magnitude (MaxDiagonalValue)) {
719 // if (STS::real (DV[i]) < STM::zero ()) {
720 // DV[i] = -MinDiagonalValue;
721 // }
722 // else {
723 // DV[i] = MinDiagonalValue;
724 // }
725 // }
726 // else
727  {
728  int lapackInfo = 0;
729  for (int k = 0; k < blockSize_; ++k) {
730  ipiv[k] = 0;
731  }
732 
733  Tpetra::Experimental::GETF2 (dmat, ipiv, lapackInfo);
734  //lapack.GETRF(blockSize_, blockSize_, d_raw, blockSize_, ipiv.getRawPtr(), &lapackInfo);
736  lapackInfo != 0, std::runtime_error, "Ifpack2::Experimental::RBILUK::compute: "
737  "lapackInfo = " << lapackInfo << " which indicates an error in the factorization GETRF.");
738 
739  Tpetra::Experimental::GETRI (dmat, ipiv, work, lapackInfo);
740  //lapack.GETRI(blockSize_, d_raw, blockSize_, ipiv.getRawPtr(), work.getRawPtr(), lwork, &lapackInfo);
742  lapackInfo != 0, std::runtime_error, "Ifpack2::Experimental::RBILUK::compute: "
743  "lapackInfo = " << lapackInfo << " which indicates an error in the matrix inverse GETRI.");
744  }
745 
746  for (local_ordinal_type j = 0; j < NumU; ++j) {
747  little_block_type currentVal((typename little_block_type::value_type*) &InV[(NumL+1+j)*blockMatSize], blockSize_, rowStride); // current_mults++;
748  // scale U by the diagonal inverse
749 #ifndef IFPACK2_RBILUK_INITIAL_NOKK
750  KokkosBatched::Experimental::SerialGemm
751  <KokkosBatched::Experimental::Trans::NoTranspose,
752  KokkosBatched::Experimental::Trans::NoTranspose,
753  KokkosBatched::Experimental::Algo::Gemm::Blocked>::invoke
754  (STS::one (), dmat, currentVal, STS::zero (), matTmp);
755 #else
756  Tpetra::Experimental::GEMM ("N", "N", STS::one (), dmat, currentVal,
757  STS::zero (), matTmp);
758 #endif
759  //blockMatOpts.square_matrix_matrix_multiply(reinterpret_cast<impl_scalar_type*>(dmat.data ()), reinterpret_cast<impl_scalar_type*>(currentVal.data ()), reinterpret_cast<impl_scalar_type*>(matTmp.data ()), blockSize_);
760  //currentVal.assign(matTmp);
761  Tpetra::Experimental::COPY (matTmp, currentVal);
762  }
763 
764  if (NumU) {
765  // Replace current row of L and U
766  U_block_->replaceLocalValues (local_row, &InI[NumL+1], &InV[blockMatSize*(NumL+1)], NumU);
767  }
768 
769 #ifndef IFPACK2_RBILUK_INITIAL
770  // Reset column flags
771  for (size_t j = 0; j < NumIn; ++j) {
772  colflag[InI[j]] = -1;
773  }
774 #else
775  // Reset column flags
776  for (size_t j = 0; j < num_cols; ++j) {
777  colflag[j] = -1;
778  }
779 #endif
780  }
781 
782  } // Stop timing
783 
784  // Sync everything back to device, for efficient solves.
785  {
786  typedef typename block_crs_matrix_type::device_type device_type;
787  if (! A_block_.is_null ()) {
789  Teuchos::rcp_const_cast<block_crs_matrix_type> (A_block_);
790  A_nc->template sync<device_type> ();
791  }
792  L_block_->template sync<device_type> ();
793  U_block_->template sync<device_type> ();
794  D_block_->template sync<device_type> ();
795  }
796 
797  this->isComputed_ = true;
798  this->numCompute_ += 1;
799  this->computeTime_ += timer.totalElapsedTime ();
800 }
801 
802 
803 template<class MatrixType>
804 void
806 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
807  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
808  Teuchos::ETransp mode,
809  scalar_type alpha,
810  scalar_type beta) const
811 {
812  using Teuchos::RCP;
813  typedef Tpetra::Experimental::BlockMultiVector<scalar_type,
815  typedef Tpetra::MultiVector<scalar_type,
816  local_ordinal_type, global_ordinal_type, node_type> MV;
817 
819  A_block_.is_null (), std::runtime_error, "Ifpack2::Experimental::RBILUK::apply: The matrix is "
820  "null. Please call setMatrix() with a nonnull input, then initialize() "
821  "and compute(), before calling this method.");
823  ! this->isComputed (), std::runtime_error,
824  "Ifpack2::Experimental::RBILUK::apply: If you have not yet called compute(), "
825  "you must call compute() before calling this method.");
826  TEUCHOS_TEST_FOR_EXCEPTION(
827  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
828  "Ifpack2::Experimental::RBILUK::apply: X and Y do not have the same number of columns. "
829  "X.getNumVectors() = " << X.getNumVectors ()
830  << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
831  TEUCHOS_TEST_FOR_EXCEPTION(
832  STS::isComplex && mode == Teuchos::CONJ_TRANS, std::logic_error,
833  "Ifpack2::Experimental::RBILUK::apply: mode = Teuchos::CONJ_TRANS is not implemented for "
834  "complex Scalar type. Please talk to the Ifpack2 developers to get this "
835  "fixed. There is a FIXME in this file about this very issue.");
836 
837  const local_ordinal_type blockMatSize = blockSize_*blockSize_;
838 
839  const local_ordinal_type rowStride = blockSize_;
840 
841  BMV yBlock (Y, * (A_block_->getGraph ()->getDomainMap ()), blockSize_);
842  const BMV xBlock (X, * (A_block_->getColMap ()), blockSize_);
843 
844  Teuchos::Array<scalar_type> lclarray(blockSize_);
845  little_vec_type lclvec((typename little_vec_type::value_type*)&lclarray[0], blockSize_);
846  const scalar_type one = STM::one ();
847  const scalar_type zero = STM::zero ();
848 
849  Teuchos::Time timer ("RBILUK::apply");
850  { // Start timing
851  Teuchos::TimeMonitor timeMon (timer);
852  if (alpha == one && beta == zero) {
853  if (mode == Teuchos::NO_TRANS) { // Solve L (D (U Y)) = X for Y.
854  // Start by solving L C = X for C. C must have the same Map
855  // as D. We have to use a temp multivector, since our
856  // implementation of triangular solves does not allow its
857  // input and output to alias one another.
858  //
859  // FIXME (mfh 24 Jan 2014) Cache this temp multivector.
860  const local_ordinal_type numVectors = xBlock.getNumVectors();
861  BMV cBlock (* (A_block_->getGraph ()->getDomainMap ()), blockSize_, numVectors);
862  BMV rBlock (* (A_block_->getGraph ()->getDomainMap ()), blockSize_, numVectors);
863  for (local_ordinal_type imv = 0; imv < numVectors; ++imv)
864  {
865  for (size_t i = 0; i < D_block_->getNodeNumRows(); ++i)
866  {
867  local_ordinal_type local_row = i;
868  little_vec_type xval = xBlock.getLocalBlock(local_row,imv);
869  little_vec_type cval = cBlock.getLocalBlock(local_row,imv);
870  //cval.assign(xval);
871  Tpetra::Experimental::COPY (xval, cval);
872 
873  local_ordinal_type NumL;
874  const local_ordinal_type * colValsL;
875  scalar_type * valsL;
876 
877  L_block_->getLocalRowView(local_row, colValsL, valsL, NumL);
878 
879  for (local_ordinal_type j = 0; j < NumL; ++j)
880  {
881  local_ordinal_type col = colValsL[j];
882  little_vec_type prevVal = cBlock.getLocalBlock(col, imv);
883 
884  const local_ordinal_type matOffset = blockMatSize*j;
885  little_block_type lij((typename little_block_type::value_type*) &valsL[matOffset],blockSize_,rowStride);
886 
887  //cval.matvecUpdate(-one, lij, prevVal);
888  Tpetra::Experimental::GEMV (-one, lij, prevVal, cval);
889  }
890  }
891  }
892 
893  // Solve D R = C. Note that D has been replaced by D^{-1} at this point.
894  D_block_->applyBlock(cBlock, rBlock);
895 
896  // Solve U Y = R.
897  for (local_ordinal_type imv = 0; imv < numVectors; ++imv)
898  {
899  const local_ordinal_type numRows = D_block_->getNodeNumRows();
900  for (local_ordinal_type i = 0; i < numRows; ++i)
901  {
902  local_ordinal_type local_row = (numRows-1)-i;
903  little_vec_type rval = rBlock.getLocalBlock(local_row,imv);
904  little_vec_type yval = yBlock.getLocalBlock(local_row,imv);
905  //yval.assign(rval);
906  Tpetra::Experimental::COPY (rval, yval);
907 
908  local_ordinal_type NumU;
909  const local_ordinal_type * colValsU;
910  scalar_type * valsU;
911 
912  U_block_->getLocalRowView(local_row, colValsU, valsU, NumU);
913 
914  for (local_ordinal_type j = 0; j < NumU; ++j)
915  {
916  local_ordinal_type col = colValsU[NumU-1-j];
917  little_vec_type prevVal = yBlock.getLocalBlock(col, imv);
918 
919  const local_ordinal_type matOffset = blockMatSize*(NumU-1-j);
920  little_block_type uij((typename little_block_type::value_type*) &valsU[matOffset], blockSize_, rowStride);
921 
922  //yval.matvecUpdate(-one, uij, prevVal);
923  Tpetra::Experimental::GEMV (-one, uij, prevVal, yval);
924  }
925  }
926  }
927  }
928  else { // Solve U^P (D^P (L^P Y)) = X for Y (where P is * or T).
929  TEUCHOS_TEST_FOR_EXCEPTION(
930  true, std::runtime_error,
931  "Ifpack2::Experimental::RBILUK::apply: transpose apply is not implemented for the block algorithm. ");
932  }
933  }
934  else { // alpha != 1 or beta != 0
935  if (alpha == zero) {
936  if (beta == zero) {
937  Y.putScalar (zero);
938  } else {
939  Y.scale (beta);
940  }
941  } else { // alpha != zero
942  MV Y_tmp (Y.getMap (), Y.getNumVectors ());
943  apply (X, Y_tmp, mode);
944  Y.update (alpha, Y_tmp, beta);
945  }
946  }
947  } // Stop timing
948 
949  this->numApply_ += 1;
950  this->applyTime_ = timer.totalElapsedTime ();
951 }
952 
953 
954 template<class MatrixType>
956 {
957  std::ostringstream os;
958 
959  // Output is a valid YAML dictionary in flow style. If you don't
960  // like everything on a single line, you should call describe()
961  // instead.
962  os << "\"Ifpack2::Experimental::RBILUK\": {";
963  os << "Initialized: " << (this->isInitialized () ? "true" : "false") << ", "
964  << "Computed: " << (this->isComputed () ? "true" : "false") << ", ";
965 
966  os << "Level-of-fill: " << this->getLevelOfFill() << ", ";
967 
968  if (A_block_.is_null ()) {
969  os << "Matrix: null";
970  }
971  else {
972  os << "Global matrix dimensions: ["
973  << A_block_->getGlobalNumRows () << ", " << A_block_->getGlobalNumCols () << "]"
974  << ", Global nnz: " << A_block_->getGlobalNumEntries();
975  }
976 
977  os << "}";
978  return os.str ();
979 }
980 
981 } // namespace Experimental
982 
983 } // namespace Ifpack2
984 
985 // FIXME (mfh 26 Aug 2015) We only need to do instantiation for
986 // MatrixType = Tpetra::RowMatrix. Conversions to BlockCrsMatrix are
987 // handled internally via dynamic cast.
988 
989 #define IFPACK2_EXPERIMENTAL_RBILUK_INSTANT(S,LO,GO,N) \
990  template class Ifpack2::Experimental::RBILUK< Tpetra::BlockCrsMatrix<S, LO, GO, N> >; \
991  template class Ifpack2::Experimental::RBILUK< Tpetra::RowMatrix<S, LO, GO, N> >;
992 
993 #endif
void initialize()
Initialize by computing the symbolic incomplete factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:187
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:145
virtual ~RBILUK()
Destructor (declared virtual for memory safety).
Definition: Ifpack2_Experimental_RBILUK_def.hpp:81
Teuchos::RCP< const block_crs_matrix_type > getBlockMatrix() const
Get the input matrix.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:182
const block_crs_matrix_type & getUBlock() const
Return the U factor of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:139
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:136
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_Experimental_RBILUK_def.hpp:806
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:151
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:243
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:148
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:142
T * getRawPtr() const
const block_crs_matrix_type & getDBlock() const
Return the diagonal entries of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:125
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void compute()
Compute the (numeric) incomplete factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:458
File for utility functions.
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:97
void setMatrix(const Teuchos::RCP< const block_crs_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:86
const block_crs_matrix_type & getLBlock() const
Return the L factor of the ILU factorization.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:111
Sparse matrix (Tpetra::RowMatrix subclass) with ghost rows.
Definition: Ifpack2_OverlappingRowMatrix_decl.hpp:58
double totalElapsedTime(bool readCurrentTime=false) const
Access only local rows and columns of a sparse matrix.
Definition: Ifpack2_LocalFilter_decl.hpp:160
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:157
std::string description() const
A one-line description of this object.
Definition: Ifpack2_Experimental_RBILUK_def.hpp:955
ILU(k) factorization of a given Tpetra::BlockCrsMatrix.
Definition: Ifpack2_Experimental_RBILUK_decl.hpp:128