Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Container_decl.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_DECL_HPP
44 #define IFPACK2_CONTAINER_DECL_HPP
45 
48 
49 #include "Ifpack2_ConfigDefs.hpp"
50 #include "Tpetra_RowMatrix.hpp"
51 #include "Teuchos_Describable.hpp"
52 #include <Tpetra_Map.hpp>
53 #include <Tpetra_BlockCrsMatrix.hpp>
55 #include <iostream>
56 
57 #ifndef DOXYGEN_SHOULD_SKIP_THIS
58 namespace Teuchos {
59  // Forward declaration to avoid include.
60  class ParameterList;
61 } // namespace Teuchos
62 #endif // DOXYGEN_SHOULD_SKIP_THIS
63 
64 namespace Ifpack2 {
65 
111 template<class MatrixType>
113 {
114 protected:
115  using scalar_type = typename MatrixType::scalar_type;
116  using local_ordinal_type = typename MatrixType::local_ordinal_type;
117  using global_ordinal_type = typename MatrixType::global_ordinal_type;
118  using node_type = typename MatrixType::node_type;
119  using SC = scalar_type;
120  using LO = local_ordinal_type;
121  using GO = global_ordinal_type;
122  using NO = node_type;
123  using import_type = Tpetra::Import<LO, GO, NO>;
124  using mv_type = Tpetra::MultiVector<SC, LO, GO, NO>;
125  using vector_type = Tpetra::Vector<SC, LO, GO, NO>;
126  using map_type = Tpetra::Map<LO, GO, NO>;
127  using crs_matrix_type = Tpetra::CrsMatrix<SC, LO, GO, NO>;
128  using block_crs_matrix_type = Tpetra::BlockCrsMatrix<SC, LO, GO, NO>;
129  using row_matrix_type = Tpetra::RowMatrix<SC, LO, GO, NO>;
130 
131  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
132  "Ifpack2::Container: Please use MatrixType = Tpetra::RowMatrix.");
133 
135  using ISC = typename Kokkos::Details::ArithTraits<SC>::val_type;
136 
139  using HostView = typename mv_type::dual_view_type::t_host;
140 
141 public:
153  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
154  bool pointIndexed);
155 
157  virtual ~Container();
158 
174  Teuchos::ArrayView<const LO> getBlockRows(int blockIndex) const;
175 
184  virtual void initialize () = 0;
185 
187  void setBlockSizes(const Teuchos::Array<Teuchos::Array<LO> >& partitions);
188 
189  void getMatDiag() const;
190 
192  bool isInitialized () const;
193 
195  bool isComputed () const;
196 
205  virtual void compute () = 0;
206 
207  void DoJacobi(HostView X, HostView Y, SC dampingFactor) const;
208  void DoOverlappingJacobi(HostView X, HostView Y, HostView W, SC dampingFactor) const;
209  void DoGaussSeidel(HostView X, HostView Y, HostView Y2, SC dampingFactor) const;
210  void DoSGS(HostView X, HostView Y, HostView Y2, SC dampingFactor) const;
211 
213  virtual void setParameters (const Teuchos::ParameterList& List) = 0;
214 
227  virtual void
228  apply(HostView X,
229  HostView Y,
230  int blockIndex,
232  SC alpha = Teuchos::ScalarTraits<SC>::one(),
233  SC beta = Teuchos::ScalarTraits<SC>::zero()) const = 0;
234 
236  virtual void
238  HostView Y,
239  HostView D,
240  int blockIndex,
242  SC alpha = Teuchos::ScalarTraits<SC>::one(),
243  SC beta = Teuchos::ScalarTraits<SC>::zero()) const = 0;
244 
247  // The underlying container implements the splitting <tt>A = D + R</tt>. Only
248  // it can have efficient access to D and R, as these are constructed in the
249  // symbolic and numeric phases.
250  //
251  // This is the first performance-portable implementation of a block
252  // relaxation, and it is supported currently only by BlockTriDiContainer.
253  virtual void applyInverseJacobi (const mv_type& /* X */, mv_type& /* Y */,
254  SC dampingFactor,
255  bool /* zeroStartingSolution = false */,
256  int /* numSweeps = 1 */) const = 0;
257 
259  virtual void applyMV (mv_type& X, mv_type& Y) const;
260 
262  virtual void weightedApplyMV (mv_type& X,
263  mv_type& Y,
264  vector_type& W) const;
265 
266  virtual void clearBlocks();
267 
269  virtual std::ostream& print (std::ostream& os) const = 0;
270 
273  static std::string getName();
274 
275 protected:
276 
278  virtual void DoGSBlock(HostView X, HostView Y, HostView Y2, HostView Resid,
279  SC dampingFactor, LO i) const;
280 
283 
286 
289 
293  Teuchos::Array<LO> blockRows_; //size: total # of local rows (in all local blocks)
295  Teuchos::Array<LO> blockSizes_; //size: # of blocks
317 
322 
323  LO maxBlockSize_;
324 };
325 
326 namespace Details
327 {
328  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
330 }
331 
340 
341 template<class MatrixType, class LocalScalarType>
342 class ContainerImpl : public Container<MatrixType>
343 {
345 
346 protected:
347  using local_scalar_type = LocalScalarType;
348  using SC = typename Container<MatrixType>::scalar_type;
349  using LO = typename Container<MatrixType>::local_ordinal_type;
350  using GO = typename Container<MatrixType>::global_ordinal_type;
351  using NO = typename Container<MatrixType>::node_type;
353  using typename Container<MatrixType>::import_type;
354  using typename Container<MatrixType>::row_matrix_type;
355  using typename Container<MatrixType>::crs_matrix_type;
356  using typename Container<MatrixType>::block_crs_matrix_type;
357  using typename Container<MatrixType>::mv_type;
358  using typename Container<MatrixType>::vector_type;
359  using typename Container<MatrixType>::map_type;
360  using typename Container<MatrixType>::ISC;
362  using LSC = LocalScalarType;
363  using LISC = typename Kokkos::Details::ArithTraits<LSC>::val_type;
364 
365  using local_mv_type = Tpetra::MultiVector<LSC, LO, GO, NO>;
366 
367  using typename Container<MatrixType>::HostView;
368  using HostViewLocal = typename local_mv_type::dual_view_type::t_host;
369  using HostSubviewLocal = Kokkos::View<LISC**, Kokkos::LayoutStride, typename HostViewLocal::memory_space>;
370 
371  static_assert(std::is_same<MatrixType, row_matrix_type>::value,
372  "Ifpack2::Container: Please use MatrixType = Tpetra::RowMatrix.");
374  //
375 
382  LO translateRowToCol(LO row);
383 
386 
391  virtual void extract() = 0;
392 
393 public:
394 
396  const Teuchos::Array<Teuchos::Array<LO> >& partitions,
397  bool pointIndexed);
398 
400  virtual ~ContainerImpl();
401 
410  virtual void initialize () = 0;
411 
420  virtual void compute () = 0;
421 
423  virtual void setParameters (const Teuchos::ParameterList& List);
424 
437  virtual void
438  apply(HostView X,
439  HostView Y,
440  int blockIndex,
442  SC alpha = Teuchos::ScalarTraits<SC>::one(),
443  SC beta = Teuchos::ScalarTraits<SC>::zero()) const;
444 
446  virtual void
448  HostView Y,
449  HostView D,
450  int blockIndex,
452  SC alpha = Teuchos::ScalarTraits<SC>::one(),
453  SC beta = Teuchos::ScalarTraits<SC>::zero()) const;
454 
457  // The underlying container implements the splitting <tt>A = D + R</tt>. Only
458  // it can have efficient access to D and R, as these are constructed in the
459  // symbolic and numeric phases.
460  //
461  // This is the first performance-portable implementation of a block
462  // relaxation, and it is supported currently only by BlockTriDiContainer.
463  virtual void applyInverseJacobi (const mv_type& /* X */, mv_type& /* Y */,
464  SC dampingFactor,
465  bool /* zeroStartingSolution = false */,
466  int /* numSweeps = 1 */) const;
467 
469  void applyMV (mv_type& X, mv_type& Y) const;
470 
472  void weightedApplyMV (mv_type& X,
473  mv_type& Y,
474  vector_type& W) const;
475 
476  virtual void clearBlocks();
477 
479  virtual std::ostream& print (std::ostream& os) const = 0;
480 
483  static std::string getName();
484 
485 protected:
486  //Do Gauss-Seidel on only block i (this is used by DoGaussSeidel and DoSGS)
487  void DoGSBlock(HostView X, HostView Y, HostView Y2, HostView Resid,
488  SC dampingFactor, LO i) const;
489 
495  virtual void
496  solveBlock(HostSubviewLocal X,
497  HostSubviewLocal Y,
498  int blockIndex,
499  Teuchos::ETransp mode,
500  const LSC alpha,
501  const LSC beta) const;
502 
504  mutable HostViewLocal X_local_; //length: blockRows_.size()
505  mutable HostViewLocal Y_local_; //length: blockRows_.size()
506 
517  mutable HostViewLocal weightedApplyScratch_;
518 
520  mutable std::vector<HostSubviewLocal> X_localBlocks_;
521 
523  mutable std::vector<HostSubviewLocal> Y_localBlocks_;
524 };
525 
526 namespace Details {
532  template<typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, typename Node>
533  struct StridedRowView
534  {
535  using SC = Scalar;
536  using LO = LocalOrdinal;
538  StridedRowView(const SC* vals_, const LO* inds_, int blockSize_, size_t nnz_);
539 
542 
543  SC val(size_t i) const;
544  LO ind(size_t i) const;
545 
546  size_t size() const;
547 
548  private:
549  const SC* vals;
550  const LO* inds;
551  int blockSize;
552  size_t nnz;
553  //These arrays are only used if the inputMatrix_ doesn't support row views.
554  Teuchos::Array<SC> valsCopy;
555  Teuchos::Array<LO> indsCopy;
556  };
557 } // namespace Details
558 
559 } // namespace Ifpack2
560 
562 template <class MatrixType>
563 std::ostream& operator<<(std::ostream& os, const Ifpack2::Container<MatrixType>& obj);
564 
565 namespace Teuchos {
566 
571 template<class MatrixType>
572 class TEUCHOSCORE_LIB_DLL_EXPORT TypeNameTraits< ::Ifpack2::Container<MatrixType> >
573 {
574  public:
575  static std::string name () {
576  return std::string ("Ifpack2::Container<") +
578  }
579 
580  static std::string concreteName (const ::Ifpack2::Container<MatrixType>&) {
581  return name ();
582  }
583 };
584 
585 } // namespace Teuchos
586 
587 #endif // IFPACK2_CONTAINER_HPP
virtual void applyInverseJacobi(const mv_type &, mv_type &, SC dampingFactor, bool, int) const =0
Compute Y := (1 - a) Y + a D^{-1} (X - R*Y).
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 =0
Compute Y := alpha * diag(D) * M^{-1} (diag(D) * X) + beta*Y.
virtual void initialize()=0
Do all set-up operations that only require matrix structure.
bool isComputed() const
Whether the container has been successfully computed.
Definition: Ifpack2_Container_def.hpp:170
std::vector< HostSubviewLocal > X_localBlocks_
Views for holding pieces of X corresponding to each block.
Definition: Ifpack2_Container_decl.hpp:520
Teuchos::Array< LO > blockRows_
Local indices of the rows of the input matrix that belong to this block.
Definition: Ifpack2_Container_decl.hpp:293
void setBlockSizes(const Teuchos::Array< Teuchos::Array< LO > > &partitions)
Initialize arrays with information about block sizes.
Definition: Ifpack2_Container_def.hpp:125
bool IsInitialized_
If true, the container has been successfully initialized.
Definition: Ifpack2_Container_decl.hpp:319
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
std::vector< HostSubviewLocal > Y_localBlocks_
Views for holding pieces of Y corresponding to each block.
Definition: Ifpack2_Container_decl.hpp:523
virtual void initialize()=0
Do all set-up operations that only require matrix structure.
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
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
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
HostViewLocal weightedApplyScratch_
Definition: Ifpack2_Container_decl.hpp:517
virtual void compute()=0
Extract the local diagonal blocks and prepare the solver.
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
InverseType::scalar_type LSC
The internal representation of LocalScalarType in Kokkos::View.
Definition: Ifpack2_Container_decl.hpp:362
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
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
HostViewLocal X_local_
Scratch vectors used in apply().
Definition: Ifpack2_Container_decl.hpp:504
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
Teuchos::RCP< vector_type > Diag_
Diagonal elements.
Definition: Ifpack2_Container_decl.hpp:299
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
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 =0
Compute Y := alpha * M^{-1} X + beta*Y.
bool IsComputed_
If true, the container has been successfully computed.
Definition: Ifpack2_Container_decl.hpp:321
virtual void compute()=0
Extract the local diagonal blocks and prepare the solver.
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
virtual void extract()=0
Extract the submatrices identified by the local indices set by the constructor. BlockMatrix may be an...
virtual void setParameters(const Teuchos::ParameterList &List)=0
Set parameters, if any.
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
virtual std::ostream & print(std::ostream &os) const =0
Print basic information about the container to os.
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
Teuchos::Array< LO > blockOffsets_
Starting index in blockRows_ of local row indices for each block.
Definition: Ifpack2_Container_decl.hpp:297
virtual std::ostream & print(std::ostream &os) const =0
Print basic information about the container to os.
static std::string name()
Teuchos::RCP< const crs_matrix_type > inputCrsMatrix_
The input matrix, dynamic cast to CrsMatrix. May be null.
Definition: Ifpack2_Container_decl.hpp:285