Xpetra_CrsMatrixWrap_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

// WARNING: This code is experimental. Backwards compatibility should not be expected.

#ifndef XPETRA_CRSMATRIXWRAP_DECL_HPP
#define XPETRA_CRSMATRIXWRAP_DECL_HPP

#include <Tpetra_KokkosCompat_DefaultNode.hpp>

#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"

#include "Xpetra_MultiVector_decl.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"

#include "Xpetra_Matrix.hpp"
#include "Xpetra_TpetraOperator.hpp"
#include "Xpetra_TpetraBlockCrsMatrix.hpp"
#include "Xpetra_TpetraRowMatrix.hpp"
#include "Xpetra_TpetraCrsMatrix.hpp"

#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_Hashtable.hpp>

#include <Tpetra_Filter.hpp>

namespace Xpetra {

typedef std::string viewLabel_t;

template <class Scalar,
          class LocalOrdinal,
          class GlobalOrdinal,
          class Node = Tpetra::KokkosClassic::DefaultNode::DefaultNodeType>
class CrsMatrixWrap : public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
  typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> Map;
  typedef Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> CrsMatrix;
  typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> Matrix;
  typedef Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> CrsGraph;
#ifdef HAVE_XPETRA_TPETRA
  typedef Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpetraCrsMatrix;
#endif
  typedef Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> CrsMatrixFactory;
  typedef Xpetra::MatrixView<Scalar, LocalOrdinal, GlobalOrdinal, Node> MatrixView;
#ifdef HAVE_XPETRA_TPETRA
  using local_matrix_type        = typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type;
  using local_matrix_device_type = typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type;
  using local_matrix_host_type   = typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_host_type;
#endif

 public:


  CrsMatrixWrap(const RCP<const Map> &rowMap);

  CrsMatrixWrap(const RCP<const Map> &rowMap,
                size_t maxNumEntriesPerRow);

  CrsMatrixWrap(const RCP<const Map> &rowMap,
                const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc);

  CrsMatrixWrap(const RCP<const Map> &rowMap, const RCP<const Map> &colMap, size_t maxNumEntriesPerRow);

  CrsMatrixWrap(const RCP<const Map> &rowMap, const RCP<const Map> &colMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc);

#ifdef HAVE_XPETRA_TPETRA
  CrsMatrixWrap(const RCP<const Map> &rowMap, const RCP<const Map> &colMap, const local_matrix_type &lclMatrix, const Teuchos::RCP<Teuchos::ParameterList> &params = null);

  CrsMatrixWrap(const local_matrix_type &lclMatrix, const RCP<const Map> &rowMap, const RCP<const Map> &colMap,
                const RCP<const Map> &domainMap = Teuchos::null, const RCP<const Map> &rangeMap = Teuchos::null,
                const Teuchos::RCP<Teuchos::ParameterList> &params = null);
#else
#ifdef __GNUC__
#warning "Xpetra Kokkos interface for CrsMatrix is enabled (HAVE_XPETRA_KOKKOS_REFACTOR) but Tpetra is disabled. The Kokkos interface needs Tpetra to be enabled, too."
#endif
#endif

  CrsMatrixWrap(RCP<CrsMatrix> matrix);

  CrsMatrixWrap(const RCP<const CrsGraph> &graph, const RCP<ParameterList> &paramList = Teuchos::null);

  CrsMatrixWrap(const RCP<const CrsGraph> &graph,
                typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type::values_type &values,
                const RCP<ParameterList> &paramList = Teuchos::null);

  virtual ~CrsMatrixWrap();





  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 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);




  void resumeFill(const RCP<ParameterList> &params = null);

  void fillComplete(const RCP<const Map> &domainMap, const RCP<const Map> &rangeMap, const RCP<Teuchos::ParameterList> &params = null);

  // TODO : Get ride of "Tpetra"::OptimizeOption
  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;


  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(Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) const;

  void getLocalDiagOffsets(Teuchos::ArrayRCP<size_t> &offsets) const;

  void getLocalDiagCopy(Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag, const Teuchos::ArrayView<const size_t> &offsets) const;

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

  void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x);

  void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x);

  bool haveGlobalConstants() const;





  // TODO virtual=0 // TODO: Add default parameters ?
  //   void multiply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, Scalar alpha, Scalar beta) const {
  //      matrixData_->multiply(X, Y, trans, alpha, beta);
  //   }





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

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

  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> getDomainMap() const;

  const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> getRangeMap() const;

  const RCP<const Map> &getColMap() const;

  const RCP<const Map> &getColMap(viewLabel_t viewLabel) const;

  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map> &newMap);


  //{@

  const Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> getMap() const;

  void doImport(const Matrix &source,
                const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM);

  void doExport(const Matrix &dest,
                const Xpetra::Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM);

  void doImport(const Matrix &source,
                const Xpetra::Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM);

  void doExport(const Matrix &dest,
                const Xpetra::Export<LocalOrdinal, GlobalOrdinal, Node> &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);

