Xpetra_BlockedCrsMatrix_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//             Xpetra: A linear algebra interface package
//
// Copyright 2012 NTESS and the Xpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef XPETRA_BLOCKEDCRSMATRIX_DECL_HPP
#define XPETRA_BLOCKEDCRSMATRIX_DECL_HPP
 
#include <Tpetra_KokkosCompat_DefaultNode.hpp>
 
#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_Hashtable.hpp>
 
#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"
 
#include "Xpetra_MapFactory.hpp"
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_BlockedMultiVector.hpp"
#include "Xpetra_MultiVectorFactory.hpp"
#include "Xpetra_BlockedVector.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"
 
#include "Xpetra_MapExtractor.hpp"
#include "Xpetra_MapExtractorFactory.hpp"
 
#include "Xpetra_Matrix.hpp"
#include "Xpetra_MatrixFactory.hpp"
#include "Xpetra_CrsMatrixWrap.hpp"
 
#ifdef HAVE_XPETRA_THYRA
#include <Thyra_ProductVectorSpaceBase.hpp>
#include <Thyra_VectorSpaceBase.hpp>
#include <Thyra_LinearOpBase.hpp>
#include <Thyra_BlockedLinearOpBase.hpp>
#include <Thyra_PhysicallyBlockedLinearOpBase.hpp>
#include "Xpetra_ThyraUtils.hpp"
#endif
 
#include "Xpetra_VectorFactory.hpp"

namespace Xpetra {
 
#ifdef HAVE_XPETRA_THYRA
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
class ThyraUtils;
#endif
 
typedef std::string viewLabel_t;
 
template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
class BlockedCrsMatrix : public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
 public:
  typedef Scalar scalar_type;
  typedef LocalOrdinal local_ordinal_type;
  typedef GlobalOrdinal global_ordinal_type;
  typedef Node node_type;
 
 private:
#undef XPETRA_BLOCKEDCRSMATRIX_SHORT
#include "Xpetra_UseShortNames.hpp"
 
 public:



  BlockedCrsMatrix(const Teuchos::RCP<const BlockedMap>& rangeMaps,
                   const Teuchos::RCP<const BlockedMap>& domainMaps,
                   size_t numEntriesPerRow);


  BlockedCrsMatrix(Teuchos::RCP<const MapExtractor>& rangeMapExtractor,
                   Teuchos::RCP<const MapExtractor>& domainMapExtractor,
                   size_t numEntriesPerRow);
 
#ifdef HAVE_XPETRA_THYRA

  BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& /* comm */);
 
 private:

  Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > >& subMaps);
 
 public:
