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