#ifdef HAVE_XPETRA_TPETRA
  virtual local_matrix_type getLocalMatrixDevice() const;
#if KOKKOS_VERSION >= 40799
  virtual typename local_matrix_type::host_mirror_type getLocalMatrixHost() const;
#else
  virtual typename local_matrix_type::HostMirror getLocalMatrixHost() const;
#endif
#else
#ifdef __GNUC__
#warning "Xpetra Kokkos interface for CrsMatrix is enabled (HAVE_XPETRA_KOKKOS_REFACTOR) but Tpetra is disabled. The Kokkos interface needs Tpetra to be enabled, too."
#endif
#endif

  // JG: Added:

  bool hasCrsGraph() const;

  RCP<const CrsGraph> getCrsGraph() const;

  RCP<CrsMatrix> getCrsMatrix() const;

  LocalOrdinal GetStorageBlockSize() const;

  void residual(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X,
                const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B,
                MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &R) const;

  void replaceCrsMatrix(RCP<CrsMatrix> &M);

 private:
  // Default view is created after fillComplete()
  // Because ColMap might not be available before fillComplete().
  void CreateDefaultView();

  // The colMap can be <tt>null</tt> until fillComplete() is called. The default view of the Matrix have to be updated when fillComplete() is called.
  // If CrsMatrix::fillComplete() have been used instead of CrsMatrixWrap::fillComplete(), the default view is updated when getColMap() is called.
  void updateDefaultView() const;
  // The boolean finalDefaultView_ keep track of the status of the default view (= already updated or not)
  // See also CrsMatrixWrap::updateDefaultView()
  mutable bool finalDefaultView_;

  // The underlying matrix object
  RCP<CrsMatrix> matrixData_;

};  // class CrsMatrixWrap

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toCrsMatrix(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &matrix) {
  return Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(matrix, true)->getCrsMatrix();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toCrsMatrix(const Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &matrix) {
  return Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(matrix, true)->getCrsMatrix();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toCrsMatrix(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &matrix) {
  return *dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(matrix).getCrsMatrix();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toCrsMatrix(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &matrix) {
  return *dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(matrix).getCrsMatrix();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toXpetra(const Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto xpTpCrs = Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
  auto xpCrs   = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpTpCrs, true);
  return rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrs));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toXpetra(Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto xpTpCrs = Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
  auto xpCrs   = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpTpCrs, true);
  return rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrs));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toXpetra(Teuchos::RCP<Tpetra::FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto xpTpCrs = Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
  auto xpCrs   = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpTpCrs, true);
  return rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrs));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toXpetra(Teuchos::RCP<const Tpetra::FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto xpTpCrs = Teuchos::rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
  auto xpCrs   = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(xpTpCrs, true);
  return rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(xpCrs));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetra(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return toTpetra(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetra(const Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return toTpetra(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetra(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return toTpetra(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetra(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return toTpetra(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraBlock(const Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A, true)->getTpetra_BlockCrsMatrix();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraBlock(const Teuchos::RCP<const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return Teuchos::rcp_dynamic_cast<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A, true)->getTpetra_BlockCrsMatrixNonConst();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetraBlock(Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return *(dynamic_cast<Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(A).getTpetra_BlockCrsMatrixNonConst());
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetraBlock(const Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return *(dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(A).getTpetra_BlockCrsMatrix());
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraBlock(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return toTpetraBlock(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraBlock(const Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  return toTpetraBlock(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetraBlock(Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return toTpetraBlock(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
const Tpetra::BlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &
toTpetraBlock(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
  return toTpetraBlock(toCrsMatrix(A));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraRowMatrix(const Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &mat) {
  auto crsMat  = toCrsMatrix(mat);
  auto tmp_Crs = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
  if (!tmp_Crs.is_null()) {
    return tmp_Crs->getTpetra_CrsMatrixNonConst();
  } else {
    auto tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
    if (tmp_BlockCrs.is_null())
      throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
    return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
  }
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraRowMatrix(const Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &mat) {
  auto crsMat  = toCrsMatrix(mat);
  auto tmp_Crs = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
  if (!tmp_Crs.is_null()) {
    return tmp_Crs->getTpetra_CrsMatrix();
  } else {
    auto tmp_BlockCrs = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsMat);
    if (tmp_BlockCrs.is_null())
      throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
    return tmp_BlockCrs->getTpetra_BlockCrsMatrix();
  }
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraRowMatrix(const Teuchos::RCP<Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto mat  = Teuchos::rcp_dynamic_cast<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A);
  auto rmat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A);
  if (!mat.is_null()) {
    return toTpetraRowMatrix(mat);
  } else if (!rmat.is_null()) {
    return rmat->getTpetra_RowMatrixNonConst();
  } else {
    auto tpOp = Teuchos::rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A, true);
    auto tOp  = tpOp->getOperator();
    auto tRow = Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tOp, true);
    return tRow;
  }
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
toTpetraRowMatrix(const Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A) {
  auto mat  = Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A);
  auto rmat = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A);
  if (!mat.is_null()) {
    return toTpetraRowMatrix(mat);
  } else if (!rmat.is_null()) {
    return rmat->getTpetra_RowMatrix();
  } else {
    auto tpOp = Teuchos::rcp_dynamic_cast<const Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A, true);
    auto tOp  = tpOp->getOperator();
    auto tRow = Teuchos::rcp_dynamic_cast<const Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tOp, true);
    return tRow;
  }
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_vals(const RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &crsA, const filter_type &filter) {
  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsFiltered;
#ifdef HAVE_XPETRA_TPETRA
  auto tpCrsA = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsA);
  if (!tpCrsA.is_null()) {
    auto tpCrsFiltered = Tpetra::applyFilter_vals(*tpCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsFiltered);
  }
#endif
#ifdef HAVE_XPETRA_EPETRA
  auto epCrsA = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(crsA);
  if (!epCrsA.is_null()) {
    auto epCrsFiltered = Tpetra::applyFilter_vals(*epCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(epCrsFiltered);
  }
#endif
  return crsFiltered;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_vals(const RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A, const filter_type &filter) {
  auto crsA        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A)->getCrsMatrix();
  auto crsFiltered = applyFilter_vals(crsA, filter);
  auto filtered    = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(crsFiltered));
  return filtered;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_LID(const RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &crsA, const filter_type &filter) {
  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsFiltered;
#ifdef HAVE_XPETRA_TPETRA
  auto tpCrsA = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsA);
  if (!tpCrsA.is_null()) {
    auto tpCrsFiltered = Tpetra::applyFilter_LID(*tpCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsFiltered);
  }
#endif
#ifdef HAVE_XPETRA_EPETRA
  auto epCrsA = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(crsA);
  if (!epCrsA.is_null()) {
    auto epCrsFiltered = Tpetra::applyFilter_LID(*epCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(epCrsFiltered);
  }
#endif
  return crsFiltered;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_LID(const RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A, const filter_type &filter) {
  auto crsA        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A)->getCrsMatrix();
  auto crsFiltered = applyFilter_LID(crsA, filter);
  auto filtered    = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(crsFiltered));
  return filtered;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_GID(const RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &crsA, const filter_type &filter) {
  Teuchos::RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> crsFiltered;
#ifdef HAVE_XPETRA_TPETRA
  auto tpCrsA = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(crsA);
  if (!tpCrsA.is_null()) {
    auto tpCrsFiltered = Tpetra::applyFilter_GID(*tpCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(tpCrsFiltered);
  }
#endif
#ifdef HAVE_XPETRA_EPETRA
  auto epCrsA = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node>>(crsA);
  if (!epCrsA.is_null()) {
    auto epCrsFiltered = Tpetra::applyFilter_GID(*epCrsA, filter);
    crsFiltered        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(epCrsFiltered);
  }
#endif
  return crsFiltered;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node, class filter_type>
RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> applyFilter_GID(const RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>> &A, const filter_type &filter) {
  auto crsA        = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>>(A)->getCrsMatrix();
  auto crsFiltered = applyFilter_GID(crsA, filter);
  auto filtered    = rcp(new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(crsFiltered));
  return filtered;
}

}  // namespace Xpetra

#define XPETRA_CRSMATRIXWRAP_SHORT
#endif  // XPETRA_CRSMATRIXWRAP_DECL_HPP

// NOTE: if CrsMatrix and VbrMatrix share a common interface for fillComplete() etc, I can move some stuff in Xpetra_Matrix.hpp
// TODO: getUnderlyingMatrix() method