Ifpack2 Templated Preconditioning Package  Version 1.0
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>
55  const Teuchos::RCP<const row_matrix_type>& matrix,
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;
74  hasBlockCrs_ = !inputBlockMatrix_.is_null();
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  {
90  Teuchos::ArrayView<const LO> blockRows = getBlockRows(i);
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  }
99  TEUCHOS_TEST_FOR_EXCEPTION(
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>
117 Teuchos::ArrayView<const typename MatrixType::local_ordinal_type>
119 {
120  return Teuchos::ArrayView<const LO>
121  (&blockRows_[blockOffsets_[blockIndex]], blockSizes_[blockIndex]);
122 }
123 
124 template<class MatrixType>
125 void Container<MatrixType>::setBlockSizes(const Teuchos::Array<Teuchos::Array<LO> >& partitions)
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(const 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(const 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>
196 void Container<MatrixType>::DoGSBlock(ConstHostView X, HostView Y, HostView Y2, HostView Resid,
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(ConstHostView 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  auto diagView = Diag_->getLocalViewHost(Tpetra::Access::ReadOnly);
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(ConstHostView X, HostView Y, ConstHostView 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  ConstHostView 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  using inds_type = typename block_crs_matrix_type::local_inds_host_view_type;
268  using vals_type = typename block_crs_matrix_type::values_host_view_type;
269  for(size_t j = 0; j < localNumRows; j++)
270  {
271  LO row = blockRows[j]; // Containers_[i]->ID (j);
272  vals_type values;
273  inds_type colinds;
274  this->inputBlockMatrix_->getLocalRowView(row, colinds, values);
275  LO numEntries = (LO) colinds.size();
276  for(size_t m = 0; m < numVecs; m++)
277  {
278  for (int localR = 0; localR < this->bcrsBlockSize_; localR++)
279  Resid(row * this->bcrsBlockSize_ + localR, m) = X(row * this->bcrsBlockSize_ + localR, m);
280  for (LO k = 0; k < numEntries; ++k)
281  {
282  const LO col = colinds[k];
283  for(int localR = 0; localR < this->bcrsBlockSize_; localR++)
284  {
285  for(int localC = 0; localC < this->bcrsBlockSize_; localC++)
286  {
287  Resid(row * this->bcrsBlockSize_ + localR, m) -=
288  values[k * this->bcrsBlockSize_ * this->bcrsBlockSize_ + localR + localC * this->bcrsBlockSize_]
289  * Y2(col * this->bcrsBlockSize_ + localC, m); }
290  }
291  }
292  }
293  }
294  // solve with this block
295  //
296  // Note: I'm abusing the ordering information, knowing that X/Y
297  // and Y2 have the same ordering for on-proc unknowns.
298  //
299  // Note: Add flop counts for inverse apply
300  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
301  }
302  else if(!this->hasBlockCrs_ && this->blockSizes_[i] == 1)
303  {
304  // singleton, can't access Containers_[i] as it was never filled and may be null.
305  // a singleton calculation (just using matrix diagonal) is exact, all residuals should be zero.
306  LO LRID = this->blockOffsets_[i]; // by definition, a singleton 1 row in block.
307  ConstHostView diagView = this->Diag_->getLocalViewHost(Tpetra::Access::ReadOnly);
308  ISC d = one / diagView(LRID, 0);
309  for(size_t m = 0; m < numVecs; m++)
310  {
311  ISC x = X(LRID, m);
312  ISC newy = x * d;
313  Y2(LRID, m) = newy;
314  }
315  }
316  else if(!this->inputCrsMatrix_.is_null() &&
317  std::is_same<typename crs_matrix_type::device_type::memory_space, Kokkos::HostSpace>::value)
318  {
319  //Use the KokkosSparse internal matrix for low-overhead values/indices access
320  //But, can only do this if the matrix is accessible directly from host, since it's not a DualView
322  container_exec_space().fence();
323  auto localA = this->inputCrsMatrix_->getLocalMatrixHost();
324  using size_type = typename crs_matrix_type::local_matrix_host_type::size_type;
325  const auto& rowmap = localA.graph.row_map;
326  const auto& entries = localA.graph.entries;
327  const auto& values = localA.values;
328  ArrayView<const LO> blockRows = this->getBlockRows(i);
329  for(size_t j = 0; j < size_t(blockRows.size()); j++)
330  {
331  const LO row = blockRows[j];
332  for(size_t m = 0; m < numVecs; m++)
333  {
334  ISC r = X(row, m);
335  for(size_type k = rowmap(row); k < rowmap(row + 1); k++)
336  {
337  const LO col = entries(k);
338  r -= values(k) * Y2(col, m);
339  }
340  Resid(row, m) = r;
341  }
342  }
343  // solve with this block
344  //
345  // Note: I'm abusing the ordering information, knowing that X/Y
346  // and Y2 have the same ordering for on-proc unknowns.
347  //
348  // Note: Add flop counts for inverse apply
349  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
350  }
351  else
352  {
353  //Either a point-indexed block matrix, or a normal row matrix
354  //that doesn't support getLocalMatrix
355  ArrayView<const LO> blockRows = this->getBlockRows(i);
356  for(size_t j = 0; j < size_t(blockRows.size()); j++)
357  {
358  const LO row = blockRows[j];
359  auto rowView = getInputRowView(row);
360  for(size_t m = 0; m < numVecs; m++)
361  {
362  Resid(row, m) = X(row, m);
363  for (size_t k = 0; k < rowView.size(); ++k)
364  {
365  const LO col = rowView.ind(k);
366  Resid(row, m) -= rowView.val(k) * Y2(col, m);
367  }
368  }
369  }
370  // solve with this block
371  //
372  // Note: I'm abusing the ordering information, knowing that X/Y
373  // and Y2 have the same ordering for on-proc unknowns.
374  //
375  // Note: Add flop counts for inverse apply
376  apply(Resid, Y2, i, Teuchos::NO_TRANS, dampingFactor, one);
377  }
378 }
379 
380 template<class MatrixType>
382 DoGaussSeidel(ConstHostView X, HostView Y, HostView Y2, SC dampingFactor) const
383 {
384  using Teuchos::Array;
385  using Teuchos::ArrayRCP;
386  using Teuchos::ArrayView;
387  using Teuchos::Ptr;
388  using Teuchos::RCP;
389  using Teuchos::rcp;
390  using Teuchos::rcpFromRef;
391  //This function just extracts the diagonal if it hasn't already.
392  getMatDiag();
393  auto numVecs = X.extent(1);
394  // X = RHS, Y = initial guess
395  HostView Resid("", X.extent(0), X.extent(1));
396  for(LO i = 0; i < numBlocks_; i++)
397  {
398  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
399  }
400  if(IsParallel_)
401  {
402  auto numMyRows = inputMatrix_->getNodeNumRows();
403  for (size_t m = 0; m < numVecs; ++m)
404  {
405  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
406  {
407  Y(i, m) = Y2(i, m);
408  }
409  }
410  }
411 }
412 
413 template<class MatrixType>
414 void Container<MatrixType>::
415 DoSGS(ConstHostView X, HostView Y, HostView Y2, SC dampingFactor) const
416 {
417  // X = RHS, Y = initial guess
418  using Teuchos::Array;
419  using Teuchos::ArrayRCP;
420  using Teuchos::ArrayView;
421  using Teuchos::Ptr;
422  using Teuchos::ptr;
423  using Teuchos::RCP;
424  using Teuchos::rcp;
425  using Teuchos::rcpFromRef;
426  auto numVecs = X.extent(1);
427  HostView Resid("", X.extent(0), X.extent(1));
428  // Forward Sweep
429  for(LO i = 0; i < numBlocks_; i++)
430  {
431  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
432  }
433  static_assert(std::is_signed<LO>::value,
434  "Local ordinal must be signed (unsigned breaks reverse iteration to 0)");
435  // Reverse Sweep
436  for(LO i = numBlocks_ - 1; i >= 0; --i)
437  {
438  DoGSBlock(X, Y, Y2, Resid, dampingFactor, i);
439  }
440  if(IsParallel_)
441  {
442  auto numMyRows = inputMatrix_->getNodeNumRows();
443  for (size_t m = 0; m < numVecs; ++m)
444  {
445  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
446  {
447  Y(i, m) = Y2(i, m);
448  }
449  }
450  }
451 }
452 
453 template<class MatrixType>
454 void Container<MatrixType>::
455 clearBlocks()
456 {
457  numBlocks_ = 0;
458  blockRows_.clear();
459  blockSizes_.clear();
460  blockOffsets_.clear();
461  Diag_ = Teuchos::null; //Diag_ will be recreated if needed
462 }
463 
464 //Implementation of Ifpack2::ContainerImpl
465 
466 template<class MatrixType, class LocalScalarType>
467 ContainerImpl<MatrixType, LocalScalarType>::
468 ContainerImpl(
469  const Teuchos::RCP<const row_matrix_type>& matrix,
470  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
471  bool pointIndexed)
472  : Container<MatrixType>(matrix, partitions, pointIndexed) {}
473 
474 template<class MatrixType, class LocalScalarType>
475 ContainerImpl<MatrixType, LocalScalarType>::
476 ~ContainerImpl() {}
477 
478 template<class MatrixType, class LocalScalarType>
480 setParameters (const Teuchos::ParameterList& List) {}
481 
482 template<class MatrixType, class LocalScalarType>
484 applyInverseJacobi (const mv_type& /* X */, mv_type& /* Y */,
485  SC dampingFactor,
486  bool /* zeroStartingSolution = false */,
487  int /* numSweeps = 1 */) const
488 {
489  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
490 }
491 
492 template<class MatrixType, class LocalScalarType>
494 applyMV (const mv_type& X, mv_type& Y) const
495 {
496  ConstHostView XView = X.getLocalViewHost(Tpetra::Access::ReadOnly);
497  HostView YView = Y.getLocalViewHost(Tpetra::Access::ReadWrite);
498  this->apply (XView, YView, 0);
499 }
500 
501 template<class MatrixType, class LocalScalarType>
503 weightedApplyMV (const mv_type& X,
504  mv_type& Y,
505  vector_type& W) const
506 {
507  ConstHostView XView = X.getLocalViewHost(Tpetra::Access::ReadOnly);
508  HostView YView = Y.getLocalViewHost(Tpetra::Access::ReadWrite);
509  ConstHostView WView = W.getLocalViewHost(Tpetra::Access::ReadOnly);
510  weightedApply (XView, YView, WView, 0);
511 }
512 
513 template<class MatrixType, class LocalScalarType>
516 {
517  return "Generic";
518 }
519 
520 template<class MatrixType, class LocalScalarType>
522 solveBlock(ConstHostSubviewLocal X,
523  HostSubviewLocal Y,
524  int blockIndex,
525  Teuchos::ETransp mode,
526  const LSC alpha,
527  const LSC beta) const
528 {
529  TEUCHOS_TEST_FOR_EXCEPT_MSG(true, "Not implemented.");
530 }
531 
532 template<class MatrixType, class LocalScalarType>
533 typename ContainerImpl<MatrixType, LocalScalarType>::LO
536 {
537  LO LO_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
538  GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();
539  const map_type& globalRowMap = *(this->inputMatrix_->getRowMap());
540  const map_type& globalColMap = *(this->inputMatrix_->getColMap());
541  LO rowLID = row;
542  LO dofOffset = 0;
543  if(this->pointIndexed_)
544  {
545  rowLID = row / this->bcrsBlockSize_;
546  dofOffset = row % this->bcrsBlockSize_;
547  }
548  GO diagGID = globalRowMap.getGlobalElement(rowLID);
549  TEUCHOS_TEST_FOR_EXCEPTION(
550  diagGID == GO_INVALID,
551  std::runtime_error, "Ifpack2::Container::translateRowToCol: "
552  "On process " << this->inputMatrix_->getRowMap()->getComm()->getRank() <<
553  ", at least one row index in the set of local "
554  "row indices given to the constructor is not a valid local row index in "
555  "the input matrix's row Map on this process. This should be impossible "
556  "because the constructor checks for this case. Here is the complete set "
557  "of invalid local row indices: " << rowLID << ". "
558  "Please report this bug to the Ifpack2 developers.");
559  //now, can translate diagGID (both a global row AND global col ID) to local column
560  LO colLID = globalColMap.getLocalElement(diagGID);
561  TEUCHOS_TEST_FOR_EXCEPTION(
562  colLID == LO_INVALID,
563  std::runtime_error, "Ifpack2::Container::translateRowToCol: "
564  "On process " << this->inputMatrix_->getRowMap()->getComm()->getRank() << ", "
565  "at least one row index in the set of row indices given to the constructor "
566  "does not have a corresponding column index in the input matrix's column "
567  "Map. This probably means that the column(s) in question is/are empty on "
568  "this process, which would make the submatrix to extract structurally "
569  "singular. The invalid global column index is " << diagGID << ".");
570  //colLID could identify a block column - translate to split column if needed
571  if(this->pointIndexed_)
572  return colLID * this->bcrsBlockSize_ + dofOffset;
573  return colLID;
574 }
575 
576 template<class MatrixType, class LocalScalarType>
578 apply (ConstHostView X,
579  HostView Y,
580  int blockIndex,
581  Teuchos::ETransp mode,
582  SC alpha,
583  SC beta) const
584 {
585  using Teuchos::ArrayView;
586  using Teuchos::as;
587  using Teuchos::RCP;
588  using Teuchos::rcp;
589 
590  // The local operator might have a different Scalar type than
591  // MatrixType. This means that we might have to convert X and Y to
592  // the Tpetra::MultiVector specialization that the local operator
593  // wants. This class' X_ and Y_ internal fields are of the right
594  // type for the local operator, so we can use those as targets.
595 
597 
598  TEUCHOS_TEST_FOR_EXCEPTION(
599  ! this->IsComputed_, std::runtime_error, "Ifpack2::Container::apply: "
600  "You must have called the compute() method before you may call apply(). "
601  "You may call the apply() method as many times as you want after calling "
602  "compute() once, but you must have called compute() at least once.");
603 
604  const size_t numVecs = X.extent(1);
605 
606  if(numVecs == 0) {
607  return; // done! nothing to do
608  }
609 
610  // The local operator works on a permuted subset of the local parts
611  // of X and Y. The subset and permutation are defined by the index
612  // array returned by getBlockRows(). If the permutation is trivial
613  // and the subset is exactly equal to the local indices, then we
614  // could use the local parts of X and Y exactly, without needing to
615  // permute. Otherwise, we have to use temporary storage to permute
616  // X and Y. For now, we always use temporary storage.
617  //
618  // Create temporary permuted versions of the input and output.
619  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
620  // store the permuted versions of X resp. Y. Note that X_local has
621  // the domain Map of the operator, which may be a permuted subset of
622  // the local Map corresponding to X.getMap(). Similarly, Y_local
623  // has the range Map of the operator, which may be a permuted subset
624  // of the local Map corresponding to Y.getMap(). numRows_ here
625  // gives the number of rows in the row Map of the local Inverse_
626  // operator.
627  //
628  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
629  // here that the row Map and the range Map of that operator are
630  // the same.
631  //
632  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
633  // really belongs in Tpetra.
634 
635  if(X_localBlocks_.size() == 0 || X.extent(1) != X_local_.extent(1))
636  {
637  //need to resize (or create for the first time) the three scratch arrays
638  X_localBlocks_.clear();
639  Y_localBlocks_.clear();
640  X_localBlocks_.reserve(this->numBlocks_);
641  Y_localBlocks_.reserve(this->numBlocks_);
642 
643  X_local_ = HostViewLocal("X_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
644  Y_local_ = HostViewLocal("Y_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
645 
646  //create all X_local and Y_local managed Views at once, are
647  //reused in subsequent apply() calls
648  for(int i = 0; i < this->numBlocks_; i++)
649  {
650  auto blockBounds = std::make_pair(this->blockOffsets_[i] * this->scalarsPerRow_,
651  (this->blockOffsets_[i] + this->blockSizes_[i]) * this->scalarsPerRow_);
652  X_localBlocks_.emplace_back(X_local_, blockBounds, Kokkos::ALL());
653  Y_localBlocks_.emplace_back(Y_local_, blockBounds, Kokkos::ALL());
654  }
655  }
656 
657  const ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
658 
659  if(this->scalarsPerRow_ == 1)
660  mvgs.gatherViewToView (X_localBlocks_[blockIndex], X, blockRows);
661  else
662  mvgs.gatherViewToViewBlock (X_localBlocks_[blockIndex], X, blockRows, this->scalarsPerRow_);
663 
664  // We must gather the contents of the output multivector Y even on
665  // input to solveBlock(), since the inverse operator might use it as
666  // an initial guess for a linear solve. We have no way of knowing
667  // whether it does or does not.
668 
669  if(this->scalarsPerRow_ == 1)
670  mvgs.gatherViewToView (Y_localBlocks_[blockIndex], Y, blockRows);
671  else
672  mvgs.gatherViewToViewBlock (Y_localBlocks_[blockIndex], Y, blockRows, this->scalarsPerRow_);
673 
674  // Apply the local operator:
675  // Y_local := beta*Y_local + alpha*M^{-1}*X_local
676  this->solveBlock (X_localBlocks_[blockIndex], Y_localBlocks_[blockIndex], blockIndex, mode,
677  LSC(alpha), LSC(beta));
678 
679  // Scatter the permuted subset output vector Y_local back into the
680  // original output multivector Y.
681  if(this->scalarsPerRow_ == 1)
682  mvgs.scatterViewToView (Y, Y_localBlocks_[blockIndex], blockRows);
683  else
684  mvgs.scatterViewToViewBlock (Y, Y_localBlocks_[blockIndex], blockRows, this->scalarsPerRow_);
685 }
686 
687 template<class MatrixType, class LocalScalarType>
689 weightedApply(ConstHostView X,
690  HostView Y,
691  ConstHostView D,
692  int blockIndex,
693  Teuchos::ETransp mode,
694  SC alpha,
695  SC beta) const
696 {
697  using Teuchos::ArrayRCP;
698  using Teuchos::ArrayView;
699  using Teuchos::Range1D;
700  using Teuchos::Ptr;
701  using Teuchos::ptr;
702  using Teuchos::RCP;
703  using Teuchos::rcp;
704  using Teuchos::rcp_const_cast;
705  using std::endl;
706  using STS = Teuchos::ScalarTraits<SC>;
707 
708  // The local operator template parameter might have a different
709  // Scalar type than MatrixType. This means that we might have to
710  // convert X and Y to the Tpetra::MultiVector specialization that
711  // the local operator wants. This class' X_ and Y_ internal fields
712  // are of the right type for the local operator, so we can use those
713  // as targets.
714 
715  const char prefix[] = "Ifpack2::Container::weightedApply: ";
716  TEUCHOS_TEST_FOR_EXCEPTION(
717  ! this->IsComputed_, std::runtime_error, prefix << "You must have called the "
718  "compute() method before you may call this method. You may call "
719  "weightedApply() as many times as you want after calling compute() once, "
720  "but you must have called compute() at least once first.");
721 
722  //bmk 7-2019: BlockRelaxation already checked this, but if that changes...
723  TEUCHOS_TEST_FOR_EXCEPTION(
724  this->scalarsPerRow_ > 1, std::logic_error, prefix << "Use of block rows isn't allowed "
725  "in overlapping Jacobi (the only method that uses weightedApply");
726 
727  const size_t numVecs = X.extent(1);
728 
729  TEUCHOS_TEST_FOR_EXCEPTION(
730  X.extent(1) != Y.extent(1), std::runtime_error,
731  prefix << "X and Y have different numbers of vectors (columns). X has "
732  << X.extent(1) << ", but Y has " << Y.extent(1) << ".");
733 
734  if(numVecs == 0) {
735  return; // done! nothing to do
736  }
737 
738  const size_t numRows = this->blockSizes_[blockIndex];
739 
740  // The local operator works on a permuted subset of the local parts
741  // of X and Y. The subset and permutation are defined by the index
742  // array returned by getBlockRows(). If the permutation is trivial
743  // and the subset is exactly equal to the local indices, then we
744  // could use the local parts of X and Y exactly, without needing to
745  // permute. Otherwise, we have to use temporary storage to permute
746  // X and Y. For now, we always use temporary storage.
747  //
748  // Create temporary permuted versions of the input and output.
749  // (Re)allocate X_ and/or Y_ only if necessary. We'll use them to
750  // store the permuted versions of X resp. Y. Note that X_local has
751  // the domain Map of the operator, which may be a permuted subset of
752  // the local Map corresponding to X.getMap(). Similarly, Y_local
753  // has the range Map of the operator, which may be a permuted subset
754  // of the local Map corresponding to Y.getMap(). numRows_ here
755  // gives the number of rows in the row Map of the local operator.
756  //
757  // FIXME (mfh 20 Aug 2013) There might be an implicit assumption
758  // here that the row Map and the range Map of that operator are
759  // the same.
760  //
761  // FIXME (mfh 20 Aug 2013) This "local permutation" functionality
762  // really belongs in Tpetra.
763  if(X_localBlocks_.size() == 0 || X.extent(1) != X_local_.extent(1))
764  {
765  //need to resize (or create for the first time) the three scratch arrays
766  X_localBlocks_.clear();
767  Y_localBlocks_.clear();
768  X_localBlocks_.reserve(this->numBlocks_);
769  Y_localBlocks_.reserve(this->numBlocks_);
770 
771  X_local_ = HostViewLocal("X_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
772  Y_local_ = HostViewLocal("Y_local", this->blockRows_.size() * this->scalarsPerRow_, numVecs);
773 
774  //create all X_local and Y_local managed Views at once, are
775  //reused in subsequent apply() calls
776  for(int i = 0; i < this->numBlocks_; i++)
777  {
778  auto blockBounds = std::make_pair(this->blockOffsets_[i] * this->scalarsPerRow_,
779  (this->blockOffsets_[i] + this->blockSizes_[i]) * this->scalarsPerRow_);
780  X_localBlocks_.emplace_back(X_local_, blockBounds, Kokkos::ALL());
781  Y_localBlocks_.emplace_back(Y_local_, blockBounds, Kokkos::ALL());
782  }
783  }
784  if(int(weightedApplyScratch_.extent(0)) != 3 * this->maxBlockSize_ ||
785  weightedApplyScratch_.extent(1) != numVecs)
786  {
787  weightedApplyScratch_ = HostViewLocal("weightedApply scratch", 3 * this->maxBlockSize_, numVecs);
788  }
789 
790  ArrayView<const LO> blockRows = this->getBlockRows(blockIndex);
791 
793 
794  //note: BlockCrs w/ weighted Jacobi isn't allowed, so no need to use block gather/scatter
795  mvgs.gatherViewToView (X_localBlocks_[blockIndex], X, blockRows);
796  // We must gather the output multivector Y even on input to
797  // solveBlock(), since the local operator might use it as an initial
798  // guess for a linear solve. We have no way of knowing whether it
799  // does or does not.
800 
801  mvgs.gatherViewToView (Y_localBlocks_[blockIndex], Y, blockRows);
802 
803  // Apply the diagonal scaling D to the input X. It's our choice
804  // whether the result has the original input Map of X, or the
805  // permuted subset Map of X_local. If the latter, we also need to
806  // gather D into the permuted subset Map. We choose the latter, to
807  // save memory and computation. Thus, we do the following:
808  //
809  // 1. Gather D into a temporary vector D_local.
810  // 2. Create a temporary X_scaled to hold diag(D_local) * X_local.
811  // 3. Compute X_scaled := diag(D_loca) * X_local.
812  auto maxBS = this->maxBlockSize_;
813  auto bs = this->blockSizes_[blockIndex] * this->scalarsPerRow_;
814 
815  HostSubviewLocal D_local(weightedApplyScratch_, std::make_pair(0, bs), std::make_pair(0, 1));
816  mvgs.gatherViewToView (D_local, D, blockRows);
817  HostSubviewLocal X_scaled(weightedApplyScratch_, std::make_pair(maxBS, maxBS + bs), Kokkos::ALL());
818  for(size_t j = 0; j < numVecs; j++)
819  for(size_t i = 0; i < numRows; i++)
820  X_scaled(i, j) = X_localBlocks_[blockIndex](i, j) * D_local(i, 0);
821 
822  HostSubviewLocal Y_temp(weightedApplyScratch_, std::make_pair(maxBS * 2, maxBS * 2 + bs), Kokkos::ALL());
823  // Apply the local operator: Y_temp := M^{-1} * X_scaled
824  this->solveBlock (X_scaled, Y_temp, blockIndex, mode, STS::one(), STS::zero());
825  // Y_local := beta * Y_local + alpha * diag(D_local) * Y_temp.
826  //
827  // Note that we still use the permuted subset scaling D_local here,
828  // because Y_temp has the same permuted subset Map. That's good, in
829  // fact, because it's a subset: less data to read and multiply.
830  LISC a = alpha;
831  LISC b = beta;
832  for(size_t j = 0; j < numVecs; j++)
833  for(size_t i = 0; i < numRows; i++)
834  Y_localBlocks_[blockIndex](i, j) = b * Y_localBlocks_[blockIndex](i, j) + a * Y_temp(i, j) * D_local(i, 0);
835 
836  // Copy the permuted subset output vector Y_local into the original
837  // output multivector Y.
838  mvgs.scatterViewToView (Y, Y_localBlocks_[blockIndex], blockRows);
839 }
840 
841 template<class MatrixType, class LocalScalarType>
843  typename ContainerImpl<MatrixType, LocalScalarType>::SC,
844  typename ContainerImpl<MatrixType, LocalScalarType>::LO,
845  typename ContainerImpl<MatrixType, LocalScalarType>::GO,
846  typename ContainerImpl<MatrixType, LocalScalarType>::NO>
848 getInputRowView(LO row) const
849 {
850 
851  typedef typename MatrixType::nonconst_local_inds_host_view_type nonconst_local_inds_host_view_type;
852  typedef typename MatrixType::nonconst_values_host_view_type nonconst_values_host_view_type;
853 
854  typedef typename MatrixType::local_inds_host_view_type local_inds_host_view_type;
855  typedef typename MatrixType::values_host_view_type values_host_view_type;
856  using IST = typename row_matrix_type::impl_scalar_type;
857 
858  if(this->hasBlockCrs_)
859  {
860  using h_inds_type = typename block_crs_matrix_type::local_inds_host_view_type;
861  using h_vals_type = typename block_crs_matrix_type::values_host_view_type;
862  h_inds_type colinds;
863  h_vals_type values;
864  this->inputBlockMatrix_->getLocalRowView(row / this->bcrsBlockSize_, colinds, values);
865  size_t numEntries = colinds.size();
866  // CMS: Can't say I understand what this really does
867  //return StridedRowView(values + row % this->bcrsBlockSize_, colinds, this->bcrsBlockSize_, numEntries * this->bcrsBlockSize_);
868  h_vals_type subvals = Kokkos::subview(values,std::pair<size_t,size_t>(row % this->bcrsBlockSize_,values.size()));
869  return StridedRowView(subvals, colinds, this->bcrsBlockSize_, numEntries * this->bcrsBlockSize_);
870  }
871  else if(!this->inputMatrix_->supportsRowViews())
872  {
873  size_t maxEntries = this->inputMatrix_->getNodeMaxNumRowEntries();
874  Teuchos::Array<LO> inds(maxEntries);
875  Teuchos::Array<SC> vals(maxEntries);
876  nonconst_local_inds_host_view_type inds_v(inds.data(),maxEntries);
877  nonconst_values_host_view_type vals_v(reinterpret_cast<IST*>(vals.data()),maxEntries);
878  size_t numEntries;
879  this->inputMatrix_->getLocalRowCopy(row, inds_v, vals_v, numEntries);
880  vals.resize(numEntries); inds.resize(numEntries);
881  return StridedRowView(vals, inds);
882  }
883  else
884  {
885  // CMS - This is dangerous and might not work.
886  local_inds_host_view_type colinds;
887  values_host_view_type values;
888  this->inputMatrix_->getLocalRowView(row, colinds, values);
889  return StridedRowView(values, colinds, 1, colinds.size());
890  }
891 }
892 
893 template<class MatrixType, class LocalScalarType>
895 clearBlocks()
896 {
897  X_localBlocks_.clear();
898  Y_localBlocks_.clear();
899  X_local_ = HostViewLocal();
900  Y_local_ = HostViewLocal();
902 }
903 
904 namespace Details {
905 
906 //Implementation of Ifpack2::Details::StridedRowView
907 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
908 StridedRowView<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
909 StridedRowView(h_vals_type vals_, h_inds_type inds_, int blockSize_, size_t nnz_)
910  : vals(vals_), inds(inds_), blockSize(blockSize_), nnz(nnz_)
911 {}
912 
913 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
915 StridedRowView(Teuchos::Array<SC>& vals_, Teuchos::Array<LO>& inds_)
916  : vals(), inds(), blockSize(1), nnz(vals_.size())
917 {
918  valsCopy.swap(vals_);
919  indsCopy.swap(inds_);
920 }
921 
922 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
924 val(size_t i) const
925 {
926  #ifdef HAVE_IFPACK2_DEBUG
927  TEUCHOS_TEST_FOR_EXCEPTION(i >= nnz, std::runtime_error,
928  "Out-of-bounds access into Ifpack2::Container::StridedRowView");
929  #endif
930  if(vals.size() > 0)
931  {
932  if(blockSize == 1)
933  return vals[i];
934  else
935  return vals[i * blockSize];
936  }
937  else
938  return valsCopy[i];
939 }
940 
941 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
942 LocalOrdinal StridedRowView<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
943 ind(size_t i) const
944 {
945  #ifdef HAVE_IFPACK2_DEBUG
946  TEUCHOS_TEST_FOR_EXCEPTION(i >= nnz, std::runtime_error,
947  "Out-of-bounds access into Ifpack2::Container::StridedRowView");
948  #endif
949  //inds is smaller than vals by a factor of the block size (dofs/node)
950  if(inds.size() > 0)
951  {
952  if(blockSize == 1)
953  return inds[i];
954  else
955  return inds[i / blockSize] * blockSize + i % blockSize;
956  }
957  else
958  return indsCopy[i];
959 }
960 
961 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
962 size_t StridedRowView<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
963 size() const
964 {
965  return nnz;
966 }
967 }
968 
969 }
970 
971 template <class MatrixType>
972 std::ostream& operator<<(std::ostream& os, const Ifpack2::Container<MatrixType>& obj)
973 {
974  return obj.print(os);
975 }
976 
977 #define IFPACK2_CONTAINER_INSTANT(S,LO,GO,N) \
978  template class Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N>>; \
979  template class Ifpack2::ContainerImpl<Tpetra::RowMatrix<S, LO, GO, N>, S>; \
980  template class Ifpack2::Details::StridedRowView<S, LO, GO, N>; \
981  template std::ostream& operator<< <Tpetra::RowMatrix<S, LO, GO, N>>( \
982  std::ostream& os, const Ifpack2::Container<Tpetra::RowMatrix<S, LO, GO, N>>& obj);
983 
984 #endif
985 
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:296
StridedRowView(h_vals_type vals_, h_inds_type inds_, int blockSize_, size_t nnz_)
Constructor for row views (preferred)
Definition: Ifpack2_Container_def.hpp:909
bool pointIndexed_
(If hasBlockCrs_) Whether the blocks are described using sub-block row indices instead of full block ...
Definition: Ifpack2_Container_decl.hpp:314
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container_decl.hpp:292
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
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:343
LO scalarsPerRow_
Definition: Ifpack2_Container_decl.hpp:317
virtual ~Container()
Destructor.
Definition: Ifpack2_Container_def.hpp:114
bool hasBlockCrs_
Whether the input matrix is a BlockCRS matrix.
Definition: Ifpack2_Container_decl.hpp:310
InverseType::scalar_type LSC
The internal representation of LocalScalarType in Kokkos::View.
Definition: Ifpack2_Container_decl.hpp:363
Ifpack2 implementation details.
bool isInitialized() const
Whether the container has been successfully initialized.
Definition: Ifpack2_Container_def.hpp:165
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:312
Teuchos::RCP< const block_crs_matrix_type > inputBlockMatrix_
The input matrix, dynamic cast to BlockCrsMatrix. May be null.
Definition: Ifpack2_Container_decl.hpp:289
Teuchos::RCP< const row_matrix_type > inputMatrix_
The input matrix to the constructor.
Definition: Ifpack2_Container_decl.hpp:283
typename mv_type::dual_view_type::t_host HostView
Definition: Ifpack2_Container_decl.hpp:139
void DoGSBlock(ConstHostView 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
bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_Container_def.hpp:170
bool IsParallel_
Whether the problem is distributed across multiple MPI processes.
Definition: Ifpack2_Container_decl.hpp:302
virtual void weightedApplyMV(const mv_type &X, mv_type &Y, vector_type &W) const
Wrapper for weightedApply with MultiVector.
Definition: Ifpack2_Container_def.hpp:183
Definition: Ifpack2_Container_decl.hpp:576
Interface for creating and solving a set of local linear problems.
Definition: Ifpack2_Container_decl.hpp:112
virtual void applyMV(const mv_type &X, mv_type &Y) const
Wrapper for apply with MultiVector.
Definition: Ifpack2_Container_def.hpp:176
Structure for read-only views of general matrix rows.
Definition: Ifpack2_Container_decl.hpp:330
virtual void DoGSBlock(ConstHostView 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
static std::string getName()
Definition: Ifpack2_Container_def.hpp:190
GO NumGlobalRows_
Number of global rows in input matrix.
Definition: Ifpack2_Container_decl.hpp:306
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
GO NumGlobalNonzeros_
Number of nonzeros in input matrix.
Definition: Ifpack2_Container_decl.hpp:308
LO NumLocalRows_
Number of local rows in input matrix.
Definition: Ifpack2_Container_decl.hpp:304
Implementation detail of Ifpack2::Container subclasses.
Definition: Ifpack2_Details_MultiVectorLocalGatherScatter.hpp:85