Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Container_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_CONTAINER_DEF_HPP
44 #define IFPACK2_CONTAINER_DEF_HPP
45 
47 #include <Teuchos_Time.hpp>
48 
49 namespace Ifpack2 {
50 
51 //Implementation of Ifpack2::Container
52 
53 template<class MatrixType>
56  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
57  bool pointIndexed) :
58  inputMatrix_ (matrix),
59  inputCrsMatrix_ (Teuchos::rcp_dynamic_cast<const crs_matrix_type>(inputMatrix_)),
60  inputBlockMatrix_ (Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(inputMatrix_)),
61  pointIndexed_(pointIndexed),
62  IsInitialized_(false),
63  IsComputed_(false)
64 {
65  using Teuchos::Ptr;
66  using Teuchos::RCP;
67  using Teuchos::Array;
68  using Teuchos::ArrayView;
69  using Teuchos::Comm;
70  NumLocalRows_ = inputMatrix_->getNodeNumRows();
71  NumGlobalRows_ = inputMatrix_->getGlobalNumRows();
72  NumGlobalNonzeros_ = inputMatrix_->getGlobalNumEntries();
73  IsParallel_ = inputMatrix_->getRowMap()->getComm()->getSize() != 1;
75  if(hasBlockCrs_)
76  bcrsBlockSize_ = inputBlockMatrix_->getBlockSize();
77  else
78  bcrsBlockSize_ = 1;
81  else
82  scalarsPerRow_ = 1;
83  setBlockSizes(partitions);
84  //Sanity check the partitions
85  #ifdef HAVE_IFPACK2_DEBUG
86  // Check whether the input set of local row indices is correct.
87  const map_type& rowMap = *inputMatrix_->getRowMap();
88  for(int i = 0; i < numBlocks_; i++)
89  {
91  for(LO j = 0; j < blockSizes_[i]; j++)
92  {
93  LO row = blockRows[j];
94  if(pointIndexed)
95  {
96  //convert the point row to the corresponding block row
97  row /= bcrsBlockSize_;
98  }
100  !rowMap.isNodeLocalElement(row),
101  std::invalid_argument, "Ifpack2::Container: "
102  "On process " << rowMap.getComm()->getRank() << " of "
103  << rowMap.getComm()->getSize() << ", in the given set of local row "
104  "indices blockRows = " << Teuchos::toString(blockRows) << ", the following "
105  "entries is not valid local row index on the calling process: "
106  << row << ".");
107  }
108  }
109  #endif
110 }
111 
112 template<class MatrixType>
115 
116 template<class MatrixType>
119 {
121  (&blockRows_[blockOffsets_[blockIndex]], blockSizes_[blockIndex]);
122 }
123 
124 template<class MatrixType>
126 {
127  //First, create a grand total of all the rows in all the blocks
128  //Note: If partitioner allowed overlap, this could be greater than the # of local rows
129  LO totalBlockRows = 0;
130  numBlocks_ = partitions.size();
131  blockSizes_.resize(numBlocks_);
132  blockOffsets_.resize(numBlocks_);
133  maxBlockSize_ = 0;
134  for(int i = 0; i < numBlocks_; i++)
135  {
136  LO rowsInBlock = partitions[i].size();
137  blockSizes_[i] = rowsInBlock;
138  blockOffsets_[i] = totalBlockRows;
139  totalBlockRows += rowsInBlock;
140  maxBlockSize_ = std::max(maxBlockSize_, rowsInBlock * scalarsPerRow_);
141  }
142  blockRows_.resize(totalBlockRows);
143  //set blockRows_: each entry is the partition/block of the row
144  LO iter = 0;
145  for(int i = 0; i < numBlocks_; i++)
146  {
147  for(int j = 0; j < blockSizes_[i]; j++)
148  {
149  blockRows_[iter++] = partitions[i][j];
150  }
151  }
152 }
153 
154 template<class MatrixType>
156 {
157  if(Diag_.is_null())
158  {
159  Diag_ = rcp(new vector_type(inputMatrix_->getDomainMap()));
160  inputMatrix_->getLocalDiagCopy(*Diag_);
161  }
162 }
163 
164 template<class MatrixType>
166  return IsInitialized_;
167 }
168 
169 template<class MatrixType>
171  return IsComputed_;
172 }
173 
174 template<class MatrixType>
176 applyMV(mv_type& X, mv_type& Y) const
177 {
178  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
179 }
180 
181 template<class MatrixType>
183 weightedApplyMV(mv_type& X, mv_type& Y, vector_type& W) const
184 {
185  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
186 }
187 
188 template<class MatrixType>
189 std::string Container<MatrixType>::
191 {
192  return "Generic";
193 }
194 
195 template<class MatrixType>
197  SC dampingFactor, LO i) const
198 {
199  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
200 }
201 
202 template <class MatrixType>
203 void Container<MatrixType>::DoJacobi(HostView X, HostView Y, SC dampingFactor) const
204 {
205  using STS = Teuchos::ScalarTraits<ISC>;
206  const ISC one = STS::one();
207  // use blockRows_ and blockSizes_
208  size_t numVecs = X.extent(1);
209  // Non-overlapping Jacobi
210  for (LO i = 0; i < numBlocks_; i++)
211  {
212  // may happen that a partition is empty
213  if(blockSizes_[i] != 1 || hasBlockCrs_)
214  {
215  if(blockSizes_[i] == 0 )
216  continue;
217  apply(X, Y, i, Teuchos::NO_TRANS, dampingFactor, one);
218  }
219  else // singleton, can't access Containers_[i] as it was never filled and may be null.
220  {
221  LO LRID = blockRows_[blockOffsets_[i]];
222  getMatDiag();
223  HostView diagView = Diag_->getLocalViewHost();
224  ISC d = one / diagView(LRID, 0);
225  for(size_t nv = 0; nv < numVecs; nv++)
226  {
227  ISC x = X(LRID, nv);
228  Y(LRID, nv) = x * d;
229  }
230  }
231  }
232 }
233 
234 template <class MatrixType>
235 void Container<MatrixType>::DoOverlappingJacobi(HostView X, HostView Y, HostView W, SC dampingFactor) const
236 {
237  using STS = Teuchos::ScalarTraits<SC>;
238  // Overlapping Jacobi
239  for(LO i = 0; i < numBlocks_; i++)
240  {
241  // may happen that a partition is empty
242  if(blockSizes_[i] == 0)
243  continue;
244  if(blockSizes_[i] != 1)
245  weightedApply(X, Y, W, i, Teuchos::NO_TRANS, dampingFactor, STS::one());
246  }
247 }
248 
249 //Do Gauss-Seidel with just block i
250 //This is used 3 times: once in DoGaussSeidel and twice in DoSGS
251 template<class MatrixType, typename LocalScalarType>
253  HostView X, HostView Y, HostView Y2, HostView Resid,
254  SC dampingFactor, LO i) const
255 {
256  using Teuchos::ArrayView;
257  using STS = Teuchos::ScalarTraits<ISC>;
258  size_t numVecs = X.extent(1);
259  const ISC one = STS::one();
260  if(this->blockSizes_[i] == 0)
261  return; // Skip empty partitions
262  if(this->hasBlockCrs_ && !this->pointIndexed_)
263  {
264  //Use efficient blocked version
265  ArrayView<const LO> blockRows = this->getBlockRows(i);
266  const size_t localNumRows = this->blockSizes_[i];
267  for(size_t j = 0; j < localNumRows; j++)
268  {
269  LO row = blockRows[j]; // Containers_[i]->ID (j);
270  LO numEntries;
271  SC* values;
272  const LO* colinds;
273  this->inputBlockMatrix_->getLocalRowView(row, colinds, values, numEntries);
274  for(size_t m = 0; m < numVecs; m++)
275  {
276  for (int localR = 0; localR < this->bcrsBlockSize_; localR++)
277  Resid(row * this->bcrsBlockSize_ + localR, m) = X(row * this->bcrsBlockSize_ + localR, m);
278  for (LO k = 0; k < numEntries; ++k)
279  {
280  const LO col = colinds[k];
281  for(int localR = 0; localR < this->bcrsBlockSize_; localR++)
282  {
283  for(int localC = 0; localC < this->bcrsBlockSize_; localC++)
284  {
285  Resid(row * this->bcrsBlockSize_ + localR, m) -=
286  values[k * this->bcrsBlockSize_ * this->bcrsBlockSize_ + localR + localC * this->bcrsBlockSize_]
287  * Y2(col * this->bcrsBlockSize_ + localC, m); }
288  }
289  }
290  }
291  }
292  // solve with this block
293  //
294  // Note: I'm abusing the ordering information, knowing that X/Y
295  // and Y2 have the same ordering for on-proc unknowns.
296  //
297  // Note: Add flop counts for inverse apply
298  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
299  }
300  else if(!this->hasBlockCrs_ && this->blockSizes_[i] == 1)
301  {
302  // singleton, can't access Containers_[i] as it was never filled and may be null.
303  // a singleton calculation (just using matrix diagonal) is exact, all residuals should be zero.
304  LO LRID = this->blockOffsets_[i]; // by definition, a singleton 1 row in block.
305  HostView diagView = this->Diag_->getLocalViewHost();
306  ISC d = one / diagView(LRID, 0);
307  for(size_t m = 0; m < numVecs; m++)
308  {
309  ISC x = X(LRID, m);
310  ISC newy = x * d;
311  Y2(LRID, m) = newy;
312  }
313  }
314  else if(!this->inputCrsMatrix_.is_null() &&
315  std::is_same<typename crs_matrix_type::device_type::memory_space, Kokkos::HostSpace>::value)
316  {
317  //Use the KokkosSparse internal matrix for low-overhead values/indices access
318  //But, can only do this if the matrix is accessible directly from host, since it's not a DualView
320  container_exec_space().fence();
321  auto localA = this->inputCrsMatrix_->getLocalMatrix();
322  using size_type = typename crs_matrix_type::local_matrix_type::size_type;
323  const auto& rowmap = localA.graph.row_map;
324  const auto& entries = localA.graph.entries;
325  const auto& values = localA.values;
326  ArrayView<const LO> blockRows = this->getBlockRows(i);
327  for(size_t j = 0; j < size_t(blockRows.size()); j++)
328  {
329  const LO row = blockRows[j];
330  for(size_t m = 0; m < numVecs; m++)
331  {
332  ISC r = X(row, m);
333  for(size_type k = rowmap(row); k < rowmap(row + 1); k++)
334  {
335  const LO col = entries(k);
336  r -= values(k) * Y2(col, m);
337  }
338  Resid(row, m) = r;
339  }
340  }
341  // solve with this block
342  //
343  // Note: I'm abusing the ordering information, knowing that X/Y
344  // and Y2 have the same ordering for on-proc unknowns.
345  //
346  // Note: Add flop counts for inverse apply
347  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
348  }
349  else
350  {
351  //Either a point-indexed block matrix, or a normal row matrix
352  //that doesn't support getLocalMatrix
353  ArrayView<const LO> blockRows = this->getBlockRows(i);
354  for(size_t j = 0; j < size_t(blockRows.size()); j++)
355  {
356  const LO row = blockRows[j];
357  auto rowView = getInputRowView(row);
358  for(size_t m = 0; m < numVecs; m++)
359  {
360  Resid(row, m) = X(row, m);
361  for (size_t k = 0; k < rowView.size(); ++k)
362  {
363  const LO col = rowView.ind(k);
364  Resid(row, m) -= rowView.val(k) * Y2(col, m);
365  }
366  }
367  }
368  // solve with this block
369  //
370  // Note: I'm abusing the ordering information, knowing that X/Y
371  // and Y2 have the same ordering for on-proc unknowns.
372  //
373  // Note: Add flop counts for inverse apply
374  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
375  }
376 }
377 
378 template<class MatrixType>
380 DoGaussSeidel(HostView X, HostView Y, HostView Y2, SC dampingFactor) const
381 {
382  using Teuchos::Array;
383  using Teuchos::ArrayRCP;
384  using Teuchos::ArrayView;
385  using Teuchos::Ptr;
386  using Teuchos::RCP;
387  using Teuchos::rcp;
388  using Teuchos::rcpFromRef;
389  //This function just extracts the diagonal if it hasn't already.
390  getMatDiag();
391  auto numVecs = X.extent(1);
392  // X = RHS, Y = initial guess
393  HostView Resid("", X.extent(0), X.extent(1));
394  for(LO i = 0; i < numBlocks_; i++)
395  {
396  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
397  }
398  if(IsParallel_)
399  {
400  auto numMyRows = inputMatrix_->getNodeNumRows();
401  for (size_t m = 0; m < numVecs; ++m)
402  {
403  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
404  {
405  Y(i, m) = Y2(i, m);
406  }
407  }
408  }
409 }
410 
411 template<class MatrixType>
412 void Container<MatrixType>::
413 DoSGS(HostView X, HostView Y, HostView Y2, SC dampingFactor) const
414 {
415  // X = RHS, Y = initial guess
416  using Teuchos::Array;
417  using Teuchos::ArrayRCP;
418  using Teuchos::ArrayView;
419  using Teuchos::Ptr;
420  using Teuchos::ptr;
421  using Teuchos::RCP;
422  using Teuchos::rcp;
423  using Teuchos::rcpFromRef;
424  auto numVecs = X.extent(1);
425  HostView Resid("", X.extent(0), X.extent(1));
426  // Forward Sweep
427  for(LO i = 0; i < numBlocks_; i++)
428  {
429  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
430  }
431  static_assert(std::is_signed<LO>::value,
432  "Local ordinal must be signed (unsigned breaks reverse iteration to 0)");
433  // Reverse Sweep
434  for(LO i = numBlocks_ - 1; i >= 0; --i)
435  {
436  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
437  }
438  if(IsParallel_)
439  {
440  auto numMyRows = inputMatrix_->getNodeNumRows();
441  for (size_t m = 0; m < numVecs; ++m)
442  {
443  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
444  {
445  Y(i, m) = Y2(i, m);
446  }
447  }
448  }
449 }
450 
451 template<class MatrixType>
452 void Container<MatrixType>::
453 clearBlocks()
454 {
455  numBlocks_ = 0;
456  blockRows_.clear();
457  blockSizes_.clear();
458  blockOffsets_.clear();
459  Diag_ = Teuchos::null; //Diag_ will be recreated if needed
460 }
461 
462 //Implementation of Ifpack2::ContainerImpl
463 
464 template<class MatrixType, class LocalScalarType>
465 ContainerImpl<MatrixType, LocalScalarType>::
466 ContainerImpl(
468  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
469  bool pointIndexed)
470  : Container<MatrixType>(matrix, partitions, pointIndexed) {}
471 
472 template<class MatrixType, class LocalScalarType>
475 
476 template<class MatrixType, class LocalScalarType>
479 
480 template<class MatrixType, class LocalScalarType>
482 applyInverseJacobi (const mv_type& /* X */, mv_type& /* Y */,
483  SC dampingFactor,
484  bool /* zeroStartingSolution = false */,
485  int /* numSweeps = 1 */) const
486 {
487  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
488 }
489 
490 template<class MatrixType, class LocalScalarType>
492 applyMV (mv_type& X, mv_type& Y) const
493 {
494  HostView XView = X.getLocalViewHost();
495  HostView YView = Y.getLocalViewHost();
496  this->apply (XView, YView, 0);
497 }
498 
499 template<class MatrixType, class LocalScalarType>
501 weightedApplyMV (mv_type& X,
502  mv_type& Y,
503  vector_type& W) const
504 {
505  HostView XView = X.getLocalViewHost();
506  HostView YView = Y.getLocalViewHost();
507  HostView WView = W.getLocalViewHost();
508  weightedApply (XView, YView, WView, 0);
509 }
510 
511 template<class MatrixType, class LocalScalarType>
514 {
515  return "Generic";
516 }
517 
518 template<class MatrixType, class LocalScalarType>
520 solveBlock(HostSubviewLocal X,
521  HostSubviewLocal Y,
522  int blockIndex,
523  Teuchos::ETransp mode,
524  const LSC alpha,
525  const LSC beta) const
526 {
527  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
528 }
529 
530 template<class MatrixType, class LocalScalarType>
531 typename ContainerImpl<MatrixType, LocalScalarType>::LO
534 {
535  LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
536  GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
537  const map_type& globalRowMap = *(this->inputMatrix_->getRowMap());
538  const map_type& globalColMap = *(this->inputMatrix_->getColMap());
539  LO rowLID = row;
540  LO dofOffset = 0;
541  if(this->pointIndexed_)
542  {
543  rowLID = row / this->bcrsBlockSize_;
544  dofOffset = row % this->bcrsBlockSize_;
545  }
546  GO diagGID = globalRowMap.getGlobalElement(rowLID);
548  diagGID == GO_INVALID,
549  std::runtime_error, "Ifpack2::Container::translateRowToCol: "
550  "On process " << this->inputMatrix_->getRowMap()->getComm()->getRank() <<
551  ", at least one row index in the set of local "
552  "row indices given to the constructor is not a valid local row index in "
553  "the input matrix's row Map on this process. This should be impossible "
554  "because the constructor checks for this case. Here is the complete set "
555  "of invalid local row indices: " << rowLID << ". "
556  "Please report this bug to the Ifpack2 developers.");
557  //now, can translate diagGID (both a global row AND global col ID) to local column
558  LO colLID = globalColMap.getLocalElement(diagGID);
560  colLID == LO_INVALID,
561  std::runtime_error, "Ifpack2::Container::translateRowToCol: "
562  "On process " << this->inputMatrix_->getRowMap()->getComm()->getRank() << ", "
563  "at least one row index in the set of row indices given to the constructor "
564  "does not have a corresponding column index in the input matrix's column "
565  "Map. This probably means that the column(s) in question is/are empty on "
566  "this process, which would make the submatrix to extract structurally "
567  "singular. The invalid global column index is " << diagGID << ".");
568  //colLID could identify a block column - translate to split column if needed
569  if(this->pointIndexed_)
570  return colLID * this->bcrsBlockSize_ + dofOffset;
571  return colLID;
572 }
573 
574 template<class MatrixType, class LocalScalarType>
577  HostView Y,
578  int blockIndex,
579  Teuchos::ETransp mode,
580  SC alpha,
581  SC beta) const
582 {
583  using Teuchos::ArrayView;
584  using Teuchos::as;
585  using Teuchos::RCP;
586  using Teuchos::rcp;
587 
588  // The local operator might have a different Scalar type than
589  // MatrixType. This means that we might have to convert X and Y to
590  // the Tpetra::MultiVector specialization that the local operator
591  // wants. This class' X_ and Y_ internal fields are of the right
592  // type for the local operator, so we can use those as targets.
593 
595 
597  ! this->IsComputed_, std::runtime_error, "Ifpack2::Container::apply: "
598  "You must have called the compute() method before you may call apply(). "
599  "You may call the apply() method as many times as you want after calling "
600  "compute() once, but you must have called compute() at least once.");
601 
602  const size_t numVecs = X.extent(1);
603 
604  if(numVecs == 0) {
605  return; // done! nothing to do
606  }
607 
608  // The local operator works on a permuted subset of the local parts
609  // of X and Y. The subset and permutation are defined by the index
610  // array returned by getBlockRows(). If the permutation is trivial
611  // and the subset is exactly equal to the local indices, then we
612  // could use the local parts of X and Y exactly, without needing to
613  // permute. Otherwise, we have to use temporary storage to permute
614  // X and Y. For now, we always use temporary storage.
615  //
616  // Create temporary permuted versions of the input and output.
617  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
618  // store the permuted versions of X resp. Y. Note that X_local has
619  // the domain Map of the operator, which may be a permuted subset of
620  // the local Map corresponding to X.getMap(). Similarly, Y_local
621  // has the range Map of the operator, which may be a permuted subset
622  // of the local Map corresponding to Y.getMap(). numRows_ here
623  // gives the number of rows in the row Map of the local Inverse_
624  // operator.
625  //
626  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
627  // here that the row Map and the range Map of that operator are
628  // the same.
629  //
630  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
631  // really belongs in Tpetra.
632 
633  if(X_localBlocks_.size() == 0 || X.extent(1) != X_local_.extent(1))
634  {
635  //need to resize (or create for the first time) the three scratch arrays
636  X_localBlocks_.clear();
637  Y_localBlocks_.clear();
638  X_localBlocks_.reserve(this->numBlocks_);
639  Y_localBlocks_.reserve(this->numBlocks_);
640 
641  X_local_ = HostViewLocal("X_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
642  Y_local_ = HostViewLocal("Y_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
643 
644  //create all X_local and Y_local managed Views at once, are
645  //reused in subsequent apply() calls
646  for(int i = 0; i < this->numBlocks_; i++)
647  {
648  auto blockBounds = std::make_pair(this->blockOffsets_[i] * this->scalarsPerRow_,
649  (this->blockOffsets_[i] + this->blockSizes_[i]) * this->scalarsPerRow_);
650  X_localBlocks_.emplace_back(X_local_, blockBounds, Kokkos::ALL());
651  Y_localBlocks_.emplace_back(Y_local_, blockBounds, Kokkos::ALL());
652  }
653  }
654 
655  const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
656 
657  if(this->scalarsPerRow_ == 1)
658  mvgs.gatherViewToView (X_localBlocks_[blockIndex], X, blockRows);
659  else
660  mvgs.gatherViewToViewBlock (X_localBlocks_[blockIndex], X, blockRows, this->scalarsPerRow_);
661 
662  // We must gather the contents of the output multivector Y even on
663  // input to solveBlock(), since the inverse operator might use it as
664  // an initial guess for a linear solve. We have no way of knowing
665  // whether it does or does not.
666 
667  if(this->scalarsPerRow_ == 1)
668  mvgs.gatherViewToView (Y_localBlocks_[blockIndex], Y, blockRows);
669  else
670  mvgs.gatherViewToViewBlock (Y_localBlocks_[blockIndex], Y, blockRows, this->scalarsPerRow_);
671 
672  // Apply the local operator:
673  // Y_local := beta*Y_local + alpha*M^{-1}*X_local
674  this->solveBlock (X_localBlocks_[blockIndex], Y_localBlocks_[blockIndex], blockIndex, mode,
675  LSC(alpha), LSC(beta));
676 
677  // Scatter the permuted subset output vector Y_local back into the
678  // original output multivector Y.
679  if(this->scalarsPerRow_ == 1)
680  mvgs.scatterViewToView (Y, Y_localBlocks_[blockIndex], blockRows);
681  else
682  mvgs.scatterViewToViewBlock (Y, Y_localBlocks_[blockIndex], blockRows, this->scalarsPerRow_);
683 }
684 
685 template<class MatrixType, class LocalScalarType>
688  HostView Y,
689  HostView D,
690  int blockIndex,
691  Teuchos::ETransp mode,
692  SC alpha,
693  SC beta) const
694 {
695  using Teuchos::ArrayRCP;
696  using Teuchos::ArrayView;
697  using Teuchos::Range1D;
698  using Teuchos::Ptr;
699  using Teuchos::ptr;
700  using Teuchos::RCP;
701  using Teuchos::rcp;
702  using Teuchos::rcp_const_cast;
703  using std::endl;
704  using STS = Teuchos::ScalarTraits<SC>;
705 
706  // The local operator template parameter might have a different
707  // Scalar type than MatrixType. This means that we might have to
708  // convert X and Y to the Tpetra::MultiVector specialization that
709  // the local operator wants. This class' X_ and Y_ internal fields
710  // are of the right type for the local operator, so we can use those
711  // as targets.
712 
713  const char prefix[] = "Ifpack2::Container::weightedApply: ";
715  ! this->IsComputed_, std::runtime_error, prefix << "You must have called the "
716  "compute() method before you may call this method. You may call "
717  "weightedApply() as many times as you want after calling compute() once, "
718  "but you must have called compute() at least once first.");
719 
720  //bmk 7-2019: BlockRelaxation already checked this, but if that changes...
722  this->scalarsPerRow_ > 1, std::logic_error, prefix << "Use of block rows isn't allowed "
723  "in overlapping Jacobi (the only method that uses weightedApply");
724 
725  const size_t numVecs = X.extent(1);
726 
728  X.extent(1) != Y.extent(1), std::runtime_error,
729  prefix << "X and Y have different numbers of vectors (columns). X has "
730  << X.extent(1) << ", but Y has " << Y.extent(1) << ".");
731 
732  if(numVecs == 0) {
733  return; // done! nothing to do
734  }
735 
736  const size_t numRows = this->blockSizes_[blockIndex];
737 
738  // The local operator works on a permuted subset of the local parts
739  // of X and Y. The subset and permutation are defined by the index
740  // array returned by getBlockRows(). If the permutation is trivial
741  // and the subset is exactly equal to the local indices, then we
742  // could use the local parts of X and Y exactly, without needing to
743  // permute. Otherwise, we have to use temporary storage to permute
744  // X and Y. For now, we always use temporary storage.
745  //
746  // Create temporary permuted versions of the input and output.
747  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
748  // store the permuted versions of X resp. Y. Note that X_local has
749  // the domain Map of the operator, which may be a permuted subset of
750  // the local Map corresponding to X.getMap(). Similarly, Y_local
751  // has the range Map of the operator, which may be a permuted subset
752  // of the local Map corresponding to Y.getMap(). numRows_ here
753  // gives the number of rows in the row Map of the local operator.
754  //
755  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
756  // here that the row Map and the range Map of that operator are
757  // the same.
758  //
759  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
760  // really belongs in Tpetra.
761  if(X_localBlocks_.size() == 0 || X.extent(1) != X_local_.extent(1))
762  {
763  //need to resize (or create for the first time) the three scratch arrays
764  X_localBlocks_.clear();
765  Y_localBlocks_.clear();
766  X_localBlocks_.reserve(this->numBlocks_);
767  Y_localBlocks_.reserve(this->numBlocks_);
768 
769  X_local_ = HostViewLocal("X_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
770  Y_local_ = HostViewLocal("Y_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
771 
772  //create all X_local and Y_local managed Views at once, are
773  //reused in subsequent apply() calls
774  for(int i = 0; i < this->numBlocks_; i++)
775  {
776  auto blockBounds = std::make_pair(this->blockOffsets_[i] * this->scalarsPerRow_,
777  (this->blockOffsets_[i] + this->blockSizes_[i]) * this->scalarsPerRow_);
778  X_localBlocks_.emplace_back(X_local_, blockBounds, Kokkos::ALL());
779  Y_localBlocks_.emplace_back(Y_local_, blockBounds, Kokkos::ALL());
780  }
781  }
782  if(int(weightedApplyScratch_.extent(0)) != 3 * this->maxBlockSize_ ||
783  weightedApplyScratch_.extent(1) != numVecs)
784  {
785  weightedApplyScratch_ = HostViewLocal("weightedApply scratch", 3 * this->maxBlockSize_, numVecs);
786  }
787 
788  ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
789 
791 
792  //note: BlockCrs w/ weighted Jacobi isn't allowed, so no need to use block gather/scatter
793  mvgs.gatherViewToView (X_localBlocks_[blockIndex], X, blockRows);
794  // We must gather the output multivector Y even on input to
795  // solveBlock(), since the local operator might use it as an initial
796  // guess for a linear solve. We have no way of knowing whether it
797  // does or does not.
798 
799  mvgs.gatherViewToView (Y_localBlocks_[blockIndex], Y, blockRows);
800 
801  // Apply the diagonal scaling D to the input X. It's our choice
802  // whether the result has the original input Map of X, or the
803  // permuted subset Map of X_local. If the latter, we also need to
804  // gather D into the permuted subset Map. We choose the latter, to
805  // save memory and computation. Thus, we do the following:
806  //
807  // 1. Gather D into a temporary vector D_local.
808  // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
809  // 3. Compute X_scaled := diag(D_loca) * X_local.
810  auto maxBS = this->maxBlockSize_;
811  auto bs = this->blockSizes_[blockIndex] * this->scalarsPerRow_;
812 
813  HostSubviewLocal D_local(weightedApplyScratch_, std::make_pair(0, bs), std::make_pair(0, 1));
814  mvgs.gatherViewToView (D_local, D, blockRows);
815  HostSubviewLocal X_scaled(weightedApplyScratch_, std::make_pair(maxBS, maxBS + bs), Kokkos::ALL());
816  for(size_t j = 0; j < numVecs; j++)
817  for(size_t i = 0; i < numRows; i++)
818  X_scaled(i, j) = X_localBlocks_[blockIndex](i, j) * D_local(i, 0);
819 
820  HostSubviewLocal Y_temp(weightedApplyScratch_, std::make_pair(maxBS * 2, maxBS * 2 + bs), Kokkos::ALL());
821  // Apply the local operator: Y_temp := M^{-1} * X_scaled
822  this->solveBlock (X_scaled, Y_temp, blockIndex, mode, STS::one(), STS::zero());
823  // Y_local := beta * Y_local + alpha * diag(D_local) * Y_temp.
824  //
825  // Note that we still use the permuted subset scaling D_local here,
826  // because Y_temp has the same permuted subset Map. That's good, in
827  // fact, because it's a subset: less data to read and multiply.
828  LISC a = alpha;
829  LISC b = beta;
830  for(size_t j = 0; j < numVecs; j++)
831  for(size_t i = 0; i < numRows; i++)
832  Y_localBlocks_[blockIndex](i, j) = b * Y_localBlocks_[blockIndex](i, j) + a * Y_temp(i, j) * D_local(i, 0);
833 
834  // Copy the permuted subset output vector Y_local into the original
835  // output multivector Y.
836  mvgs.scatterViewToView (Y, Y_localBlocks_[blockIndex], blockRows);
837 }
838 
839 template<class MatrixType, class LocalScalarType>
841  typename ContainerImpl<MatrixType, LocalScalarType>::SC,
842  typename ContainerImpl<MatrixType, LocalScalarType>::LO,
843  typename ContainerImpl<MatrixType, LocalScalarType>::GO,
844  typename ContainerImpl<MatrixType, LocalScalarType>::NO>
846 getInputRowView(LO row) const
847 {
848  if(this->hasBlockCrs_)
849  {
850  const LO* colinds;
851  SC* values;
852  LO numEntries;
853  this->inputBlockMatrix_->getLocalRowView(row / this->bcrsBlockSize_, colinds, values, numEntries);
854  return StridedRowView(values + row % this->bcrsBlockSize_, colinds, this->bcrsBlockSize_, numEntries * this->bcrsBlockSize_);
855  }
856  else if(!this->inputMatrix_->supportsRowViews())
857  {
858  size_t maxEntries = this->inputMatrix_->getNodeMaxNumRowEntries();
859  Teuchos::Array<LO> indsCopy(maxEntries);
860  Teuchos::Array<SC> valsCopy(maxEntries);
861  size_t numEntries;
862  this->inputMatrix_->getLocalRowCopy(row, indsCopy, valsCopy, numEntries);
863  indsCopy.resize(numEntries);
864  valsCopy.resize(numEntries);
865  return StridedRowView(valsCopy, indsCopy);
866  }
867  else
868  {
869  const LO* colinds;
870  const SC* values;
871  LO numEntries;
872  this->inputMatrix_->getLocalRowViewRaw(row, numEntries, colinds, values);
873  return StridedRowView(values, colinds, 1, numEntries);
874  }
875 }
876 
877 template<class MatrixType, class LocalScalarType>
879 clearBlocks()
880 {
881  X_localBlocks_.clear();
882  Y_localBlocks_.clear();
883  X_local_ = HostViewLocal();
884  Y_local_ = HostViewLocal();
886 }
887 
888 namespace Details {
889 
890 //Implementation of Ifpack2::Details::StridedRowView
891 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
893 StridedRowView(const SC* vals_, const LO* inds_, int blockSize_, size_t nnz_)
894  : vals(vals_), inds(inds_), blockSize(blockSize_), nnz(nnz_)
895 {}
896 
897 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
900  : vals(nullptr), inds(nullptr), blockSize(1), nnz(vals_.size())
901 {
902  valsCopy.swap(vals_);
903  indsCopy.swap(inds_);
904 }
905 
906 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
908 val(size_t i) const
909 {
910  #ifdef HAVE_IFPACK2_DEBUG
911  TEUCHOS_TEST_FOR_EXCEPTION(i >= nnz, std::runtime_error,
912  "Out-of-bounds access into Ifpack2::Container::StridedRowView");
913  #endif
914  if(vals)
915  {
916  if(blockSize == 1)
917  return vals[i];
918  else
919  return vals[i * blockSize];
920  }
921  else
922  return valsCopy[i];
923 }
924 
925 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
926 LocalOrdinal StridedRowView<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
927 ind(size_t i) const
928 {
929  #ifdef HAVE_IFPACK2_DEBUG
930  TEUCHOS_TEST_FOR_EXCEPTION(i >= nnz, std::runtime_error,
931  "Out-of-bounds access into Ifpack2::Container::StridedRowView");
932  #endif
933  //inds is smaller than vals by a factor of the block size (dofs/node)
934  if(inds)
935  {
936  if(blockSize == 1)
937  return inds[i];
938  else
939  return inds[i / blockSize] * blockSize + i % blockSize;
940  }
941  else
942  return indsCopy[i];
943 }
944 
945 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
946 size_t StridedRowView<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
947 size() const
948 {
949  return nnz;
950 }
951 }
952 
953 }
954 
955 template <class MatrixType>
956 std::ostream& operator<<(std::ostream& os, const Ifpack2::Container<MatrixType>& obj)
957 {
958  return obj.print(os);
959 }
960 
961 #define IFPACK2_CONTAINER_INSTANT(S,LO,GO,N) \
962  template class Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N>>; \
963  template class Ifpack2::ContainerImpl<Tpetra::RowMatrix<S, LO, GO, N>, S>; \
964  template class Ifpack2::Details::StridedRowView<S, LO, GO, N>; \
965  template std::ostream& operator<< <Tpetra::RowMatrix<S, LO, GO, N>>( \
966  std::ostream& os, const Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N>>& obj);
967 
968 #endif
969 
bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_Container_def.hpp:170
void setBlockSizes(const Teuchos::Array< Teuchos::Array< LO > > &partitions)
Initialize arrays with information about block sizes.
Definition: Ifpack2_Container_def.hpp:125
Teuchos::Array< LO > blockSizes_
Number of rows in each block.
Definition: Ifpack2_Container_decl.hpp:295
void applyMV(mv_type &X, mv_type &Y) const
Wrapper for apply with MVs, used in unit tests (never called by BlockRelaxation)
Definition: Ifpack2_Container_def.hpp:492
bool pointIndexed_
(If hasBlockCrs_) Whether the blocks are described using sub-block row indices instead of full block ...
Definition: Ifpack2_Container_decl.hpp:313
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container_decl.hpp:291
void swap(Array &x)
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 setParameters(const Teuchos::ParameterList &List)
Set parameters, if any.
Definition: Ifpack2_Container_def.hpp:478
bool isInitialized() const
Whether the container has been successfully initialized.
Definition: Ifpack2_Container_def.hpp:165
LO scalarsPerRow_
Definition: Ifpack2_Container_decl.hpp:316
virtual ~Container()
Destructor.
Definition: Ifpack2_Container_def.hpp:114
bool hasBlockCrs_
Whether the input matrix is a BlockCRS matrix.
Definition: Ifpack2_Container_decl.hpp:309
Details::StridedRowView< SC, LO, GO, NO > getInputRowView(LO row) const
View a row of the input matrix.
Definition: Ifpack2_Container_def.hpp:846
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_Container_def.hpp:576
LocalScalarType LSC
The internal representation of LocalScalarType in Kokkos::View.
Definition: Ifpack2_Container_decl.hpp:362
Declaration and definition of the Ifpack2::Details::MultiVectorLocalGatherScatter class...
typename Kokkos::Details::ArithTraits< SC >::val_type ISC
Internal representation of Scalar in Kokkos::View.
Definition: Ifpack2_Container_decl.hpp:135
Container(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< LO > > &partitions, bool pointIndexed)
Constructor.
Definition: Ifpack2_Container_def.hpp:54
int bcrsBlockSize_
If hasBlockCrs_, the number of DOFs per vertex. Otherwise 1.
Definition: Ifpack2_Container_decl.hpp:311
virtual void DoGSBlock(HostView X, HostView Y, HostView Y2, HostView Resid, SC dampingFactor, LO i) const
Do one step of Gauss-Seidel on block i (used by DoGaussSeidel and DoSGS)
Definition: Ifpack2_Container_def.hpp:196
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
StridedRowView(const SC *vals_, const LO *inds_, int blockSize_, size_t nnz_)
Constructor for row views (preferred)
Definition: Ifpack2_Container_def.hpp:893
LO translateRowToCol(LO row)
Definition: Ifpack2_Container_def.hpp:533
Teuchos::RCP< const block_crs_matrix_type > inputBlockMatrix_
The input matrix, dynamic cast to BlockCrsMatrix. May be null.
Definition: Ifpack2_Container_decl.hpp:288
Teuchos::RCP< const row_matrix_type > inputMatrix_
The input matrix to the constructor.
Definition: Ifpack2_Container_decl.hpp:282
virtual void applyInverseJacobi(const mv_type &, mv_type &, SC dampingFactor, bool, int) const
Compute Y := (1 - a) Y + a D^{-1} (X - R*Y).
Definition: Ifpack2_Container_def.hpp:482
typename mv_type::dual_view_type::t_host HostView
Definition: Ifpack2_Container_decl.hpp:139
virtual void weightedApplyMV(mv_type &X, mv_type &Y, vector_type &W) const
Wrapper for weightedApply with MultiVector.
Definition: Ifpack2_Container_def.hpp:183
virtual void solveBlock(HostSubviewLocal X, HostSubviewLocal Y, int blockIndex, Teuchos::ETransp mode, const LSC alpha, const LSC beta) const
Definition: Ifpack2_Container_def.hpp:520
void resize(size_type new_size, const value_type &x=value_type())
virtual void applyMV(mv_type &X, mv_type &Y) const
Wrapper for apply with MultiVector.
Definition: Ifpack2_Container_def.hpp:176
virtual ~ContainerImpl()
Destructor.
Definition: Ifpack2_Container_def.hpp:474
bool IsParallel_
Whether the problem is distributed across multiple MPI processes.
Definition: Ifpack2_Container_decl.hpp:301
void DoGSBlock(HostView X, HostView Y, HostView Y2, HostView Resid, SC dampingFactor, LO i) const
Do one step of Gauss-Seidel on block i (used by DoGaussSeidel and DoSGS)
Definition: Ifpack2_Container_def.hpp:252
virtual void weightedApply(HostView X, HostView Y, HostView D, 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_Container_def.hpp:687
static std::string getName()
Definition: Ifpack2_Container_def.hpp:513
TypeTo as(const TypeFrom &t)
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:112
Teuchos::ArrayView< const LO > getBlockRows(int blockIndex) const
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container_def.hpp:118
Structure for read-only views of general matrix rows.
Definition: Ifpack2_Container_decl.hpp:329
void weightedApplyMV(mv_type &X, mv_type &Y, vector_type &W) const
Wrapper for weightedApply with MVs, used in unit tests (never called by BlockRelaxation) ...
Definition: Ifpack2_Container_def.hpp:501
static std::string getName()
Definition: Ifpack2_Container_def.hpp:190
GO NumGlobalRows_
Number of global rows in input matrix.
Definition: Ifpack2_Container_decl.hpp:305
GO NumGlobalNonzeros_
Number of nonzeros in input matrix.
Definition: Ifpack2_Container_decl.hpp:307
LO NumLocalRows_
Number of local rows in input matrix.
Definition: Ifpack2_Container_decl.hpp:303
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85
std::string toString(const T &t)
bool is_null() const