#endif

  virtual ~BlockedCrsMatrix();





    \post <tt>isGloballyIndexed() == true</tt>
 
    \note If \c globalRow does not belong to the matrix on this node, then it
    will be communicated to the appropriate node when globalAssemble() is
    called (which will, at the latest, occur during the next call to
    fillComplete().) Otherwise, the entries will be inserted in the local
    matrix.
 
    \note If the matrix row already contains values at the indices
    corresponding to values in \c cols, then the new values will be summed
    with the old values; this may happen at insertion or during the next call
    to fillComplete().
 
    \note If <tt>hasColMap() == true</tt>, only (cols[i],vals[i]) where
    cols[i] belongs to the column map on this node will be inserted into the
    matrix.
    */
  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals);


  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals);
 
  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map>& newMap);



  void replaceGlobalValues(GlobalOrdinal globalRow,
                           const ArrayView<const GlobalOrdinal>& cols,
                           const ArrayView<const Scalar>& vals);


  void replaceLocalValues(LocalOrdinal localRow,
                          const ArrayView<const LocalOrdinal>& cols,
                          const ArrayView<const Scalar>& vals);

  virtual void setAllToScalar(const Scalar& alpha);

  void scale(const Scalar& alpha);




    After calling fillComplete(), resumeFill() must be called before initiating any changes to the matrix.
 
    For BlockedCrsMatrix objects we call the routine iteratively for all sub-blocks.
 
    resumeFill() may be called repeatedly.
    */
  void resumeFill(const RCP<ParameterList>& params = null);

 
    Note: for blocked operators the specified domain and range maps have no meaning.
          We just call fillComplete for all underlying blocks
    */
  void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null);


    \note This method calls fillComplete( getRowMap(), getRowMap(), os ).
 
    \pre  <tt>isFillActive() == true<tt>
    \pre <tt>isFillComplete()() == false<tt>
 
    \post <tt>isFillActive() == false<tt>
    \post <tt>isFillComplete() == true<tt>
    \post if <tt>os == DoOptimizeStorage<tt>, then <tt>isStorageOptimized() == true</tt>
    */
  void fillComplete(const RCP<ParameterList>& params = null);




  global_size_t getGlobalNumRows() const;


  global_size_t getGlobalNumCols() const;

  size_t getLocalNumRows() const;

  global_size_t getGlobalNumEntries() const;

  size_t getLocalNumEntries() const;


  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;


  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;


  size_t getGlobalMaxNumRowEntries() const;


  size_t getLocalMaxNumRowEntries() const;


  bool isLocallyIndexed() const;


  bool isGloballyIndexed() const;

  bool isFillComplete() const;


  virtual void getLocalRowCopy(LocalOrdinal LocalRow,
                               const ArrayView<LocalOrdinal>& Indices,
                               const ArrayView<Scalar>& Values,
                               size_t& NumEntries) const;


  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const;


  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const;


  void getLocalDiagCopy(Vector& diag) const;

  void leftScale(const Vector& x);

  void rightScale(const Vector& x);

  virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;

  virtual bool haveGlobalConstants() const;









  virtual void apply(const MultiVector& X, MultiVector& Y, Teuchos::ETransp mode, Scalar alpha, Scalar beta, bool sumInterfaceValues,
                     const RCP<Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> >& regionInterfaceImporter,
                     const Teuchos::ArrayRCP<LocalOrdinal>& regionInterfaceLIDs) const;


  virtual void apply(const MultiVector& X, MultiVector& Y,
                     Teuchos::ETransp mode = Teuchos::NO_TRANS,
                     Scalar alpha          = ScalarTraits<Scalar>::one(),
                     Scalar beta           = ScalarTraits<Scalar>::zero()) const;

  RCP<const Map> getFullDomainMap() const;

  RCP<const BlockedMap> getBlockedDomainMap() const;

  const RCP<const Map> getDomainMap() const;

  RCP<const Map> getDomainMap(size_t i) const;

  RCP<const Map> getDomainMap(size_t i, bool bThyraMode) const;

  RCP<const Map> getFullRangeMap() const;

  RCP<const BlockedMap> getBlockedRangeMap() const;

  const RCP<const Map> getRangeMap() const;

  RCP<const Map> getRangeMap(size_t i) const;

  RCP<const Map> getRangeMap(size_t i, bool bThyraMode) const;

  RCP<const MapExtractor> getRangeMapExtractor() const;

  RCP<const MapExtractor> getDomainMapExtractor() const;


  //{@

  virtual void bgs_apply(
      const MultiVector& X,                                 
      MultiVector& Y,                                       
      size_t row,                                           
      Teuchos::ETransp mode = Teuchos::NO_TRANS,            
      Scalar alpha          = ScalarTraits<Scalar>::one(),  
      Scalar beta           = ScalarTraits<Scalar>::zero()  
  ) const;


  //{@

  const Teuchos::RCP<const Map> getMap() const;

  void doImport(const Matrix& source, const Import& importer, CombineMode CM);

  void doExport(const Matrix& dest, const Import& importer, CombineMode CM);

  void doImport(const Matrix& source, const Export& exporter, CombineMode CM);

  void doExport(const Matrix& dest, const Export& exporter, CombineMode CM);
 
  // @}



  std::string description() const;

  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const;


  void setObjectLabel(const std::string& objectLabel);

  bool hasCrsGraph() const;

  RCP<const CrsGraph> getCrsGraph() const;



 
  virtual bool isDiagonal() const;

  virtual size_t Rows() const;

  virtual size_t Cols() const;

  Teuchos::RCP<Matrix> getCrsMatrix() const;

  Teuchos::RCP<Matrix> getInnermostCrsMatrix();

  Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const;

  // void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
  void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat);

  // NOTE: This is a rather expensive operation, since all blocks are copied
  // into a new big CrsMatrix
  Teuchos::RCP<Matrix> Merge() const;
 
  typedef typename CrsMatrix::local_matrix_type local_matrix_type;
  local_matrix_type getLocalMatrixDevice() const;
  typename local_matrix_type::HostMirror getLocalMatrixHost() const;
 
#ifdef HAVE_XPETRA_THYRA
  Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > getThyraOperator();
#endif
  LocalOrdinal GetStorageBlockSize() const;

  void residual(const MultiVector& X,
                const MultiVector& B,
                MultiVector& R) const;
 
 private:


  void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const;

 
  // Default view is created after fillComplete()
  // Because ColMap might not be available before fillComplete().
  void CreateDefaultView();
 
 private:
  bool is_diagonal_;                             
  Teuchos::RCP<const MapExtractor> domainmaps_;  
  Teuchos::RCP<const MapExtractor> rangemaps_;   
 
  std::vector<Teuchos::RCP<Matrix> > blocks_;  
#ifdef HAVE_XPETRA_THYRA
  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > thyraOp_;  
#endif
  bool bRangeThyraMode_;   
  bool bDomainThyraMode_;  
};
 
}  // namespace Xpetra
 
#define XPETRA_BLOCKEDCRSMATRIX_SHORT
#endif /* XPETRA_BLOCKEDCRSMATRIX_DECL_HPP */