Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_SparseContainer_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_SPARSECONTAINER_DEF_HPP
44 #define IFPACK2_SPARSECONTAINER_DEF_HPP
45 
47 #ifdef HAVE_MPI
48 #include <mpi.h>
50 #else
51 #include "Teuchos_DefaultSerialComm.hpp"
52 #endif
54 
55 namespace Ifpack2 {
56 
57 //==============================================================================
58 template<class MatrixType, class InverseType>
62  const Teuchos::RCP<const import_type>& importer,
63  int OverlapLevel,
64  scalar_type DampingFactor) :
65  Container<MatrixType> (matrix, partitions, importer, OverlapLevel,
66  DampingFactor),
67  IsInitialized_ (false),
68  IsComputed_ (false),
69 #ifdef HAVE_MPI
70  localComm_ (Teuchos::rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF)))
71 #else
72  localComm_ (Teuchos::rcp (new Teuchos::SerialComm<int> ()))
73 #endif // HAVE_MPI
74 {
75  global_ordinal_type indexBase = 0;
76  for(int i = 0; i < this->numBlocks_; i++)
77  localMaps_.emplace_back(this->blockRows_[i], indexBase, localComm_);
78 }
79 
80 //==============================================================================
81 template<class MatrixType, class InverseType>
84  const Teuchos::Array<local_ordinal_type>& localRows) :
85  Container<MatrixType> (matrix, localRows),
86  IsInitialized_(false),
87  IsComputed_(false),
88 #ifdef HAVE_MPI
89  localComm_ (Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_SELF)))
90 #else
91  localComm_ (Teuchos::rcp(new Teuchos::SerialComm<int>()))
92 #endif // HAVE_MPI
93 {
94  global_ordinal_type indexBase = 0;
95  localMaps_.emplace_back(this->blockRows_[0], indexBase, localComm_);
96 }
97 
98 //==============================================================================
99 template<class MatrixType,class InverseType>
101 {
102  for(auto inv : Inverses_)
103  delete inv.get();
104 }
105 
106 //==============================================================================
107 template<class MatrixType, class InverseType>
109 {
110  return IsInitialized_;
111 }
112 
113 //==============================================================================
114 template<class MatrixType, class InverseType>
116 {
117  return IsComputed_;
118 }
119 
120 //==============================================================================
121 template<class MatrixType, class InverseType>
123 {
124  List_ = List;
125 }
126 
127 //==============================================================================
128 template<class MatrixType, class InverseType>
130 {
131  using Teuchos::RCP;
132  // We assume that if you called this method, you intend to recompute
133  // everything. Thus, we release references to all the internal
134  // objects. We do this first to save memory. (In an RCP
135  // assignment, the right-hand side and left-hand side coexist before
136  // the left-hand side's reference count gets updated.)
137  IsInitialized_ = false;
138  IsComputed_ = false;
139 
140  // (Re)create the CrsMatrices that will contain the
141  // local matrices to use for solves.
142  diagBlocks_.reserve(this->numBlocks_);
143  Inverses_.reserve(this->numBlocks_);
144  for(int i = 0; i < this->numBlocks_; i++)
145  {
146  // Create a local map for the block, with same size as block has rows.
147  // Note: this map isn't needed elsewhere in SparseContainer, but the
148  // diagBlocks_[...] will keep it alive
149  RCP<InverseMap> tempMap(new InverseMap(this->blockRows_[i], 0, localComm_));
150  diagBlocks_.emplace_back(new InverseCrs(tempMap, 0));
151  // Create the inverse operator using the local matrix. We give it
152  // the matrix here, but don't call its initialize() or compute()
153  // methods yet, since we haven't initialized the matrix yet.
154  Inverses_.push_back(ptr(new InverseType(diagBlocks_[i])));
155  }
156  IsInitialized_ = true;
157 }
158 
159 //==============================================================================
160 template<class MatrixType, class InverseType>
162 {
163  IsComputed_ = false;
164  if (! this->isInitialized ()) {
165  this->initialize ();
166  }
167 
168  // Extract the submatrix.
169  this->extract ();
170 
171  // The inverse operator already has a pointer to the submatrix. Now
172  // that the submatrix is filled in, we can initialize and compute
173  // the inverse operator.
174  for(int i = 0; i < this->numBlocks_; i++)
175  Inverses_[i]->initialize ();
176  for(int i = 0; i < this->numBlocks_; i++)
177  Inverses_[i]->compute ();
178  IsComputed_ = true;
179 }
180 
181 //==============================================================================
182 template<class MatrixType, class InverseType>
184 {
185  for(auto inv : Inverses_)
186  delete inv.get();
187  Inverses_.clear();
188  diagBlocks_.clear();
190 }
191 
192 //==============================================================================
193 template<class MatrixType, class InverseType>
195 applyImpl (inverse_mv_type& X,
196  inverse_mv_type& Y,
197  int blockIndex,
198  int /* stride */,
199  Teuchos::ETransp mode,
200  InverseScalar alpha,
201  InverseScalar beta) const
202 {
204  Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() != X.getLocalLength(),
205  std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
206  "operator and X have incompatible dimensions (" <<
207  Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() << " resp. "
208  << X.getLocalLength() << "). Please report this bug to "
209  "the Ifpack2 developers.");
211  Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() != Y.getLocalLength(),
212  std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
213  "operator and Y have incompatible dimensions (" <<
214  Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() << " resp. "
215  << Y.getLocalLength() << "). Please report this bug to "
216  "the Ifpack2 developers.");
217  Inverses_[blockIndex]->apply(X, Y, mode, alpha, beta);
218 }
219 
220 template<class MatrixType, class InverseType>
222 apply (HostView& X,
223  HostView& Y,
224  int blockIndex,
225  int stride,
226  Teuchos::ETransp mode,
227  scalar_type alpha,
228  scalar_type beta) const
229 {
230  using Teuchos::ArrayView;
231  using Teuchos::as;
232  using Teuchos::RCP;
233  using Teuchos::rcp;
234 
235  // The InverseType template parameter might have different template
236  // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
237  // example, MatrixType (a global object) might use a bigger GO
238  // (global ordinal) type than InverseType (which deals with the
239  // diagonal block, a local object). This means that we might have
240  // to convert X and Y to the Tpetra::MultiVector specialization that
241  // InverseType wants. This class' X_ and Y_ internal fields are of
242  // the right type for InverseType, so we can use those as targets.
243 
244  // Tpetra::MultiVector specialization corresponding to InverseType.
246  size_t numVecs = X.extent(1);
247 
249  ! IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply: "
250  "You must have called the compute() method before you may call apply(). "
251  "You may call the apply() method as many times as you want after calling "
252  "compute() once, but you must have called compute() at least once.");
254  X.extent(1) != Y.extent(1), std::runtime_error,
255  "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
256  "vectors. X has " << X.extent(1)
257  << ", but Y has " << Y.extent(1) << ".");
258 
259  if (numVecs == 0) {
260  return; // done! nothing to do
261  }
262 
263  const local_ordinal_type numRows_ = this->blockRows_[blockIndex];
264 
265  // The operator Inverse_ works on a permuted subset of the local
266  // parts of X and Y. The subset and permutation are defined by the
267  // index array returned by getLocalRows(). If the permutation is
268  // trivial and the subset is exactly equal to the local indices,
269  // then we could use the local parts of X and Y exactly, without
270  // needing to permute. Otherwise, we have to use temporary storage
271  // to permute X and Y. For now, we always use temporary storage.
272  //
273  // Create temporary permuted versions of the input and output.
274  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
275  // store the permuted versions of X resp. Y. Note that X_local has
276  // the domain Map of the operator, which may be a permuted subset of
277  // the local Map corresponding to X.getMap(). Similarly, Y_local
278  // has the range Map of the operator, which may be a permuted subset
279  // of the local Map corresponding to Y.getMap(). numRows_ here
280  // gives the number of rows in the row Map of the local Inverse_
281  // operator.
282  //
283  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
284  // here that the row Map and the range Map of that operator are
285  // the same.
286  //
287  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
288  // really belongs in Tpetra.
289 
290  if(invX.size() == 0)
291  {
292  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
293  invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
294  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
295  invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
296  }
297  inverse_mv_type& X_local = invX[blockIndex];
299  X_local.getLocalLength() != (size_t) numRows_, std::logic_error,
300  "Ifpack2::SparseContainer::apply: "
301  "X_local has length " << X_local.getLocalLength() << ", which does "
302  "not match numRows_ = " << numRows_ << ". Please report this bug to "
303  "the Ifpack2 developers.");
304  const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
305  mvgs.gatherMVtoView(X_local, X, localRows);
306 
307  // We must gather the output multivector Y even on input to
308  // Inverse_->apply(), since the Inverse_ operator might use it as an
309  // initial guess for a linear solve. We have no way of knowing
310  // whether it does or does not.
311 
312  inverse_mv_type& Y_local = invY[blockIndex];
314  Y_local.getLocalLength () != (size_t) numRows_, std::logic_error,
315  "Ifpack2::SparseContainer::apply: "
316  "Y_local has length " << Y_local.getLocalLength () << ", which does "
317  "not match numRows_ = " << numRows_ << ". Please report this bug to "
318  "the Ifpack2 developers.");
319  mvgs.gatherMVtoView(Y_local, Y, localRows);
320 
321  // Apply the local operator:
322  // Y_local := beta*Y_local + alpha*M^{-1}*X_local
323  this->applyImpl(X_local, Y_local, blockIndex, stride, mode, as<InverseScalar>(alpha),
324  as<InverseScalar>(beta));
325 
326  // Scatter the permuted subset output vector Y_local back into the
327  // original output multivector Y.
328  mvgs.scatterMVtoView(Y, Y_local, localRows);
329 }
330 
331 //==============================================================================
332 template<class MatrixType, class InverseType>
334 weightedApply (HostView& X,
335  HostView& Y,
336  HostView& D,
337  int blockIndex,
338  int /* stride */,
339  Teuchos::ETransp mode,
340  scalar_type alpha,
341  scalar_type beta) const
342 {
343  using Teuchos::ArrayRCP;
344  using Teuchos::ArrayView;
345  using Teuchos::Range1D;
346  using Teuchos::Ptr;
347  using Teuchos::ptr;
348  using Teuchos::RCP;
349  using Teuchos::rcp;
350  using Teuchos::rcp_const_cast;
351  using std::cerr;
352  using std::endl;
354 
355  // The InverseType template parameter might have different template
356  // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
357  // example, MatrixType (a global object) might use a bigger GO
358  // (global ordinal) type than InverseType (which deals with the
359  // diagonal block, a local object). This means that we might have
360  // to convert X and Y to the Tpetra::MultiVector specialization that
361  // InverseType wants. This class' X_ and Y_ internal fields are of
362  // the right type for InverseType, so we can use those as targets.
363 
364  // Tpetra::Vector specialization corresponding to InverseType.
365  typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal, InverseGlobalOrdinal, InverseNode> inverse_vector_type;
366 
368  const size_t numVecs = X.extent(1);
369 
371  ! IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::"
372  "weightedApply: You must have called the compute() method before you may "
373  "call apply(). You may call the apply() method as many times as you want "
374  "after calling compute() once, but you must have called compute() at least "
375  "once.");
377  X.extent(1) != Y.extent(1), std::runtime_error,
378  "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
379  "of vectors. X has " << X.extent(1) << ", but Y has "
380  << Y.extent(1) << ".");
381 
382  if (numVecs == 0) {
383  return; // done! nothing to do
384  }
385 
386  // The operator Inverse_ works on a permuted subset of the local
387  // parts of X and Y. The subset and permutation are defined by the
388  // index array returned by getLocalRows(). If the permutation is
389  // trivial and the subset is exactly equal to the local indices,
390  // then we could use the local parts of X and Y exactly, without
391  // needing to permute. Otherwise, we have to use temporary storage
392  // to permute X and Y. For now, we always use temporary storage.
393  //
394  // Create temporary permuted versions of the input and output.
395  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
396  // store the permuted versions of X resp. Y. Note that X_local has
397  // the domain Map of the operator, which may be a permuted subset of
398  // the local Map corresponding to X.getMap(). Similarly, Y_local
399  // has the range Map of the operator, which may be a permuted subset
400  // of the local Map corresponding to Y.getMap(). numRows_ here
401  // gives the number of rows in the row Map of the local Inverse_
402  // operator.
403  //
404  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
405  // here that the row Map and the range Map of that operator are
406  // the same.
407  //
408  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
409  // really belongs in Tpetra.
410  const local_ordinal_type numRows = this->blockRows_[blockIndex];
411 
412  if(invX.size() == 0)
413  {
414  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
415  invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
416  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
417  invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
418  }
419  inverse_mv_type& X_local = invX[blockIndex];
420  const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(blockIndex);
421  mvgs.gatherMVtoView(X_local, X, localRows);
422 
423  // We must gather the output multivector Y even on input to
424  // Inverse_->apply(), since the Inverse_ operator might use it as an
425  // initial guess for a linear solve. We have no way of knowing
426  // whether it does or does not.
427 
428  inverse_mv_type Y_local = invY[blockIndex];
430  Y_local.getLocalLength() != (size_t) numRows, std::logic_error,
431  "Ifpack2::SparseContainer::weightedApply: "
432  "Y_local has length " << X_local.getLocalLength() << ", which does "
433  "not match numRows_ = " << numRows << ". Please report this bug to "
434  "the Ifpack2 developers.");
435  mvgs.gatherMVtoView(Y_local, Y, localRows);
436 
437  // Apply the diagonal scaling D to the input X. It's our choice
438  // whether the result has the original input Map of X, or the
439  // permuted subset Map of X_local. If the latter, we also need to
440  // gather D into the permuted subset Map. We choose the latter, to
441  // save memory and computation. Thus, we do the following:
442  //
443  // 1. Gather D into a temporary vector D_local.
444  // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
445  // 3. Compute X_scaled := diag(D_loca) * X_local.
446 
447  inverse_vector_type D_local(Inverses_[blockIndex]->getDomainMap());
449  D_local.getLocalLength() != (size_t) this->blockRows_[blockIndex], std::logic_error,
450  "Ifpack2::SparseContainer::weightedApply: "
451  "D_local has length " << X_local.getLocalLength () << ", which does "
452  "not match numRows_ = " << this->blockRows_[blockIndex] << ". Please report this bug to "
453  "the Ifpack2 developers.");
454  mvgs.gatherMVtoView(D_local, D, localRows);
455  inverse_mv_type X_scaled(Inverses_[blockIndex]->getDomainMap(), numVecs);
456  X_scaled.elementWiseMultiply(STS::one(), D_local, X_local, STS::zero());
457 
458  // Y_temp will hold the result of M^{-1}*X_scaled. If beta == 0, we
459  // can write the result of Inverse_->apply() directly to Y_local, so
460  // Y_temp may alias Y_local. Otherwise, if beta != 0, we need
461  // temporary storage for M^{-1}*X_scaled, so Y_temp must be
462  // different than Y_local.
463  Ptr<inverse_mv_type> Y_temp;
464  bool deleteYT = false;
465  if (beta == STS::zero ()) {
466  Y_temp = ptr(&Y_local);
467  } else {
468  Y_temp = ptr(new inverse_mv_type(Inverses_[blockIndex]->getRangeMap(), numVecs));
469  deleteYT = true;
470  }
471  // Apply the local operator: Y_temp := M^{-1} * X_scaled
472  Inverses_[blockIndex]->apply(X_scaled, *Y_temp, mode);
473  // Y_local := beta * Y_local + alpha * diag(D_local) * Y_tmp.
474  //
475  // Note that we still use the permuted subset scaling D_local here,
476  // because Y_temp has the same permuted subset Map. That's good, in
477  // fact, because it's a subset: less data to read and multiply.
478  Y_local.elementWiseMultiply(alpha, D_local, *Y_temp, beta);
479  if(deleteYT)
480  delete Y_temp.get();
481  // Copy the permuted subset output vector Y_local into the original
482  // output multivector Y.
483  mvgs.scatterMVtoView(Y, Y_local, localRows);
484 }
485 
486 //==============================================================================
487 template<class MatrixType, class InverseType>
488 std::ostream& SparseContainer<MatrixType,InverseType>::print(std::ostream& os) const
489 {
490  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
491  fos.setOutputToRootOnly(0);
492  describe(fos);
493  return(os);
494 }
495 
496 //==============================================================================
497 template<class MatrixType, class InverseType>
499 {
500  std::ostringstream oss;
501  oss << "\"Ifpack2::SparseContainer\": {";
502  if (isInitialized()) {
503  if (isComputed()) {
504  oss << "status = initialized, computed";
505  }
506  else {
507  oss << "status = initialized, not computed";
508  }
509  }
510  else {
511  oss << "status = not initialized, not computed";
512  }
513  for(int i = 0; i < this->numBlocks_; i++)
514  {
515  oss << ", Block Inverse " << i << ": {";
516  oss << Inverses_[i]->description();
517  oss << "}";
518  }
519  oss << "}";
520  return oss.str();
521 }
522 
523 //==============================================================================
524 template<class MatrixType, class InverseType>
526 {
527  using std::endl;
528  if(verbLevel==Teuchos::VERB_NONE) return;
529  os << "================================================================================" << endl;
530  os << "Ifpack2::SparseContainer" << endl;
531  for(int i = 0; i < this->numBlocks_; i++)
532  {
533  os << "Block " << i << " rows: = " << this->blockRows_[i] << endl;
534  }
535  os << "isInitialized() = " << IsInitialized_ << endl;
536  os << "isComputed() = " << IsComputed_ << endl;
537  os << "================================================================================" << endl;
538  os << endl;
539 }
540 
541 //==============================================================================
542 template<class MatrixType, class InverseType>
544 extract ()
545 {
546  using Teuchos::Array;
547  using Teuchos::ArrayView;
548 
549  auto& A = *this->inputMatrix_;
550  const size_t MatrixInNumRows = A.getNodeNumRows ();
551 
552  // Sanity checking
553  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
554  {
555  const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
556  for (local_ordinal_type j = 0; j < localRows.size(); j++)
557  {
559  localRows[j] < 0 ||
560  static_cast<size_t> (localRows[j]) >= MatrixInNumRows,
561  std::runtime_error, "Ifpack2::SparseContainer::extract: localRows[j="
562  << j << "] = " << localRows[j] << ", which is out of the valid range. "
563  "This probably means that compute() has not yet been called.");
564  }
565  }
566 
567  const size_t maxNumEntriesInRow = A.getNodeMaxNumRowEntries();
568  Array<scalar_type> Values;
569  Array<local_ordinal_type> Indices;
570  Array<InverseScalar> Values_insert;
571  Array<InverseGlobalOrdinal> Indices_insert;
572 
573  Values.resize (maxNumEntriesInRow);
574  Indices.resize (maxNumEntriesInRow);
575  Values_insert.resize (maxNumEntriesInRow);
576  Indices_insert.resize (maxNumEntriesInRow);
577 
578  const InverseLocalOrdinal INVALID = Teuchos::OrdinalTraits<InverseLocalOrdinal>::invalid ();
579  for(local_ordinal_type i = 0; i < this->numBlocks_; i++)
580  {
581  const local_ordinal_type numRows_ = this->blockRows_[i];
582  const ArrayView<const local_ordinal_type> localRows = this->getLocalRows(i);
583  for (local_ordinal_type j = 0; j < numRows_; j++)
584  {
585  const local_ordinal_type localRow = this->partitions_[this->partitionIndices_[i] + j];
586  size_t numEntries;
587  A.getLocalRowCopy(localRow, Indices(), Values(), numEntries);
588  size_t num_entries_found = 0;
589  for(size_t k = 0; k < numEntries; ++k)
590  {
591  const local_ordinal_type localCol = Indices[k];
592  // Skip off-process elements
593  //
594  // FIXME (mfh 24 Aug 2013) This assumes the following:
595  //
596  // 1. The column and row Maps begin with the same set of
597  // on-process entries, in the same order. That is,
598  // on-process row and column indices are the same.
599  // 2. All off-process indices in the column Map of the input
600  // matrix occur after that initial set.
601  if(static_cast<size_t> (localCol) >= MatrixInNumRows)
602  continue;
603  // for local column IDs, look for each ID in the list
604  // of columns hosted by this object
605  InverseLocalOrdinal jj = INVALID;
606  for(local_ordinal_type kk = 0; kk < numRows_; kk++)
607  {
608  if(localRows[kk] == localCol)
609  jj = kk;
610  }
611  if (jj != INVALID)
612  {
613  Indices_insert[num_entries_found] = localMaps_[i].getGlobalElement(jj);
614  Values_insert[num_entries_found] = Values[k];
615  num_entries_found++;
616  }
617  }
618  diagBlocks_[i]->insertGlobalValues(j, Indices_insert (0, num_entries_found),
619  Values_insert (0, num_entries_found));
620  }
621  // FIXME (mfh 24 Aug 2013) If we generalize the current set of
622  // assumptions on the column and row Maps (see note above), we may
623  // need to supply custom domain and range Maps to fillComplete().
624  diagBlocks_[i]->fillComplete ();
625  }
626 }
627 
628 template<typename MatrixType, typename InverseType>
630 {
632 #ifdef HAVE_IFPACK2_AMESOS2
634  if(std::is_same<InverseType, ILUTInverse>::value)
635  {
636  return "SparseILUT";
637  }
638  else if(std::is_same<InverseType, AmesosInverse>::value)
639  {
640  return "SparseAmesos";
641  }
642  else
643  {
644  throw std::logic_error("InverseType for SparseContainer must be Ifpack2::ILUT or Details::Amesos2Wrapper");
645  }
646 #else
647  // Macros can't have commas in their arguments, so we have to
648  // compute the bool first argument separately.
649  constexpr bool inverseTypeIsILUT = std::is_same<InverseType, ILUTInverse>::value;
650  TEUCHOS_TEST_FOR_EXCEPTION(! inverseTypeIsILUT, std::logic_error,
651  "InverseType for SparseContainer must be Ifpack2::ILUT<ROW>");
652  return "SparseILUT"; //the only supported sparse container specialization if no Amesos2
653 #endif
654 }
655 
656 } // namespace Ifpack2
657 
658 // For ETI
659 #include "Ifpack2_ILUT.hpp"
660 #ifdef HAVE_IFPACK2_AMESOS2
661 #include "Ifpack2_Details_Amesos2Wrapper.hpp"
662 #endif
663 
664 // There's no need to instantiate for CrsMatrix too. All Ifpack2
665 // preconditioners can and should do dynamic casts if they need a type
666 // more specific than RowMatrix.
667 
668 #ifdef HAVE_IFPACK2_AMESOS2
669 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
670  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
671  Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; \
672  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
673  Ifpack2::Details::Amesos2Wrapper<Tpetra::RowMatrix<S,LO,GO,N> > >;
674 #else
675 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
676  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S,LO,GO,N>, \
677  Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> > >;
678 #endif
679 #endif // IFPACK2_SPARSECONTAINER_DEF_HPP
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_SparseContainer_def.hpp:629
virtual bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_SparseContainer_def.hpp:115
Store and solve a local sparse linear problem.
Definition: Ifpack2_SparseContainer_decl.hpp:134
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container.hpp:421
Teuchos::Array< local_ordinal_type > blockRows_
Number of rows in each block.
Definition: Ifpack2_Container.hpp:425
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to the given FancyOStream.
Definition: Ifpack2_SparseContainer_def.hpp:525
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void initialize()
Do all set-up operations that only require matrix structure.
Definition: Ifpack2_SparseContainer_def.hpp:129
ILUT (incomplete LU factorization with threshold) of a Tpetra sparse matrix.
Definition: Ifpack2_ILUT_decl.hpp:91
virtual std::ostream & print(std::ostream &os) const
Print information about this object to the given output stream.
Definition: Ifpack2_SparseContainer_def.hpp:488
virtual void weightedApply(HostView &X, HostView &Y, HostView &W, int blockIndex, int stride, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:334
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_SparseContainer_def.hpp:498
virtual bool isInitialized() const
Whether the container has been successfully initialized.
Definition: Ifpack2_SparseContainer_def.hpp:108
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Wrapper class for direct solvers in Amesos2.
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:102
SparseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions, const Teuchos::RCP< const import_type > &importer, int OverlapLevel, scalar_type DampingFactor)
Constructor.
Definition: Ifpack2_SparseContainer_def.hpp:60
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual void apply(HostView &X, HostView &Y, int blockIndex, int stride, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Compute Y := alpha * M^{-1} X + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:222
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_SparseContainer_def.hpp:100
TypeTo as(const TypeFrom &t)
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition: Ifpack2_SparseContainer_def.hpp:122
virtual void compute()
Initialize and compute all blocks.
Definition: Ifpack2_SparseContainer_def.hpp:161
Ifpack2::SparseContainer class declaration.
void clearBlocks()
Definition: Ifpack2_SparseContainer_def.hpp:183
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85