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>
61  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
63  bool pointIndexed) :
64  ContainerImpl<MatrixType, InverseScalar> (matrix, partitions, pointIndexed),
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>
76 
77 //==============================================================================
78 template<class MatrixType, class InverseType>
80 {
81  List_ = List;
82 }
83 
84 //==============================================================================
85 template<class MatrixType, class InverseType>
87 {
88  // We assume that if you called this method, you intend to recompute
89  // everything. Thus, we release references to all the internal
90  // objects. We do this first to save memory. (In an RCP
91  // assignment, the right-hand side and left-hand side coexist before
92  // the left-hand side's reference count gets updated.)
93 
94  //Will create the diagonal blocks and their inverses
95  //in extract()
96  diagBlocks_.assign(this->numBlocks_, Teuchos::null);
97  Inverses_.assign(this->numBlocks_, Teuchos::null);
98 
99  this->IsInitialized_ = true;
100  this->IsComputed_ = false;
101 }
102 
103 //==============================================================================
104 template<class MatrixType, class InverseType>
106 {
107  this->IsComputed_ = false;
108  if (!this->isInitialized ()) {
109  this->initialize ();
110  }
111 
112  // Extract the submatrices.
113  this->extract();
114 
115  // The inverse operator already has a pointer to the submatrix. Now
116  // that the submatrix is filled in, we can initialize and compute
117  // the inverse operator.
118  for(int i = 0; i < this->numBlocks_; i++)
119  {
120  Inverses_[i]->setParameters(List_);
121  Inverses_[i]->initialize ();
122  }
123  for(int i = 0; i < this->numBlocks_; i++)
124  Inverses_[i]->compute ();
125  this->IsComputed_ = true;
126 }
127 
128 //==============================================================================
129 template<class MatrixType, class InverseType>
131 {
132  for(auto inv : Inverses_)
133  delete inv.get();
134  Inverses_.clear();
135  diagBlocks_.clear();
137 }
138 
139 //==============================================================================
140 template<class MatrixType, class InverseType>
142 solveBlockMV(inverse_mv_type& X,
143  inverse_mv_type& Y,
144  int blockIndex,
145  Teuchos::ETransp mode,
146  InverseScalar alpha,
147  InverseScalar beta) const
148 {
150  Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() != X.getLocalLength(),
151  std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
152  "operator and X have incompatible dimensions (" <<
153  Inverses_[blockIndex]->getDomainMap()->getNodeNumElements() << " resp. "
154  << X.getLocalLength() << "). Please report this bug to "
155  "the Ifpack2 developers.");
157  Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() != Y.getLocalLength(),
158  std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
159  "operator and Y have incompatible dimensions (" <<
160  Inverses_[blockIndex]->getRangeMap()->getNodeNumElements() << " resp. "
161  << Y.getLocalLength() << "). Please report this bug to "
162  "the Ifpack2 developers.");
163  Inverses_[blockIndex]->apply(X, Y, mode, alpha, beta);
164 }
165 
166 template<class MatrixType, class InverseType>
169  HostView Y,
170  int blockIndex,
171  Teuchos::ETransp mode,
172  SC alpha,
173  SC beta) const
174 {
175  using Teuchos::ArrayView;
176  using Teuchos::as;
177 
178  // The InverseType template parameter might have different template
179  // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
180  // example, MatrixType (a global object) might use a bigger GO
181  // (global ordinal) type than InverseType (which deals with the
182  // diagonal block, a local object). This means that we might have
183  // to convert X and Y to the Tpetra::MultiVector specialization that
184  // InverseType wants. This class' X_ and Y_ internal fields are of
185  // the right type for InverseType, so we can use those as targets.
186 
187  // Tpetra::MultiVector specialization corresponding to InverseType.
189  size_t numVecs = X.extent(1);
190 
192  ! this->IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::apply: "
193  "You must have called the compute() method before you may call apply(). "
194  "You may call the apply() method as many times as you want after calling "
195  "compute() once, but you must have called compute() at least once.");
197  X.extent(1) != Y.extent(1), std::runtime_error,
198  "Ifpack2::SparseContainer::apply: X and Y have different numbers of "
199  "vectors. X has " << X.extent(1)
200  << ", but Y has " << Y.extent(1) << ".");
201 
202  if (numVecs == 0) {
203  return; // done! nothing to do
204  }
205 
206  const LO numRows = this->blockSizes_[blockIndex];
207 
208  // The operator Inverse_ works on a permuted subset of the local
209  // parts of X and Y. The subset and permutation are defined by the
210  // index array returned by getBlockRows(). If the permutation is
211  // trivial and the subset is exactly equal to the local indices,
212  // then we could use the local parts of X and Y exactly, without
213  // needing to permute. Otherwise, we have to use temporary storage
214  // to permute X and Y. For now, we always use temporary storage.
215  //
216  // Create temporary permuted versions of the input and output.
217  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
218  // store the permuted versions of X resp. Y. Note that X_local has
219  // the domain Map of the operator, which may be a permuted subset of
220  // the local Map corresponding to X.getMap(). Similarly, Y_local
221  // has the range Map of the operator, which may be a permuted subset
222  // of the local Map corresponding to Y.getMap(). numRows here
223  // gives the number of rows in the row Map of the local Inverse_
224  // operator.
225  //
226  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
227  // here that the row Map and the range Map of that operator are
228  // the same.
229  //
230  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
231  // really belongs in Tpetra.
232 
233  if(invX.size() == 0)
234  {
235  for(LO i = 0; i < this->numBlocks_; i++)
236  invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
237  for(LO i = 0; i < this->numBlocks_; i++)
238  invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
239  }
240  inverse_mv_type& X_local = invX[blockIndex];
242  X_local.getLocalLength() != size_t(numRows * this->scalarsPerRow_), std::logic_error,
243  "Ifpack2::SparseContainer::apply: "
244  "X_local has length " << X_local.getLocalLength() << ", which does "
245  "not match numRows = " << numRows * this->scalarsPerRow_ << ". Please report this bug to "
246  "the Ifpack2 developers.");
247  const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
248  if(this->scalarsPerRow_ == 1)
249  mvgs.gatherMVtoView(X_local, X, blockRows);
250  else
251  mvgs.gatherMVtoViewBlock(X_local, X, blockRows, this->scalarsPerRow_);
252 
253  // We must gather the output multivector Y even on input to
254  // Inverse_->apply(), since the Inverse_ operator might use it as an
255  // initial guess for a linear solve. We have no way of knowing
256  // whether it does or does not.
257 
258  inverse_mv_type& Y_local = invY[blockIndex];
260  Y_local.getLocalLength () != size_t(numRows * this->scalarsPerRow_), std::logic_error,
261  "Ifpack2::SparseContainer::apply: "
262  "Y_local has length " << Y_local.getLocalLength () << ", which does "
263  "not match numRows = " << numRows * this->scalarsPerRow_ << ". Please report this bug to "
264  "the Ifpack2 developers.");
265 
266  if(this->scalarsPerRow_ == 1)
267  mvgs.gatherMVtoView(Y_local, Y, blockRows);
268  else
269  mvgs.gatherMVtoViewBlock(Y_local, Y, blockRows, this->scalarsPerRow_);
270 
271  // Apply the local operator:
272  // Y_local := beta*Y_local + alpha*M^{-1}*X_local
273  this->solveBlockMV(X_local, Y_local, blockIndex, mode,
274  InverseScalar(alpha), InverseScalar(beta));
275 
276 
277  // Scatter the permuted subset output vector Y_local back into the
278  // original output multivector Y.
279  if(this->scalarsPerRow_ == 1)
280  mvgs.scatterMVtoView(Y, Y_local, blockRows);
281  else
282  mvgs.scatterMVtoViewBlock(Y, Y_local, blockRows, this->scalarsPerRow_);
283 }
284 
285 //==============================================================================
286 template<class MatrixType, class InverseType>
289  HostView Y,
290  HostView D,
291  int blockIndex,
292  Teuchos::ETransp mode,
293  SC alpha,
294  SC beta) const
295 {
296  using Teuchos::ArrayView;
297  using Teuchos::Range1D;
298  using std::cerr;
299  using std::endl;
301 
302  // The InverseType template parameter might have different template
303  // parameters (Scalar, LO, GO, and/or Node) than MatrixType. For
304  // example, MatrixType (a global object) might use a bigger GO
305  // (global ordinal) type than InverseType (which deals with the
306  // diagonal block, a local object). This means that we might have
307  // to convert X and Y to the Tpetra::MultiVector specialization that
308  // InverseType wants. This class' X_ and Y_ internal fields are of
309  // the right type for InverseType, so we can use those as targets.
310 
311  // Tpetra::Vector specialization corresponding to InverseType.
312  typedef Tpetra::Vector<InverseScalar, InverseLocalOrdinal, InverseGlobalOrdinal, InverseNode> inverse_vector_type;
313 
315  const size_t numVecs = X.extent(1);
316 
318  ! this->IsComputed_, std::runtime_error, "Ifpack2::SparseContainer::"
319  "weightedApply: You must have called the compute() method before you may "
320  "call apply(). You may call the apply() method as many times as you want "
321  "after calling compute() once, but you must have called compute() at least "
322  "once.");
324  X.extent(1) != Y.extent(1), std::runtime_error,
325  "Ifpack2::SparseContainer::weightedApply: X and Y have different numbers "
326  "of vectors. X has " << X.extent(1) << ", but Y has "
327  << Y.extent(1) << ".");
328 
329  //bmk 7-2019: BlockRelaxation already checked this, but if that changes...
331  this->scalarsPerRow_ > 1, std::logic_error, "Ifpack2::SparseContainer::weightedApply: "
332  "Use of block rows isn't allowed in overlapping Jacobi (the only method that uses weightedApply");
333 
334  if (numVecs == 0) {
335  return; // done! nothing to do
336  }
337 
338  // The operator Inverse_ works on a permuted subset of the local
339  // parts of X and Y. The subset and permutation are defined by the
340  // index array returned by getBlockRows(). If the permutation is
341  // trivial and the subset is exactly equal to the local indices,
342  // then we could use the local parts of X and Y exactly, without
343  // needing to permute. Otherwise, we have to use temporary storage
344  // to permute X and Y. For now, we always use temporary storage.
345  //
346  // Create temporary permuted versions of the input and output.
347  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
348  // store the permuted versions of X resp. Y. Note that X_local has
349  // the domain Map of the operator, which may be a permuted subset of
350  // the local Map corresponding to X.getMap(). Similarly, Y_local
351  // has the range Map of the operator, which may be a permuted subset
352  // of the local Map corresponding to Y.getMap(). numRows here
353  // gives the number of rows in the row Map of the local Inverse_
354  // operator.
355  //
356  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
357  // here that the row Map and the range Map of that operator are
358  // the same.
359  //
360  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
361  // really belongs in Tpetra.
362  const LO numRows = this->blockSizes_[blockIndex];
363 
364  if(invX.size() == 0)
365  {
366  for(LO i = 0; i < this->numBlocks_; i++)
367  invX.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
368  for(LO i = 0; i < this->numBlocks_; i++)
369  invY.emplace_back(Inverses_[i]->getDomainMap(), numVecs);
370  }
371  inverse_mv_type& X_local = invX[blockIndex];
372  const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
373  mvgs.gatherMVtoView(X_local, X, blockRows);
374 
375  // We must gather the output multivector Y even on input to
376  // Inverse_->apply(), since the Inverse_ operator might use it as an
377  // initial guess for a linear solve. We have no way of knowing
378  // whether it does or does not.
379 
380  inverse_mv_type Y_local = invY[blockIndex];
381  TEUCHOS_TEST_FOR_EXCEPTION(
382  Y_local.getLocalLength() != size_t(numRows), std::logic_error,
383  "Ifpack2::SparseContainer::weightedApply: "
384  "Y_local has length " << X_local.getLocalLength() << ", which does "
385  "not match numRows = " << numRows << ". Please report this bug to "
386  "the Ifpack2 developers.");
387  mvgs.gatherMVtoView(Y_local, Y, blockRows);
388 
389  // Apply the diagonal scaling D to the input X. It's our choice
390  // whether the result has the original input Map of X, or the
391  // permuted subset Map of X_local. If the latter, we also need to
392  // gather D into the permuted subset Map. We choose the latter, to
393  // save memory and computation. Thus, we do the following:
394  //
395  // 1. Gather D into a temporary vector D_local.
396  // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
397  // 3. Compute X_scaled := diag(D_loca) * X_local.
398 
399  inverse_vector_type D_local(Inverses_[blockIndex]->getDomainMap());
400  TEUCHOS_TEST_FOR_EXCEPTION(
401  D_local.getLocalLength() != size_t(this->blockSizes_[blockIndex]), std::logic_error,
402  "Ifpack2::SparseContainer::weightedApply: "
403  "D_local has length " << X_local.getLocalLength () << ", which does "
404  "not match numRows = " << this->blockSizes_[blockIndex] << ". Please report this bug to "
405  "the Ifpack2 developers.");
406  mvgs.gatherMVtoView(D_local, D, blockRows);
407  inverse_mv_type X_scaled(Inverses_[blockIndex]->getDomainMap(), numVecs);
408  X_scaled.elementWiseMultiply(STS::one(), D_local, X_local, STS::zero());
409 
410  // Y_temp will hold the result of M^{-1}*X_scaled. If beta == 0, we
411  // can write the result of Inverse_->apply() directly to Y_local, so
412  // Y_temp may alias Y_local. Otherwise, if beta != 0, we need
413  // temporary storage for M^{-1}*X_scaled, so Y_temp must be
414  // different than Y_local.
415  inverse_mv_type* Y_temp;
416  if (InverseScalar(beta) == STS::zero ()) {
417  Y_temp = &Y_local;
418  } else {
419  Y_temp = new inverse_mv_type(Inverses_[blockIndex]->getRangeMap(), numVecs);
420  }
421  // Apply the local operator: Y_temp := M^{-1} * X_scaled
422  Inverses_[blockIndex]->apply(X_scaled, *Y_temp, mode);
423  // Y_local := beta * Y_local + alpha * diag(D_local) * Y_tmp.
424  //
425  // Note that we still use the permuted subset scaling D_local here,
426  // because Y_temp has the same permuted subset Map. That's good, in
427  // fact, because it's a subset: less data to read and multiply.
428  Y_local.elementWiseMultiply(alpha, D_local, *Y_temp, beta);
429  if(Y_temp != &Y_local)
430  delete Y_temp;
431  // Copy the permuted subset output vector Y_local into the original
432  // output multivector Y.
433  mvgs.scatterMVtoView(Y, Y_local, blockRows);
434 }
435 
436 //==============================================================================
437 template<class MatrixType, class InverseType>
438 std::ostream& SparseContainer<MatrixType,InverseType>::print(std::ostream& os) const
439 {
440  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
441  fos.setOutputToRootOnly(0);
442  describe(fos);
443  return(os);
444 }
445 
446 //==============================================================================
447 template<class MatrixType, class InverseType>
449 {
450  std::ostringstream oss;
451  oss << "\"Ifpack2::SparseContainer\": {";
452  if (this->isInitialized()) {
453  if (this->isComputed()) {
454  oss << "status = initialized, computed";
455  }
456  else {
457  oss << "status = initialized, not computed";
458  }
459  }
460  else {
461  oss << "status = not initialized, not computed";
462  }
463  for(int i = 0; i < this->numBlocks_; i++)
464  {
465  oss << ", Block Inverse " << i << ": {";
466  oss << Inverses_[i]->description();
467  oss << "}";
468  }
469  oss << "}";
470  return oss.str();
471 }
472 
473 //==============================================================================
474 template<class MatrixType, class InverseType>
476 {
477  using std::endl;
478  if(verbLevel==Teuchos::VERB_NONE) return;
479  os << "================================================================================" << endl;
480  os << "Ifpack2::SparseContainer" << endl;
481  for(int i = 0; i < this->numBlocks_; i++)
482  {
483  os << "Block " << i << " rows: = " << this->blockSizes_[i] << endl;
484  }
485  os << "isInitialized() = " << this->IsInitialized_ << endl;
486  os << "isComputed() = " << this->IsComputed_ << endl;
487  os << "================================================================================" << endl;
488  os << endl;
489 }
490 
491 //==============================================================================
492 template<class MatrixType, class InverseType>
494 extract ()
495 {
496  using Teuchos::Array;
497  using Teuchos::ArrayView;
498  using Teuchos::RCP;
499  const LO INVALID = Teuchos::OrdinalTraits<LO>::invalid();
500  //To extract diagonal blocks, need to translate local rows to local columns.
501  //Strategy: make a lookup table that translates local cols in the matrix to offsets in blockRows_:
502  //blockOffsets_[b] <= offset < blockOffsets_[b+1]: tests whether the column is in block b.
503  //offset - blockOffsets_[b]: gives the column within block b.
504  //
505  //This provides the block and col within a block in O(1).
506  Teuchos::Array<InverseGlobalOrdinal> indicesToInsert;
507  Teuchos::Array<InverseScalar> valuesToInsert;
508  if(this->scalarsPerRow_ > 1)
509  {
510  Array<LO> colToBlockOffset(this->inputBlockMatrix_->getNodeNumCols(), INVALID);
511  for(int i = 0; i < this->numBlocks_; i++)
512  {
513  //Get the interval where block i is defined in blockRows_
514  LO blockStart = this->blockOffsets_[i];
515  LO blockSize = this->blockSizes_[i];
516  LO blockPointSize = this->blockSizes_[i] * this->scalarsPerRow_;
517  LO blockEnd = blockStart + blockSize;
518  ArrayView<const LO> blockRows = this->getBlockRows(i);
519  //Set the lookup table entries for the columns appearing in block i.
520  //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
521  //this is OK. The values updated here are only needed to process block i's entries.
522  for(LO j = 0; j < blockSize; j++)
523  {
524  LO localCol = this->translateRowToCol(blockRows[j]);
525  colToBlockOffset[localCol] = blockStart + j;
526  }
527  //First, count the number of entries in each row of block i
528  //(to allocate it with StaticProfile)
529  Array<size_t> rowEntryCounts(blockPointSize, 0);
530  //blockRow counts the BlockCrs LIDs that are going into this block
531  //Rows are inserted into the CrsMatrix in sequential order
532  for(LO blockRow = 0; blockRow < blockSize; blockRow++)
533  {
534  //get a raw view of the whole block row
535  const LO* indices;
536  SC* values;
537  LO numEntries;
538  LO inputRow = this->blockRows_[blockStart + blockRow];
539  this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values, numEntries);
540  for(LO br = 0; br < this->bcrsBlockSize_; br++)
541  {
542  for(LO k = 0; k < numEntries; k++)
543  {
544  LO colOffset = colToBlockOffset[indices[k]];
545  if(blockStart <= colOffset && colOffset < blockEnd)
546  {
547  rowEntryCounts[blockRow * this->bcrsBlockSize_ + br] += this->bcrsBlockSize_;
548  }
549  }
550  }
551  }
552  //now that row sizes are known, can allocate the diagonal matrix
553  RCP<InverseMap> tempMap(new InverseMap(blockPointSize, 0, this->localComm_));
554  diagBlocks_[i] = rcp(new InverseCrs(tempMap, rowEntryCounts, Tpetra::StaticProfile));
555  Inverses_[i] = rcp(new InverseType(diagBlocks_[i]));
556  //insert the actual entries, one row at a time
557  for(LO blockRow = 0; blockRow < blockSize; blockRow++)
558  {
559  //get a raw view of the whole block row
560  const LO* indices;
561  SC* values;
562  LO numEntries;
563  LO inputRow = this->blockRows_[blockStart + blockRow];
564  this->inputBlockMatrix_->getLocalRowView(inputRow, indices, values, numEntries);
565  for(LO br = 0; br < this->bcrsBlockSize_; br++)
566  {
567  indicesToInsert.clear();
568  valuesToInsert.clear();
569  for(LO k = 0; k < numEntries; k++)
570  {
571  LO colOffset = colToBlockOffset[indices[k]];
572  if(blockStart <= colOffset && colOffset < blockEnd)
573  {
574  LO blockCol = colOffset - blockStart;
575  //bc iterates over the columns in (block) entry k
576  for(LO bc = 0; bc < this->bcrsBlockSize_; bc++)
577  {
578  indicesToInsert.push_back(blockCol * this->bcrsBlockSize_ + bc);
579  valuesToInsert.push_back(values[k * this->bcrsBlockSize_ * this->bcrsBlockSize_ + bc * this->bcrsBlockSize_ + br]);
580  }
581  }
582  }
583  InverseGlobalOrdinal rowToInsert = blockRow * this->bcrsBlockSize_ + br;
584  if(indicesToInsert.size())
585  diagBlocks_[i]->insertGlobalValues(rowToInsert, indicesToInsert(), valuesToInsert());
586  }
587  }
588  diagBlocks_[i]->fillComplete();
589  }
590  }
591  else
592  {
593  //get the mapping from point-indexed matrix columns to offsets in blockRows_
594  //(this includes regular CrsMatrix columns, in which case scalarsPerRow_ == 1)
595  Array<LO> colToBlockOffset(this->inputMatrix_->getNodeNumCols() * this->bcrsBlockSize_, INVALID);
596  for(int i = 0; i < this->numBlocks_; i++)
597  {
598  //Get the interval where block i is defined in blockRows_
599  LO blockStart = this->blockOffsets_[i];
600  LO blockSize = this->blockSizes_[i];
601  LO blockEnd = blockStart + blockSize;
602  ArrayView<const LO> blockRows = this->getBlockRows(i);
603  //Set the lookup table entries for the columns appearing in block i.
604  //If OverlapLevel_ > 0, then this may overwrite values for previous blocks, but
605  //this is OK. The values updated here are only needed to process block i's entries.
606  for(LO j = 0; j < blockSize; j++)
607  {
608  //translateRowToCol will return the corresponding split column
609  LO localCol = this->translateRowToCol(blockRows[j]);
610  colToBlockOffset[localCol] = blockStart + j;
611  }
612  Teuchos::Array<size_t> rowEntryCounts(blockSize, 0);
613  for(LO j = 0; j < blockSize; j++)
614  {
615  rowEntryCounts[j] = this->getInputRowView(this->blockRows_[blockStart + j]).size();
616  }
617  RCP<InverseMap> tempMap(new InverseMap(blockSize, 0, this->localComm_));
618  diagBlocks_[i] = rcp(new InverseCrs(tempMap, rowEntryCounts, Tpetra::StaticProfile));
619  Inverses_[i] = rcp(new InverseType(diagBlocks_[i]));
620  for(LO blockRow = 0; blockRow < blockSize; blockRow++)
621  {
622  valuesToInsert.clear();
623  indicesToInsert.clear();
624  //get a view of the split row
625  LO inputSplitRow = this->blockRows_[blockStart + blockRow];
626  auto rowView = this->getInputRowView(inputSplitRow);
627  for(size_t k = 0; k < rowView.size(); k++)
628  {
629  LO colOffset = colToBlockOffset[rowView.ind(k)];
630  if(blockStart <= colOffset && colOffset < blockEnd)
631  {
632  LO blockCol = colOffset - blockStart;
633  indicesToInsert.push_back(blockCol);
634  valuesToInsert.push_back(rowView.val(k));
635  }
636  }
637  if(indicesToInsert.size())
638  diagBlocks_[i]->insertGlobalValues(blockRow, indicesToInsert(), valuesToInsert());
639  }
640  diagBlocks_[i]->fillComplete ();
641  }
642  }
643 }
644 
645 template<typename MatrixType, typename InverseType>
647 {
648  typedef ILUT<Tpetra::RowMatrix<SC, LO, GO, NO> > ILUTInverse;
649 #ifdef HAVE_IFPACK2_AMESOS2
651  if(std::is_same<InverseType, ILUTInverse>::value)
652  {
653  return "SparseILUT";
654  }
655  else if(std::is_same<InverseType, AmesosInverse>::value)
656  {
657  return "SparseAmesos";
658  }
659  else
660  {
661  throw std::logic_error("InverseType for SparseContainer must be Ifpack2::ILUT or Details::Amesos2Wrapper");
662  }
663 #else
664  // Macros can't have commas in their arguments, so we have to
665  // compute the bool first argument separately.
666  constexpr bool inverseTypeIsILUT = std::is_same<InverseType, ILUTInverse>::value;
667  TEUCHOS_TEST_FOR_EXCEPTION(! inverseTypeIsILUT, std::logic_error,
668  "InverseType for SparseContainer must be Ifpack2::ILUT<ROW>");
669  return "SparseILUT"; //the only supported sparse container specialization if no Amesos2
670 #endif
671 }
672 
673 } // namespace Ifpack2
674 
675 // For ETI
676 #include "Ifpack2_ILUT.hpp"
677 #ifdef HAVE_IFPACK2_AMESOS2
678 #include "Ifpack2_Details_Amesos2Wrapper.hpp"
679 #endif
680 
681 // There's no need to instantiate for CrsMatrix too. All Ifpack2
682 // preconditioners can and should do dynamic casts if they need a type
683 // more specific than RowMatrix.
684 
685 #ifdef HAVE_IFPACK2_AMESOS2
686 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
687  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
688  Ifpack2::ILUT<Tpetra::RowMatrix<S,LO,GO,N> > >; \
689  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S, LO, GO, N>, \
690  Ifpack2::Details::Amesos2Wrapper<Tpetra::RowMatrix<S,LO,GO,N> > >;
691 #else
692 # define IFPACK2_SPARSECONTAINER_INSTANT(S,LO,GO,N) \
693  template class Ifpack2::SparseContainer< Tpetra::RowMatrix<S,LO,GO,N>, \
694  Ifpack2::ILUT<Tpetra::RowMatrix<S, LO, GO, N> > >;
695 #endif
696 #endif // IFPACK2_SPARSECONTAINER_DEF_HPP
static std::string getName()
Get the name of this container type for Details::constructContainer()
Definition: Ifpack2_SparseContainer_def.hpp:646
SparseContainer(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< LO > > &partitions, const Teuchos::RCP< const import_type > &importer, bool pointIndexed)
Constructor.
Definition: Ifpack2_SparseContainer_def.hpp:60
Store and solve a local sparse linear problem.
Definition: Ifpack2_SparseContainer_decl.hpp:134
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:475
The implementation of the numerical features of Container (Jacobi, Gauss-Seidel, SGS). This class allows a custom scalar type (LocalScalarType) to be used for storing blocks and solving the block systems. Hiding this template parameter from the Container interface simplifies the BlockRelaxation and ContainerFactory classes.
Definition: Ifpack2_Container_decl.hpp:342
#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:86
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:438
virtual std::string description() const
A one-line description of this object.
Definition: Ifpack2_SparseContainer_def.hpp:448
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
virtual void weightedApply(HostView X, HostView Y, HostView W, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:288
virtual void apply(HostView X, HostView Y, int blockIndex, Teuchos::ETransp mode=Teuchos::NO_TRANS, SC alpha=Teuchos::ScalarTraits< SC >::one(), SC beta=Teuchos::ScalarTraits< SC >::zero()) const
Compute Y := alpha * M^{-1} X + beta*Y.
Definition: Ifpack2_SparseContainer_def.hpp:168
typename mv_type::dual_view_type::t_host HostView
Definition: Ifpack2_Container_decl.hpp:139
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void push_back(const value_type &x)
virtual ~SparseContainer()
Destructor (declared virtual for memory safety of derived classes).
Definition: Ifpack2_SparseContainer_def.hpp:75
size_type size() const
TypeTo as(const TypeFrom &t)
virtual void setParameters(const Teuchos::ParameterList &List)
Set all necessary parameters.
Definition: Ifpack2_SparseContainer_def.hpp:79
virtual void compute()
Initialize and compute all blocks.
Definition: Ifpack2_SparseContainer_def.hpp:105
Ifpack2::SparseContainer class declaration.
void clearBlocks()
Definition: Ifpack2_SparseContainer_def.hpp:130
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85