Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Container.hpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_CONTAINER_HPP
44 #define IFPACK2_CONTAINER_HPP
45 
48 
49 #include "Ifpack2_ConfigDefs.hpp"
50 #include "Tpetra_RowMatrix.hpp"
51 #include "Teuchos_Describable.hpp"
52 #include <Ifpack2_Partitioner.hpp>
53 #include <Tpetra_Map.hpp>
54 #include <Tpetra_Experimental_BlockCrsMatrix_decl.hpp>
56 #include <Teuchos_Time.hpp>
57 #include <iostream>
58 #include <type_traits>
59 
60 #ifndef DOXYGEN_SHOULD_SKIP_THIS
61 namespace Teuchos {
62  // Forward declaration to avoid include.
63  class ParameterList;
64 } // namespace Teuchos
65 #endif // DOXYGEN_SHOULD_SKIP_THIS
66 
67 namespace Ifpack2 {
68 
113 template<class MatrixType>
116 
117 protected:
118  typedef typename MatrixType::scalar_type scalar_type;
119  typedef typename MatrixType::local_ordinal_type local_ordinal_type;
120  typedef typename MatrixType::global_ordinal_type global_ordinal_type;
121  typedef typename MatrixType::node_type node_type;
122  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> mv_type;
123  typedef Tpetra::Vector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> vector_type;
124  typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
126  typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type;
128  typedef Tpetra::Experimental::BlockCrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> block_crs_matrix_type;
129  typedef Tpetra::RowMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type> row_matrix_type;
130 
131  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
132  "Ifpack2::Container: Please use MatrixType = Tpetra::RowMatrix.");
133 
135  typedef typename Kokkos::Details::ArithTraits<scalar_type>::val_type impl_scalar_type;
137 
138 public:
139  typedef typename mv_type::dual_view_type::t_host HostView;
140 
151  const Teuchos::RCP<const import_type>& importer,
152  int OverlapLevel,
153  scalar_type DampingFactor) :
154  inputMatrix_ (matrix),
155  OverlapLevel_ (OverlapLevel),
156  DampingFactor_ (DampingFactor),
157  Importer_ (importer)
158  {
159  using Teuchos::Ptr;
160  using Teuchos::RCP;
161  using Teuchos::Array;
162  using Teuchos::ArrayView;
163  using Teuchos::Comm;
164  // typedef Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type> import_type; // unused
165  NumLocalRows_ = inputMatrix_->getNodeNumRows();
166  NumGlobalRows_ = inputMatrix_->getGlobalNumRows();
167  NumGlobalNonzeros_ = inputMatrix_->getGlobalNumEntries();
168  IsParallel_ = inputMatrix_->getRowMap()->getComm()->getSize() != 1;
169  auto global_bcrs =
170  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(inputMatrix_);
171  hasBlockCrs_ = !global_bcrs.is_null();
172  if(hasBlockCrs_)
173  bcrsBlockSize_ = global_bcrs->getBlockSize();
174  else
175  bcrsBlockSize_ = 1;
176  setBlockSizes(partitions);
177  }
178 
191  const Teuchos::Array<local_ordinal_type>& localRows) :
192  inputMatrix_ (matrix),
193  numBlocks_ (1),
194  partitions_ (localRows.size()),
195  blockRows_ (1),
196  partitionIndices_ (1),
197  OverlapLevel_ (0),
198  DampingFactor_ (STS::one()),
199  Importer_ (Teuchos::null)
200  {
201  NumLocalRows_ = inputMatrix_->getNodeNumRows();
202  NumGlobalRows_ = inputMatrix_->getGlobalNumRows();
203  NumGlobalNonzeros_ = inputMatrix_->getGlobalNumEntries();
204  IsParallel_ = inputMatrix_->getRowMap()->getComm()->getSize() > 1;
205  blockRows_[0] = localRows.size();
206  partitionIndices_[0] = 0;
207  const block_crs_matrix_type* global_bcrs =
208  Teuchos::rcp_dynamic_cast<const block_crs_matrix_type>(inputMatrix_).get();
209  hasBlockCrs_ = global_bcrs;
210  if(hasBlockCrs_)
211  bcrsBlockSize_ = global_bcrs->getBlockSize();
212  else
213  bcrsBlockSize_ = 1;
214  for(local_ordinal_type i = 0; i < localRows.size(); i++)
215  partitions_[i] = localRows[i];
216  }
217 
219  virtual ~Container() {};
220 
237  {
239  (&partitions_[partitionIndices_[blockIndex]], blockRows_[blockIndex]);
240  }
241 
250  virtual void initialize () = 0;
251 
254  {
255  //First, create a grand total of all the rows in all the blocks
256  //Note: If partitioner allowed overlap, this could be greater than the # of local rows
257  local_ordinal_type totalBlockRows = 0;
258  numBlocks_ = partitions.size();
261  for(int i = 0; i < numBlocks_; i++)
262  {
263  local_ordinal_type rowsInBlock = partitions[i].size();
264  blockRows_[i] = rowsInBlock;
265  partitionIndices_[i] = totalBlockRows;
266  totalBlockRows += rowsInBlock;
267  }
268  partitions_.resize(totalBlockRows);
269  //set partitions_: each entry is the partition/block of the row
270  local_ordinal_type iter = 0;
271  for(int i = 0; i < numBlocks_; i++)
272  {
273  for(int j = 0; j < blockRows_[i]; j++)
274  {
275  partitions_[iter++] = partitions[i][j];
276  }
277  }
278  }
279 
280  void getMatDiag() const
281  {
282  if(Diag_.is_null())
283  {
284  Diag_ = rcp(new vector_type(inputMatrix_->getDomainMap()));
285  inputMatrix_->getLocalDiagCopy(*Diag_);
286  }
287  }
288 
294  virtual void compute () = 0;
295 
296  void DoJacobi(HostView& X, HostView& Y, int stride) const;
297  void DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W, int stride) const;
298  void DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2, int stride) const;
299  void DoSGS(HostView& X, HostView& Y, HostView& Y2, int stride) const;
300 
302  virtual void setParameters (const Teuchos::ParameterList& List) = 0;
303 
305  virtual bool isInitialized () const = 0;
306 
308  virtual bool isComputed () const = 0;
309 
322  virtual void
323  apply (HostView& X,
324  HostView& Y,
325  int blockIndex,
326  int stride,
327  Teuchos::ETransp mode = Teuchos::NO_TRANS,
328  scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
329  scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const = 0;
330 
332  void applyMV (mv_type& X, mv_type& Y) const
333  {
334  HostView XView = X.template getLocalView<Kokkos::HostSpace>();
335  HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
336  this->apply (XView, YView, 0, X.getStride());
337  }
338 
359  virtual void
360  weightedApply (HostView& X,
361  HostView& Y,
362  HostView& W,
363  int blockIndex,
364  int stride,
365  Teuchos::ETransp mode = Teuchos::NO_TRANS,
366  scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one(),
367  scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero()) const = 0;
368 
370  void weightedApplyMV (mv_type& X,
371  mv_type& Y,
372  vector_type& W)
373  {
374  HostView XView = X.template getLocalView<Kokkos::HostSpace>();
375  HostView YView = Y.template getLocalView<Kokkos::HostSpace>();
376  HostView WView = W.template getLocalView<Kokkos::HostSpace>();
377  weightedApply (XView, YView, WView, 0, X.getStride());
378  }
379 
380  virtual void clearBlocks();
381 
383  virtual std::ostream& print (std::ostream& os) const = 0;
384 
387  static std::string getName()
388  {
389  return "Generic";
390  }
391 
392 protected:
395 
399  Teuchos::Array<local_ordinal_type> partitions_; //size: total # of local rows (in all local blocks)
411  scalar_type DampingFactor_;
415  local_ordinal_type NumLocalRows_;
417  global_ordinal_type NumGlobalRows_;
419  global_ordinal_type NumGlobalNonzeros_;
424 };
425 
427 template <class MatrixType>
428 inline std::ostream&
429 operator<< (std::ostream& os, const Ifpack2::Container<MatrixType>& obj)
430 {
431  return obj.print (os);
432 }
433 
434 template <class MatrixType>
435 void Container<MatrixType>::DoJacobi(HostView& X, HostView& Y, int stride) const
436 {
437  const scalar_type one = STS::one();
438  // Note: Flop counts copied naively from Ifpack.
439  // use partitions_ and blockRows_
440  size_t numVecs = X.dimension_1();
441  // Non-overlapping Jacobi
442  for (local_ordinal_type i = 0; i < numBlocks_; i++)
443  {
444  // may happen that a partition is empty
445  if(blockRows_[i] != 1 || hasBlockCrs_)
446  {
447  if(blockRows_[i] == 0 )
448  continue;
449  apply(X, Y, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
450  }
451  else // singleton, can't access Containers_[i] as it was never filled and may be null.
452  {
453  local_ordinal_type LRID = partitions_[partitionIndices_[i]];
454  getMatDiag();
455  HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
456  impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
457  for(size_t nv = 0; nv < numVecs; nv++)
458  {
459  impl_scalar_type x = X(LRID, nv);
460  Y(LRID, nv) = x * d;
461  }
462  }
463  }
464 }
465 
466 template <class MatrixType>
467 void Container<MatrixType>::DoOverlappingJacobi(HostView& X, HostView& Y, HostView& W, int stride) const
468 {
469  // Overlapping Jacobi
470  for(local_ordinal_type i = 0; i < numBlocks_; i++)
471  {
472  // may happen that a partition is empty
473  if(blockRows_[i] == 0)
474  continue;
475  if(blockRows_[i] != 1)
476  weightedApply(X, Y, W, i, stride, Teuchos::NO_TRANS, DampingFactor_, STS::one());
477  }
478 }
479 
480 template<class MatrixType>
481 void Container<MatrixType>::DoGaussSeidel(HostView& X, HostView& Y, HostView& Y2, int stride) const
482 {
483  using Teuchos::Array;
484  using Teuchos::ArrayRCP;
485  using Teuchos::ArrayView;
486  using Teuchos::Ptr;
487  using Teuchos::RCP;
488  using Teuchos::rcp;
489  using Teuchos::rcpFromRef;
490  // Note: Flop counts copied naively from Ifpack.
491  const scalar_type one = STS::one();
492  const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
493  auto numVecs = X.dimension_1();
494  Array<scalar_type> Values;
495  Array<local_ordinal_type> Indices;
496  Indices.resize(Length);
497 
498  Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length); //note: if A was not a BlockCRS, bcrsBlockSize_ is 1
499  // I think I decided I need two extra vectors:
500  // One to store the sum of the corrections (initialized to zero)
501  // One to store the temporary residual (doesn't matter if it is zeroed or not)
502  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
503  HostView Resid("", X.dimension_0(), X.dimension_1());
504  for(local_ordinal_type i = 0; i < numBlocks_; i++)
505  {
506  if(blockRows_[i] > 1 || hasBlockCrs_)
507  {
508  if (blockRows_[i] == 0)
509  continue;
510  // update from previous block
511  ArrayView<const local_ordinal_type> localRows = getLocalRows (i);
512  const size_t localNumRows = blockRows_[i];
513  for(size_t j = 0; j < localNumRows; j++)
514  {
515  const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
516  size_t NumEntries;
517  inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
518  for(size_t m = 0; m < numVecs; m++)
519  {
520  if(hasBlockCrs_)
521  {
522  for (int localR = 0; localR < bcrsBlockSize_; localR++)
523  Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
524  for (size_t k = 0; k < NumEntries; ++k)
525  {
526  const local_ordinal_type col = Indices[k];
527  for (int localR = 0; localR < bcrsBlockSize_; localR++)
528  {
529  for(int localC = 0; localC < bcrsBlockSize_; localC++)
530  {
531  Resid(LID * bcrsBlockSize_ + localR, m) -=
532  Values[k * bcrsBlockSize_ * bcrsBlockSize_ + localR + localC * bcrsBlockSize_]
533  * Y2(col * bcrsBlockSize_ + localC, m);
534  }
535  }
536  }
537  }
538  else
539  {
540  Resid(LID, m) = X(LID, m);
541  for (size_t k = 0; k < NumEntries; ++k)
542  {
543  const local_ordinal_type col = Indices[k];
544  Resid(LID, m) -= Values[k] * Y2(col, m);
545  }
546  }
547  }
548  }
549  // solve with this block
550  //
551  // Note: I'm abusing the ordering information, knowing that X/Y
552  // and Y2 have the same ordering for on-proc unknowns.
553  //
554  // Note: Add flop counts for inverse apply
555  apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
556  }
557  else if(blockRows_[i] == 1)
558  {
559  // singleton, can't access Containers_[i] as it was never filled and may be null.
560  // a singleton calculation is exact, all residuals should be zero.
561  local_ordinal_type LRID = partitionIndices_[i]; // by definition, a singleton 1 row in block.
562  getMatDiag();
563  HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
564  impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
565  for(size_t m = 0; m < numVecs; m++)
566  {
567  impl_scalar_type x = X(LRID, m);
568  impl_scalar_type newy = x * d;
569  Y2(LRID, m) = newy;
570  }
571  } // end else
572  } // end for numBlocks_
573  if(IsParallel_)
574  {
575  auto numMyRows = inputMatrix_->getNodeNumRows();
576  for (size_t m = 0; m < numVecs; ++m)
577  {
578  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
579  {
580  Y(i, m) = Y2(i, m);
581  }
582  }
583  }
584 }
585 
586 template<class MatrixType>
587 void Container<MatrixType>::DoSGS(HostView& X, HostView& Y, HostView& Y2, int stride) const
588 {
589  using Teuchos::Array;
590  using Teuchos::ArrayRCP;
591  using Teuchos::ArrayView;
592  using Teuchos::Ptr;
593  using Teuchos::ptr;
594  using Teuchos::RCP;
595  using Teuchos::rcp;
596  using Teuchos::rcpFromRef;
597  const scalar_type one = STS::one();
598  auto numVecs = X.dimension_1();
599  const size_t Length = inputMatrix_->getNodeMaxNumRowEntries();
600  Array<scalar_type> Values;
601  Array<local_ordinal_type> Indices(Length);
602  Values.resize(bcrsBlockSize_ * bcrsBlockSize_ * Length);
603  // I think I decided I need two extra vectors:
604  // One to store the sum of the corrections (initialized to zero)
605  // One to store the temporary residual (doesn't matter if it is zeroed or not)
606  // My apologies for making the names clear and meaningful. (X=RHS, Y=guess?! Nice.)
607  HostView Resid("", X.dimension_0(), X.dimension_1());
608  // Forward Sweep
609  for(local_ordinal_type i = 0; i < numBlocks_; i++)
610  {
611  if(blockRows_[i] != 1 || hasBlockCrs_)
612  {
613  if(blockRows_[i] == 0)
614  continue; // Skip empty partitions
615  // update from previous block
616  ArrayView<const local_ordinal_type> localRows = getLocalRows(i);
617  for(local_ordinal_type j = 0; j < blockRows_[i]; j++)
618  {
619  const local_ordinal_type LID = localRows[j]; // Containers_[i]->ID (j);
620  size_t NumEntries;
621  inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
622  //set tmpresid = initresid - A*correction
623  for(size_t m = 0; m < numVecs; m++)
624  {
625  if(hasBlockCrs_)
626  {
627  for(int localR = 0; localR < bcrsBlockSize_; localR++)
628  Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
629  for(size_t k = 0; k < NumEntries; ++k)
630  {
631  const local_ordinal_type col = Indices[k];
632  for (int localR = 0; localR < bcrsBlockSize_; localR++)
633  {
634  for(int localC = 0; localC < bcrsBlockSize_; localC++)
635  Resid(LID * bcrsBlockSize_ + localR, m) -=
636  Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
637  * Y2(col * bcrsBlockSize_ + localC, m);
638  }
639  }
640  }
641  else
642  {
643  Resid(LID, m) = X(LID, m);
644  for(size_t k = 0; k < NumEntries; k++)
645  {
646  local_ordinal_type col = Indices[k];
647  Resid(LID, m) -= Values[k] * Y2(col, m);
648  }
649  }
650  }
651  }
652  // solve with this block
653  //
654  // Note: I'm abusing the ordering information, knowing that X/Y
655  // and Y2 have the same ordering for on-proc unknowns.
656  //
657  // Note: Add flop counts for inverse apply
658  apply(Resid, Y2, i, stride, Teuchos::NO_TRANS, DampingFactor_, one);
659  // operations for all getrow's
660  }
661  else // singleton, can't access Containers_[i] as it was never filled and may be null.
662  {
663  local_ordinal_type LRID = partitions_[partitionIndices_[i]];
664  getMatDiag();
665  HostView diagView = Diag_->template getLocalView<Kokkos::HostSpace>();
666  impl_scalar_type d = (impl_scalar_type) one / diagView(LRID, 0);
667  for(size_t m = 0; m < numVecs; m++)
668  {
669  impl_scalar_type x = X(LRID, m);
670  impl_scalar_type newy = x * d;
671  Y2(LRID, m) = newy;
672  }
673  } // end else
674  } // end forward sweep over NumLocalBlocks
675  // Reverse Sweep
676  //
677  // mfh 12 July 2013: The unusual iteration bounds, and the use of
678  // i-1 rather than i in the loop body, ensure correctness even if
679  // local_ordinal_type is unsigned. "i = numBlocks_-1; i >= 0;
680  // i--" will loop forever if local_ordinal_type is unsigned, because
681  // unsigned integers are (trivially) always nonnegative.
682  for(local_ordinal_type i = numBlocks_; i > 0; --i)
683  {
684  if(hasBlockCrs_ || blockRows_[i-1] != 1)
685  {
686  if(blockRows_[i - 1] == 0)
687  continue;
688  // update from previous block
689  ArrayView<const local_ordinal_type> localRows = getLocalRows(i-1);
690  for(local_ordinal_type j = 0; j < blockRows_[i-1]; j++)
691  {
692  const local_ordinal_type LID = localRows[j]; // Containers_[i-1]->ID (j);
693  size_t NumEntries;
694  inputMatrix_->getLocalRowCopy(LID, Indices(), Values(), NumEntries);
695  //set tmpresid = initresid - A*correction
696  for (size_t m = 0; m < numVecs; m++)
697  {
698  if(hasBlockCrs_)
699  {
700  for(int localR = 0; localR < bcrsBlockSize_; localR++)
701  Resid(LID * bcrsBlockSize_ + localR, m) = X(LID * bcrsBlockSize_ + localR, m);
702  for(size_t k = 0; k < NumEntries; ++k)
703  {
704  const local_ordinal_type col = Indices[k];
705  for(int localR = 0; localR < bcrsBlockSize_; localR++)
706  {
707  for(int localC = 0; localC < bcrsBlockSize_; localC++)
708  Resid(LID*bcrsBlockSize_+localR, m) -=
709  Values[k * (bcrsBlockSize_ * bcrsBlockSize_) + (localR + localC * bcrsBlockSize_)]
710  * Y2(col * bcrsBlockSize_ + localC, m);
711  }
712  }
713  }
714  else
715  {
716  Resid(LID, m) = X(LID, m);
717  for(size_t k = 0; k < NumEntries; ++k)
718  {
719  local_ordinal_type col = Indices[k];
720  Resid(LID, m) -= Values[k] * Y2(col, m);
721  }
722  }
723  }
724  }
725  // solve with this block
726  //
727  // Note: I'm abusing the ordering information, knowing that X/Y
728  // and Y2 have the same ordering for on-proc unknowns.
729  //
730  // Note: Add flop counts for inverse apply
731  apply(Resid, Y2, i - 1, stride, Teuchos::NO_TRANS, DampingFactor_, one);
732  // operations for all getrow's
733  } // end Partitioner_->numRowsInPart (i) != 1 )
734  // else do nothing, as by definition with a singleton, the residuals are zero.
735  } //end reverse sweep
736  if(IsParallel_)
737  {
738  auto numMyRows = inputMatrix_->getNodeNumRows();
739  for (size_t m = 0; m < numVecs; ++m)
740  {
741  for (size_t i = 0; i < numMyRows * bcrsBlockSize_; ++i)
742  {
743  Y(i, m) = Y2(i, m);
744  }
745  }
746  }
747 }
748 
749 template<class MatrixType>
750 void Container<MatrixType>::clearBlocks()
751 {
752  numBlocks_ = 0;
753  partitions_.clear();
754  blockRows_.clear();
755  partitionIndices_.clear();
756  Diag_ = Teuchos::null; //Diag_ will be recreated if needed
757 }
758 
759 } // namespace Ifpack2
760 
761 namespace Teuchos {
762 
767 template<class MatrixType>
768 class TEUCHOSCORE_LIB_DLL_EXPORT TypeNameTraits< ::Ifpack2::Container<MatrixType> >
769 {
770  public:
771  static std::string name () {
772  return std::string ("Ifpack2::Container<") +
774  }
775 
776  static std::string concreteName (const ::Ifpack2::Container<MatrixType>&) {
777  return name ();
778  }
779 };
780 
781 } // namespace Teuchos
782 
783 #endif // IFPACK2_CONTAINER_HPP
int OverlapLevel_
Number of rows of overlap for adjacent blocks.
Definition: Ifpack2_Container.hpp:409
global_ordinal_type NumGlobalNonzeros_
Number of nonzeros in input matrix.
Definition: Ifpack2_Container.hpp:419
virtual void initialize()=0
Do all set-up operations that only require matrix structure.
int numBlocks_
The number of blocks (partitions) in the container.
Definition: Ifpack2_Container.hpp:397
Teuchos::Array< local_ordinal_type > blockRows_
Number of rows in each block.
Definition: Ifpack2_Container.hpp:401
static std::string getName()
Definition: Ifpack2_Container.hpp:387
global_ordinal_type NumGlobalRows_
Number of global rows in input matrix.
Definition: Ifpack2_Container.hpp:417
bool hasBlockCrs_
Whether the input matrix is a BlockCRS matrix.
Definition: Ifpack2_Container.hpp:421
scalar_type DampingFactor_
Damping factor, passed to apply() as alpha.
Definition: Ifpack2_Container.hpp:411
bool is_null() const
Ifpack2::Partitioner:
Definition: Ifpack2_Partitioner.hpp:179
int bcrsBlockSize_
If hasBlockCrs_, the number of DOFs per vertex. Otherwise 1.
Definition: Ifpack2_Container.hpp:423
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual bool isComputed() const =0
Return true if the container has been successfully computed.
local_ordinal_type NumLocalRows_
Number of local rows in input matrix.
Definition: Ifpack2_Container.hpp:415
virtual void weightedApply(HostView &X, HostView &Y, HostView &W, int blockIndex, int stride, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const =0
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
Teuchos::RCP< const row_matrix_type > inputMatrix_
The input matrix to the constructor.
Definition: Ifpack2_Container.hpp:394
size_type size() const
Teuchos::RCP< vector_type > Diag_
Diagonal elements.
Definition: Ifpack2_Container.hpp:405
void resize(size_type new_size, const value_type &x=value_type())
Container(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions, const Teuchos::RCP< const import_type > &importer, int OverlapLevel, scalar_type DampingFactor)
Constructor.
Definition: Ifpack2_Container.hpp:149
Teuchos::RCP< const Tpetra::Import< local_ordinal_type, global_ordinal_type, node_type > > Importer_
Importer for importing off-process elements of MultiVectors.
Definition: Ifpack2_Container.hpp:413
bool IsParallel_
Whether the problem is distributed across multiple MPI processes.
Definition: Ifpack2_Container.hpp:407
void weightedApplyMV(mv_type &X, mv_type &Y, vector_type &W)
Wrapper for weightedApply with MVs, used in unit tests (never called by BlockRelaxation) ...
Definition: Ifpack2_Container.hpp:370
Teuchos::ArrayView< const local_ordinal_type > getLocalRows(int blockIndex) const
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container.hpp:236
Teuchos::Array< local_ordinal_type > partitionIndices_
Starting index in partitions_ of local row indices for each block.
Definition: Ifpack2_Container.hpp:403
virtual void compute()=0
Extract the local diagonal block and prepare the solver.
Interface for creating and solving a local linear problem.
Definition: Ifpack2_Container.hpp:114
virtual void apply(HostView &X, HostView &Y, int blockIndex, int stride, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const =0
Compute Y := alpha * M^{-1} X + beta*Y.
virtual void setParameters(const Teuchos::ParameterList &List)=0
Set parameters.
virtual std::ostream & print(std::ostream &os) const =0
Print basic information about the container to os.
virtual ~Container()
Destructor.
Definition: Ifpack2_Container.hpp:219
Teuchos::Array< local_ordinal_type > partitions_
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container.hpp:399
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
virtual bool isInitialized() const =0
Return true if the container has been successfully initialized.
void setBlockSizes(const Teuchos::Array< Teuchos::Array< local_ordinal_type > > &partitions)
Initialize arrays with information about block sizes.
Definition: Ifpack2_Container.hpp:253
Kokkos::Details::ArithTraits< scalar_type >::val_type impl_scalar_type
Internal representation of Scalar in Kokkos::View.
Definition: Ifpack2_Container.hpp:132
Container(const Teuchos::RCP< const row_matrix_type > &matrix, const Teuchos::Array< local_ordinal_type > &localRows)
Constructor for single block.
Definition: Ifpack2_Container.hpp:190
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.hpp:332
static std::string name()