Xpetra_EpetraCrsMatrix.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_EPETRACRSMATRIX_HPP
#define XPETRA_EPETRACRSMATRIX_HPP
 
/* this file is automatically generated - do not edit (see script/epetra.py) */
 
#ifdef HAVE_XPETRA_TPETRA
#include <Kokkos_Core.hpp>
#endif
 
#include "Xpetra_EpetraConfigDefs.hpp"
 
#include "Xpetra_CrsMatrix.hpp"
 
#include <Epetra_CrsMatrix.h>
#include <Epetra_Map.h>
 
#include "Xpetra_EpetraMap.hpp"
#include "Xpetra_EpetraVector.hpp"
#include "Xpetra_EpetraMultiVector.hpp"
#include "Xpetra_EpetraCrsGraph.hpp"
 
#include "Xpetra_MapFactory.hpp"
 
#include "Xpetra_Utils.hpp"
#include "Xpetra_Exceptions.hpp"
 
#if defined(XPETRA_ENABLE_DEPRECATED_CODE)
#ifdef __GNUC__
#if defined(Xpetra_SHOW_DEPRECATED_WARNINGS)
#warning "The header file Trilinos/packages/xpetra/src/CrsMatrix/Xpetra_EpetraCrsMatrix.hpp is deprecated."
#endif
#endif
#else
#error "The header file Trilinos/packages/xpetra/src/CrsMatrix/Xpetra_EpetraCrsMatrix.hpp is deprecated."
#endif
 
namespace Xpetra {
 
// general implementation: empty stub
template <class EpetraGlobalOrdinal, class Node>
class XPETRA_DEPRECATED EpetraCrsMatrixT
  : public CrsMatrix<double, int, EpetraGlobalOrdinal, Node> {
  typedef EpetraGlobalOrdinal GlobalOrdinal;
  typedef typename CrsMatrix<double, int, GlobalOrdinal, Node>::scalar_type Scalar;
  typedef typename CrsMatrix<double, int, GlobalOrdinal, Node>::local_ordinal_type LocalOrdinal;
 
#ifdef HAVE_XPETRA_TPETRA
 public:
  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;
  typedef typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::node_type node_type;
 
 private:
#endif
 
#if KOKKOS_VERSION >= 40799
  using ATS = KokkosKernels::ArithTraits<Scalar>;
#else
  using ATS = Kokkos::ArithTraits<Scalar>;
#endif
  using impl_scalar_type = typename ATS::val_type;
 
 public:
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> &params = Teuchos::null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> &params = Teuchos::null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> &params = Teuchos::null) { TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with Node=Tpetra::KokkosCompat::KokkosSerialWrapperNode."); }
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> &params = Teuchos::null) { TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with Node=Tpetra::KokkosCompat::KokkosSerialWrapperNode."); }
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > &graph, const Teuchos::RCP<Teuchos::ParameterList> &params = Teuchos::null) { TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError, "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with Node=Tpetra::KokkosCompat::KokkosSerialWrapperNode."); }
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &importer,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &RowImporter,
                   const Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > DomainImporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &RowExporter,
                   const Teuchos::RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > DomainExporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
#ifdef HAVE_XPETRA_TPETRA
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
                   const local_matrix_type &lclMatrix,
                   const Teuchos::RCP<Teuchos::ParameterList> &params = null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  EpetraCrsMatrixT(
      const local_matrix_type &lclMatrix,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
      const Teuchos::RCP<Teuchos::ParameterList> &params                           = null) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
#endif
  virtual ~EpetraCrsMatrixT() {}
 
  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) {}
  void setAllToScalar(const Scalar &alpha) {}
  void scale(const Scalar &alpha) {}
  void allocateAllValues(size_t numNonZeros, ArrayRCP<size_t> &rowptr, ArrayRCP<LocalOrdinal> &colind, ArrayRCP<Scalar> &values) {}
  void setAllValues(const ArrayRCP<size_t> &rowptr, const ArrayRCP<LocalOrdinal> &colind, const ArrayRCP<Scalar> &values) {}
  void getAllValues(ArrayRCP<const size_t> &rowptr, ArrayRCP<const LocalOrdinal> &colind, ArrayRCP<const Scalar> &values) const {}
  void getAllValues(ArrayRCP<Scalar> &values) {}
  bool haveGlobalConstants() const { return true; }
  void expertStaticFillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                                const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                                const RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &importer = Teuchos::null,
                                const RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > &exporter = Teuchos::null,
                                const RCP<ParameterList> &params                                      = Teuchos::null) {}
 
  void resumeFill(const RCP<ParameterList> &params = null) {}
  void fillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap, const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap, const RCP<ParameterList> &params = null) {}
  void fillComplete(const RCP<ParameterList> &params = null) {}
  void replaceDomainMapAndImporter(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newDomainMap, Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &newImporter) {}
  const RCP<const Comm<int> > getComm() const { return Teuchos::null; }
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRowMap() const { return Teuchos::null; }
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getColMap() const { return Teuchos::null; }
  RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > getCrsGraph() const { return Teuchos::null; }
  global_size_t getGlobalNumRows() const { return 0; }
  global_size_t getGlobalNumCols() const { return 0; }
  size_t getLocalNumRows() const { return 0; }
  size_t getLocalNumCols() const { return 0; }
  global_size_t getGlobalNumEntries() const { return 0; }
  size_t getLocalNumEntries() const { return 0; }
  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const { return 0; }
  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const { return 0; }
  size_t getGlobalMaxNumRowEntries() const { return 0; }
  size_t getLocalMaxNumRowEntries() const { return 0; }
  bool isLocallyIndexed() const { return false; }
  bool isGloballyIndexed() const { return false; }
  bool isFillComplete() const { return false; }
  bool isFillActive() const { return false; }
  typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const { return ScalarTraits<Scalar>::magnitude(ScalarTraits<Scalar>::zero()); }
  bool supportsRowViews() const { return false; }
  void getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView<LocalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {}
  void getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView<GlobalOrdinal> &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<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) const {}
  void getLocalDiagOffsets(Teuchos::ArrayRCP<size_t> &offsets) const {}
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag, const Teuchos::ArrayView<const size_t> &offsets) const {}
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag, const Kokkos::View<const size_t *, typename Node::device_type, Kokkos::MemoryUnmanaged> &offsets) const {}
  void replaceDiag(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) {}
  void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x){};
  void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x){};
 
  void apply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X, MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &Y, Teuchos::ETransp mode = Teuchos::NO_TRANS, Scalar alpha = ScalarTraits<Scalar>::one(), Scalar beta = ScalarTraits<Scalar>::zero()) const {}
  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 Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const { return Teuchos::null; }
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const { return Teuchos::null; }
 
  std::string description() const { return std::string(""); }
  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {}
  void setObjectLabel(const std::string &objectLabel) {}
 
  EpetraCrsMatrixT(const EpetraCrsMatrixT &matrix) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
 
  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const { return Teuchos::null; }
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {}
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {}
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {}
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {}
  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newMap) {}
 
  bool hasMatrix() const { return false; }
  EpetraCrsMatrixT(const Teuchos::RCP<Epetra_CrsMatrix> &mtx) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
  RCP<const Epetra_CrsMatrix> getEpetra_CrsMatrix() const { return Teuchos::null; }
  RCP<Epetra_CrsMatrix> getEpetra_CrsMatrixNonConst() const { return Teuchos::null; }  // TODO: remove
#ifdef HAVE_XPETRA_TPETRA
  local_matrix_type getLocalMatrixDevice() const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }
#if KOKKOS_VERSION >= 40799
  typename local_matrix_type::host_mirror_type getLocalMatrixHost() const {
#else
  typename local_matrix_type::HostMirror getLocalMatrixHost() const {
#endif
      TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                                 "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
}
 
void
setAllValues(const typename local_matrix_type::row_map_type &ptr,
             const typename local_matrix_type::StaticCrsGraphType::entries_type::non_const_type &ind,
             const typename local_matrix_type::values_type &val) {
  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                             "Xpetra::EpetraCrsMatrix only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
}
 
LocalOrdinal GetStorageBlockSize() const { return 1; }
 
#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
 
void residual(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X,
              const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B,
              MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &R) const {
  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one;
  apply(X, R);
  R.update(one, B, negone);
}
 
};  // namespace Xpetra
 
#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
template <>
class EpetraCrsMatrixT<int, EpetraNode>
  : public CrsMatrix<double, int, int, EpetraNode> {
  typedef int GlobalOrdinal;
  typedef EpetraNode Node;
  typedef CrsMatrix<double, int, GlobalOrdinal, Node>::scalar_type Scalar;
  typedef CrsMatrix<double, int, GlobalOrdinal, Node>::local_ordinal_type LocalOrdinal;
 
  // The following typedefs are used by the Kokkos interface
#ifdef HAVE_XPETRA_TPETRA
 public:
  typedef typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type 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;
  typedef typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::node_type node_type;
 
 private:
#endif
 
 public:
 
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), 0, false)))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), maxNumEntriesPerRow, true)))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : isFillResumed_(false) {
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), numEntriesPerRowToAlloc.getRawPtr(), true));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), maxNumEntriesPerRow, true)))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : isFillResumed_(false) {
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), numEntriesPerRowToAlloc.getRawPtr(), true));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > &graph, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(graph))))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &importer,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tImporter.getEpetra_Import(), myDomainMap, myRangeMap, restrictComm));
    if (restrictComm && mtx_->NumMyRows() == 0)
      mtx_ = Teuchos::null;
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
 
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tExporter.getEpetra_Export(), myDomainMap, myRangeMap, restrictComm));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &RowImporter,
                   const Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > DomainImporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, RowImporter, tImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
    XPETRA_RCP_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, DomainImporter, tdImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tImporter.getEpetra_Import(), tdImporter->getEpetra_Import().get(), myDomainMap, myRangeMap, restrictComm));
    if (restrictComm && mtx_->NumMyRows() == 0)
      mtx_ = Teuchos::null;
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &RowExporter,
                   const Teuchos::RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > DomainExporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, RowExporter, tExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
    XPETRA_RCP_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, DomainExporter, tdExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
 
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tExporter.getEpetra_Export(), tdExporter->getEpetra_Export().get(), myDomainMap, myRangeMap, restrictComm));
  }
 
#ifdef HAVE_XPETRA_TPETRA
  // NOTE: TPETRA means we can use C++11 here
 
  EpetraCrsMatrixT(
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
      const local_matrix_type &lclMatrix,
      const Teuchos::RCP<Teuchos::ParameterList> &params = null)
    : EpetraCrsMatrixT(lclMatrix, rowMap, colMap, Teuchos::null, Teuchos::null, params) {}
 
  EpetraCrsMatrixT(
      const local_matrix_type &lclMatrix,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
      const Teuchos::RCP<Teuchos::ParameterList> &params                           = null) {
    // local typedefs from local_matrix_type
    // typedef typename local_matrix_type::size_type size_type;
    typedef typename local_matrix_type::value_type value_type;
    typedef typename local_matrix_type::ordinal_type ordinal_type;
 
    // The number of rows in the sparse matrix.
    ordinal_type lclNumRows = lclMatrix.numRows();
    ordinal_type lclNumCols = lclMatrix.numCols();  // do we need this?
 
    // plausibility checks
    TEUCHOS_TEST_FOR_EXCEPTION(lclNumRows != Teuchos::as<ordinal_type>(rowMap->getLocalNumElements()), Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrixT: number of rows in local matrix and number of local entries in row map do not match!");
    TEUCHOS_TEST_FOR_EXCEPTION(lclNumCols != Teuchos::as<ordinal_type>(colMap->getLocalNumElements()), Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrixT: number of columns in local matrix and number of local entries in column map do not match!");
 
    Teuchos::ArrayRCP<size_t> NumEntriesPerRowToAlloc(lclNumRows);
    for (ordinal_type r = 0; r < lclNumRows; ++r) {
      // extract data from current row r
      auto rowview               = lclMatrix.row(r);
      NumEntriesPerRowToAlloc[r] = rowview.length;
    }
 
    // setup matrix
    isFillResumed_ = false;
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), numEntriesPerRowToAlloc.getRawPtr(), true));
 
    // loop over all rows and colums of local matrix and fill matrix
    for (ordinal_type r = 0; r < lclNumRows; ++r) {
      // extract data from current row r
      auto rowview = lclMatrix.row(r);
 
      // arrays for current row data
      Teuchos::ArrayRCP<ordinal_type> indout(rowview.length, Teuchos::ScalarTraits<ordinal_type>::zero());
      Teuchos::ArrayRCP<value_type> valout(rowview.length, Teuchos::ScalarTraits<value_type>::zero());
 
      for (ordinal_type c = 0; c < rowview.length; c++) {
        value_type value    = rowview.value(c);
        ordinal_type colidx = rowview.colidx(c);
 
        TEUCHOS_TEST_FOR_EXCEPTION(colMap->isNodeLocalElement(colidx) == false, Xpetra::Exceptions::RuntimeError,
                                   "Xpetra::EpetraCrsMatrixT: local matrix contains column elements which are not in the provided column map!");
 
        indout[c] = colidx;
        valout[c] = value;
      }
      insertLocalValues(r, indout.view(0, indout.size()), valout.view(0, valout.size()));
    }
 
    // call fill complete
    if (!domainMap.is_null() && !rangeMap.is_null())
      this->fillComplete(domainMap, rowMap, params);
    else
      this->fillComplete(rowMap, rowMap, params);
  }
#endif
 
  virtual ~EpetraCrsMatrixT() {}
 
 
 
 
  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
    XPETRA_MONITOR("EpetraCrsMatrixT::insertGlobalValues");
    XPETRA_ERR_CHECK(mtx_->InsertGlobalValues(globalRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
  }
 
  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
    XPETRA_MONITOR("EpetraCrsMatrixT::insertLocalValues");
    XPETRA_ERR_CHECK(mtx_->InsertMyValues(localRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
  }
 
  void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &indices, const ArrayView<const Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceGlobalValues");
 
    {
      const std::string tfecfFuncName("replaceGlobalValues");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!isFillActive(), std::runtime_error,
                                            ": Fill must be active in order to call this method.  If you have already "
                                            "called fillComplete(), you need to call resumeFill() before you can "
                                            "replace values.");
 
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
                                            std::runtime_error, ": values.size() must equal indices.size().");
    }
 
    XPETRA_ERR_CHECK(mtx_->ReplaceGlobalValues(globalRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
  }
 
  void replaceLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &indices, const ArrayView<const Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceLocalValues");
 
    {
      const std::string tfecfFuncName("replaceLocalValues");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!isFillActive(), std::runtime_error,
                                            ": Fill must be active in order to call this method.  If you have already "
                                            "called fillComplete(), you need to call resumeFill() before you can "
                                            "replace values.");
 
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
                                            std::runtime_error, ": values.size() must equal indices.size().");
    }
 
    XPETRA_ERR_CHECK(mtx_->ReplaceMyValues(localRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
  }
 
  void setAllToScalar(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraCrsMatrixT::setAllToScalar");
    mtx_->PutScalar(alpha);
  }
 
  void scale(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraCrsMatrixT::scale");
    mtx_->Scale(alpha);
  }
 
  //** \warning This is an expert-only routine and should not be called from user code. */
  void allocateAllValues(size_t numNonZeros, ArrayRCP<size_t> &rowptr, ArrayRCP<LocalOrdinal> &colind, ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::allocateAllValues");
 
    // Row offsets
    // Unfortunately, we cannot do this in the same manner as column indices
    // and values (see below).  The problem is that Tpetra insists on using
    // size_t, and Epetra uses int internally.  So we only resize here, and
    // will need to copy in setAllValues
    rowptr.resize(getLocalNumRows() + 1);
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    // Column indices
    // Extract, resize, set colind
    Epetra_IntSerialDenseVector &myColind = mtx_->ExpertExtractIndices();
    myColind.Resize(numNonZeros);
    colind = Teuchos::arcp(myColind.Values(), lowerOffset, numNonZeros, ownMemory);
 
    // Values
    // Extract, reallocate, set values
    double *&myValues = mtx_->ExpertExtractValues();
    delete[] myValues;
    myValues = new double[numNonZeros];
    values   = Teuchos::arcp(myValues, lowerOffset, numNonZeros, ownMemory);
  }
 
  void setAllValues(const ArrayRCP<size_t> &rowptr, const ArrayRCP<LocalOrdinal> &colind, const ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::setAllValues");
 
    // Check sizes
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(rowptr.size()) != getLocalNumRows() + 1, Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that the size of your rowptr array is incorrect.");
    TEUCHOS_TEST_FOR_EXCEPTION(values.size() != colind.size(), Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that you mismatched your pointers.");
 
    // Check pointers
    if (values.size() > 0) {
      TEUCHOS_TEST_FOR_EXCEPTION(colind.getRawPtr() != mtx_->ExpertExtractIndices().Values(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
      TEUCHOS_TEST_FOR_EXCEPTION(values.getRawPtr() != mtx_->ExpertExtractValues(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
    }
 
    // We have to make a copy here, it is unavoidable
    // See comments in allocateAllValues
    const size_t N = getLocalNumRows();
 
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    myRowptr.Resize(N + 1);
    for (size_t i = 0; i < N + 1; i++)
      myRowptr[i] = Teuchos::as<int>(rowptr[i]);
  }
 
  void getAllValues(ArrayRCP<const size_t> &rowptr, ArrayRCP<const LocalOrdinal> &colind, ArrayRCP<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getAllValues");
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    const size_t n   = getLocalNumRows();
    const size_t nnz = getLocalNumEntries();
 
    // Row offsets
    // We have to make a copy here, it is unavoidable (see comments in allocateAllValues)
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    rowptr.resize(n + 1);
    for (size_t i = 0; i < n + 1; i++)
      (*const_cast<size_t *>(&rowptr[i])) = Teuchos::as<size_t>(myRowptr[i]);
 
    // Column indices
    colind = Teuchos::arcp(mtx_->ExpertExtractIndices().Values(), lowerOffset, nnz, ownMemory);
 
    // Values
    values = Teuchos::arcp(mtx_->ExpertExtractValues(), lowerOffset, nnz, ownMemory);
  }
 
  void getAllValues(ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::getAllValues");
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    const size_t nnz = getLocalNumEntries();
    // Values
    values = Teuchos::arcp(mtx_->ExpertExtractValues(), lowerOffset, nnz, ownMemory);
  }
 
  // Epetra always has global constants
  bool haveGlobalConstants() const { return true; }
  //** \warning This is an expert-only routine and should not be called from user code. */
  void expertStaticFillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                                const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                                const RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &importer = Teuchos::null,
                                const RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > &exporter = Teuchos::null,
                                const RCP<ParameterList> & /* params */                               = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::expertStaticFillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation.
    isFillResumed_ = false;
 
    int rv                        = 0;
    const Epetra_Import *myimport = 0;
    const Epetra_Export *myexport = 0;
 
    if (!importer.is_null()) {
      XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, *importer, eImporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
      myimport = eImporter.getEpetra_Import().getRawPtr();
    }
    if (!exporter.is_null()) {
      XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, *exporter, eExporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
      myexport = eExporter.getEpetra_Export().getRawPtr();
    }
 
    rv = mtx_->ExpertStaticFillComplete(toEpetra<GlobalOrdinal, Node>(domainMap), toEpetra<GlobalOrdinal, Node>(rangeMap), myimport, myexport);
 
    TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete FAILED!");
  }
 
 
 
  void resumeFill(const RCP<ParameterList> & /* params */ = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::resumeFill");
 
    // According to Tpetra documentation, resumeFill() may be called repeatedly.
    isFillResumed_ = true;
  }
 
  void fillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap, const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap, const RCP<ParameterList> &params = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation.
    isFillResumed_ = false;
 
    bool doOptimizeStorage = true;
    if (params != null && params->get("Optimize Storage", true) == false) doOptimizeStorage = false;
    mtx_->FillComplete(toEpetra<GlobalOrdinal, Node>(domainMap), toEpetra<GlobalOrdinal, Node>(rangeMap), doOptimizeStorage);
  }
 
  void fillComplete(const RCP<ParameterList> &params = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation. There is no need for a fillComplete after some resumeFill() operations.
    if (isFillResumed_ == true) {
      isFillResumed_ = false;
      return;
    }
 
    bool doOptimizeStorage = true;
    if (params != null && params->get("Optimize Storage", true) == false) doOptimizeStorage = false;
    mtx_->FillComplete(doOptimizeStorage);
  }
 
  void replaceDomainMapAndImporter(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newDomainMap, Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &newImporter) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceDomainMapAndImporter");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, *newImporter, eImporter, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter only accepts Xpetra::EpetraImportT.");
 
    const RCP<const Epetra_Import> &myImport = eImporter.getEpetra_Import();
    int rv                                   = 0;
    if (myImport == Teuchos::null)
      rv = mtx_->ReplaceDomainMapAndImporter(toEpetra<GlobalOrdinal, Node>(newDomainMap), 0);
    else
      rv = mtx_->ReplaceDomainMapAndImporter(toEpetra<GlobalOrdinal, Node>(newDomainMap), &*myImport);
    TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter FAILED!");
  }
 
 
 
 
  const RCP<const Comm<int> > getComm() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getComm");
    return toXpetra(mtx_->Comm());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRowMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getRowMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->RowMap());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getColMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getColMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->ColMap());
  }
 
  RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > getCrsGraph() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getCrsGraph");
    return toXpetra<GlobalOrdinal, Node>(mtx_->Graph());
  }
 
  global_size_t getGlobalNumRows() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumRows");
    return mtx_->NumGlobalRows64();
  }
 
  global_size_t getGlobalNumCols() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumCols");
    return mtx_->NumGlobalCols64();
  }
 
  size_t getLocalNumRows() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumRows");
    return mtx_->NumMyRows();
  }
 
  size_t getLocalNumCols() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumCols");
    return mtx_->NumMyCols();
  }
 
  global_size_t getGlobalNumEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumEntries");
    return mtx_->NumGlobalNonzeros64();
  }
 
  size_t getLocalNumEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumEntries");
    return mtx_->NumMyNonzeros();
  }
 
  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getNumEntriesInLocalRow");
    return mtx_->NumMyEntries(localRow);
  }
 
  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getNumEntriesInGlobalRow");
    return mtx_->NumGlobalEntries(globalRow);
  }
 
  size_t getGlobalMaxNumRowEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalMaxNumRowEntries");
    return mtx_->GlobalMaxNumEntries();
  }
 
  size_t getLocalMaxNumRowEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalMaxNumRowEntries");
    return mtx_->MaxNumEntries();
  }
 
  bool isLocallyIndexed() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isLocallyIndexed");
    return mtx_->IndicesAreLocal();
  }
 
  bool isGloballyIndexed() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isGloballyIndexed");
    return mtx_->IndicesAreGlobal();
  }
 
  bool isFillComplete() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isFillComplete");
    if (isFillResumed_)
      return false;
    else
      return mtx_->Filled();
  }
 
  bool isFillActive() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isFillActive");
    return !isFillComplete();
  }
 
  ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getFrobeniusNorm");
    return mtx_->NormFrobenius();
  }
 
  bool supportsRowViews() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::supportsRowViews");
    return true;
  }
 
  void getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView<LocalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowCopy");
 
    int numEntries = -1;
    XPETRA_ERR_CHECK(mtx_->ExtractMyRowCopy(LocalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
    NumEntries = numEntries;
  }
 
  void getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView<GlobalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowCopy");
 
    int numEntries = -1;
    XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowCopy(GlobalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
    NumEntries = numEntries;
  }
 
  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowView");
 
    int numEntries;
    double *eValues;
    GlobalOrdinal *eIndices;
 
    XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowView(GlobalRow, numEntries, eValues, eIndices));
    if (numEntries == 0) {
      eValues  = NULL;
      eIndices = NULL;
    }  // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
 
    indices = ArrayView<const GlobalOrdinal>(eIndices, numEntries);
    values  = ArrayView<const double>(eValues, numEntries);
  }
 
  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowView");
 
    int numEntries;
    double *eValues;
    int *eIndices;
 
    XPETRA_ERR_CHECK(mtx_->ExtractMyRowView(LocalRow, numEntries, eValues, eIndices));
    if (numEntries == 0) {
      eValues  = NULL;
      eIndices = NULL;
    }  // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
 
    indices = ArrayView<const int>(eIndices, numEntries);
    values  = ArrayView<const double>(eValues, numEntries);
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalDiagCopy");
    XPETRA_ERR_CHECK(mtx_->ExtractDiagonalCopy(toEpetra<GlobalOrdinal, Node>(diag)));
  }
 
  void getLocalDiagOffsets(Teuchos::ArrayRCP<size_t> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagOffsets() is not implemented or supported.");
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & /* diag */, const Teuchos::ArrayView<const size_t> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagCopy using offsets is not implemented or supported.");
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & /* diag */, const Kokkos::View<const size_t *, typename Node::device_type, Kokkos::MemoryUnmanaged> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagCopy using offsets is not implemented or supported.");
  }
 
  void replaceDiag(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) {
    mtx_->ReplaceDiagonalValues(toEpetra<GlobalOrdinal, Node>(diag));
  }
 
  void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x) {
    XPETRA_ERR_CHECK(mtx_->LeftScale(toEpetra<GlobalOrdinal, Node>(x)));
  };
  void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x) {
    XPETRA_ERR_CHECK(mtx_->RightScale(toEpetra<GlobalOrdinal, Node>(x)));
  };
 
 
 
 
  void apply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X, MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &Y, Teuchos::ETransp mode = Teuchos::NO_TRANS, Scalar alpha = ScalarTraits<Scalar>::one(), Scalar beta = ScalarTraits<Scalar>::zero()) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::apply");
 
    // TEUCHOS_TEST_FOR_EXCEPTION((alpha != 1) || (beta != 0), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.multiply() only accept alpha==1 and beta==0");
 
    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, X, eX, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, Y, eY, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
 
    TEUCHOS_TEST_FOR_EXCEPTION((mode != Teuchos::NO_TRANS) && (mode != Teuchos::TRANS), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT->apply() only accept mode == NO_TRANS or mode == TRANS");
    bool eTrans = toEpetra(mode);
 
    // /!\ UseTranspose value
    TEUCHOS_TEST_FOR_EXCEPTION(mtx_->UseTranspose(), Xpetra::Exceptions::NotImplemented, "An exception is throw to let you know that Xpetra::EpetraCrsMatrixT->apply() do not take into account the UseTranspose() parameter of Epetra_CrsMatrix.");
 
    RCP<Epetra_MultiVector> epY = eY.getEpetra_MultiVector();
 
    // helper vector: tmp = A*x
    RCP<Epetra_MultiVector> tmp = Teuchos::rcp(new Epetra_MultiVector(*epY));
    tmp->PutScalar(0.0);
    XPETRA_ERR_CHECK(mtx_->Multiply(eTrans, *eX.getEpetra_MultiVector(), *tmp));
 
    // calculate alpha * A * x + beta * y
    XPETRA_ERR_CHECK(eY.getEpetra_MultiVector()->Update(alpha, *tmp, beta));
  }
 
  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 {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.apply() is not implemented for region matricies.");
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getDomainMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->DomainMap());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getRangeMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->RangeMap());
  }
 
 
 
 
  std::string description() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::description");
 
    // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
    std::ostringstream oss;
    // TODO: oss << DistObject<char, LocalOrdinal,GlobalOrdinal>::description();
    if (isFillComplete()) {
      oss << "{status = fill complete"
          << ", global rows = " << getGlobalNumRows()
          << ", global cols = " << getGlobalNumCols()
          << ", global num entries = " << getGlobalNumEntries()
          << "}";
    } else {
      oss << "{status = fill not complete"
          << ", global rows = " << getGlobalNumRows()
          << "}";
    }
    return oss.str();
  }
 
  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::describe");
 
    // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
    using std::endl;
    using std::setw;
    using Teuchos::VERB_DEFAULT;
    using Teuchos::VERB_EXTREME;
    using Teuchos::VERB_HIGH;
    using Teuchos::VERB_LOW;
    using Teuchos::VERB_MEDIUM;
    using Teuchos::VERB_NONE;
    Teuchos::EVerbosityLevel vl = verbLevel;
    if (vl == VERB_DEFAULT) vl = VERB_LOW;
    RCP<const Comm<int> > comm = this->getComm();
    const int myImageID        = comm->getRank(),
              numImages        = comm->getSize();
    size_t width               = 1;
    for (size_t dec = 10; dec < getGlobalNumRows(); dec *= 10) {
      ++width;
    }
    width = std::max<size_t>(width, 11) + 2;
    Teuchos::OSTab tab(out);
    //    none: print nothing
    //     low: print O(1) info from node 0
    //  medium: print O(P) info, num entries per node
    //    high: print O(N) info, num entries per row
    // extreme: print O(NNZ) info: print indices and values
    //
    // for medium and higher, print constituent objects at specified verbLevel
    if (vl != VERB_NONE) {
      if (myImageID == 0) out << this->description() << std::endl;
      // O(1) globals, minus what was already printed by description()
      if (isFillComplete() && myImageID == 0) {
        out << "Global max number of entries = " << getGlobalMaxNumRowEntries() << std::endl;
      }
      // constituent objects
      if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
        if (myImageID == 0) out << "\nRow map: " << std::endl;
        getRowMap()->describe(out, vl);
        //
        if (getColMap() != null) {
          if (getColMap() == getRowMap()) {
            if (myImageID == 0) out << "\nColumn map is row map.";
          } else {
            if (myImageID == 0) out << "\nColumn map: " << std::endl;
            getColMap()->describe(out, vl);
          }
        }
        if (getDomainMap() != null) {
          if (getDomainMap() == getRowMap()) {
            if (myImageID == 0) out << "\nDomain map is row map.";
          } else if (getDomainMap() == getColMap()) {
            if (myImageID == 0) out << "\nDomain map is row map.";
          } else {
            if (myImageID == 0) out << "\nDomain map: " << std::endl;
            getDomainMap()->describe(out, vl);
          }
        }
        if (getRangeMap() != null) {
          if (getRangeMap() == getDomainMap()) {
            if (myImageID == 0) out << "\nRange map is domain map." << std::endl;
          } else if (getRangeMap() == getRowMap()) {
            if (myImageID == 0) out << "\nRange map is row map." << std::endl;
          } else {
            if (myImageID == 0) out << "\nRange map: " << std::endl;
            getRangeMap()->describe(out, vl);
          }
        }
        if (myImageID == 0) out << std::endl;
      }
      // O(P) data
      if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
        for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
          if (myImageID == imageCtr) {
            out << "Node ID = " << imageCtr << std::endl;
            // TODO: need a graph
            //               if (staticGraph_->indicesAreAllocated() == false) {
            //                 out << "Node not allocated" << std::endl;
            //               }
            //               else {
            //                 out << "Node number of allocated entries = " << staticGraph_->getLocalAllocationSize() << std::endl;
            //               }
 
            // TMP:
            //            const Epetra_CrsGraph & staticGraph_ = mtx_->Graph();
            // End of TMP
 
            out << "Node number of entries = " << getLocalNumEntries() << std::endl;
            out << "Node max number of entries = " << getLocalMaxNumRowEntries() << std::endl;
          }
          comm->barrier();
          comm->barrier();
          comm->barrier();
        }
      }
      // O(N) and O(NNZ) data
      if (vl == VERB_HIGH || vl == VERB_EXTREME) {
        for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
          if (myImageID == imageCtr) {
            out << std::setw(width) << "Node ID"
                << std::setw(width) << "Global Row"
                << std::setw(width) << "Num Entries";
            if (vl == VERB_EXTREME) {
              out << std::setw(width) << "(Index,Value)";
            }
            out << std::endl;
            for (size_t r = 0; r < getLocalNumRows(); ++r) {
              const size_t nE   = getNumEntriesInLocalRow(r);
              GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
              out << std::setw(width) << myImageID
                  << std::setw(width) << gid
                  << std::setw(width) << nE;
              if (vl == VERB_EXTREME) {
                if (isGloballyIndexed()) {
                  ArrayView<const GlobalOrdinal> rowinds;
                  ArrayView<const Scalar> rowvals;
                  getGlobalRowView(gid, rowinds, rowvals);
                  for (size_t j = 0; j < nE; ++j) {
                    out << " (" << rowinds[j]
                        << ", " << rowvals[j]
                        << ") ";
                  }
                } else if (isLocallyIndexed()) {
                  ArrayView<const LocalOrdinal> rowinds;
                  ArrayView<const Scalar> rowvals;
                  getLocalRowView(r, rowinds, rowvals);
                  for (size_t j = 0; j < nE; ++j) {
                    out << " (" << getColMap()->getGlobalElement(rowinds[j])
                        << ", " << rowvals[j]
                        << ") ";
                  }
                }
              }
              out << std::endl;
            }
          }
          comm->barrier();
          comm->barrier();
          comm->barrier();
        }
      }
    }
  }
 
  void setObjectLabel(const std::string &objectLabel) {
    Teuchos::LabeledObject::setObjectLabel(objectLabel);
    mtx_->SetLabel(objectLabel.c_str());
  }
 
  EpetraCrsMatrixT(const EpetraCrsMatrixT &matrix)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(*(matrix.mtx_))))
    , isFillResumed_(false) {}
 
  //{@
 
  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->Map());
  }
 
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
                const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
    int err                       = mtx_->Import(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
                const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
    int err                       = mtx_->Export(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
                const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
    int err                       = mtx_->Import(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
                const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
    int err                       = mtx_->Export(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Caught nonzero error code " + std::to_string(err) + " returned by Epetra.");
  }
 
  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newMap) {
    XPETRA_MONITOR("EpetraCrsMatrixT::removeEmptyProcessesInPlace");
    const Epetra_Map *newMapEpetra = (!newMap.is_null()) ? &toEpetra<GlobalOrdinal, Node>(newMap) : 0;
    int err                        = mtx_->RemoveEmptyProcessesInPlace(newMapEpetra);
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Caught nonzero error code " + std::to_string(err) + " returned by Epetra.");
  }
 
 
 
 
  bool hasMatrix() const { return !mtx_.is_null(); }
 
  EpetraCrsMatrixT(const Teuchos::RCP<Epetra_CrsMatrix> &mtx)
    : mtx_(mtx)
    , isFillResumed_(false) {}
 
  RCP<const Epetra_CrsMatrix> getEpetra_CrsMatrix() const { return mtx_; }
 
  RCP<Epetra_CrsMatrix> getEpetra_CrsMatrixNonConst() const { return mtx_; }  // TODO: remove
 
#ifdef HAVE_XPETRA_TPETRA
 
  local_matrix_type getLocalMatrixDevice() const {
#if 0
  TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented,
           "Xpetra::EpetraCrsMatrx only available on host for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  TEUCHOS_UNREACHABLE_RETURN((local_matrix_type()));
#endif
    return getLocalMatrixHost();
  }
 
#if KOKKOS_VERSION >= 40799
  typename local_matrix_type::host_mirror_type getLocalMatrixHost() const {
#else
  typename local_matrix_type::HostMirror getLocalMatrixHost() const {
#endif
    RCP<Epetra_CrsMatrix> matrix = getEpetra_CrsMatrixNonConst();
 
    const int numRows = matrix->NumMyRows();
    const int numCols = matrix->NumMyCols();
    const int nnz     = matrix->NumMyNonzeros();
 
    int *rowptr;
    int *colind;
    double *vals;
    int rv = matrix->ExtractCrsDataPointers(rowptr, colind, vals);
    TEUCHOS_TEST_FOR_EXCEPTION(rv, std::runtime_error, "Xpetra::CrsMatrix<>::getLocalMatrix: failed in ExtractCrsDataPointers");
 
    // Transform int* rowptr array to size_type* array
    typename local_matrix_type::row_map_type::non_const_type kokkosRowPtr(Kokkos::ViewAllocateWithoutInitializing("local row map"), numRows + 1);
    for (size_t i = 0; i < kokkosRowPtr.size(); i++)
      kokkosRowPtr(i) = Teuchos::asSafe<typename local_matrix_type::row_map_type::value_type>(rowptr[i]);
 
    // create Kokkos::Views
    typename local_matrix_type::index_type kokkosColind(colind, nnz);
    typename local_matrix_type::values_type kokkosVals(vals, nnz);
 
    local_matrix_type localMatrix = local_matrix_type("LocalMatrix", numRows, numCols, nnz, kokkosVals, kokkosRowPtr, kokkosColind);
 
    return localMatrix;
  }
 
  void setAllValues(const typename local_matrix_type::row_map_type &ptr,
                    const typename local_matrix_type::StaticCrsGraphType::entries_type::non_const_type &ind,
                    const typename local_matrix_type::values_type &val) {
    // Check sizes
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ptr.size()) != getLocalNumRows() + 1, Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that the size of your rowptr array is incorrect.");
    TEUCHOS_TEST_FOR_EXCEPTION(val.size() != ind.size(), Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that you mismatched your pointers.");
 
    // Check pointers
    if (val.size() > 0) {
      std::cout << ind.data() << " " << mtx_->ExpertExtractIndices().Values() << std::endl;
      TEUCHOS_TEST_FOR_EXCEPTION(ind.data() != mtx_->ExpertExtractIndices().Values(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
      TEUCHOS_TEST_FOR_EXCEPTION(val.data() != mtx_->ExpertExtractValues(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
    }
 
    // We have to make a copy here, it is unavoidable
    // See comments in allocateAllValues
    const size_t N = getLocalNumRows();
 
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    myRowptr.Resize(N + 1);
    for (size_t i = 0; i < N + 1; i++)
      myRowptr[i] = Teuchos::as<int>(ptr(i));
  }
 
  LocalOrdinal GetStorageBlockSize() const { return 1; }
 
 private:
#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
 
  void residual(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X,
                const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B,
                MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &R) const {
    Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one;
 
    apply(X, R);
    R.update(one, B, negone);
  }
 
 private:
  RCP<Epetra_CrsMatrix> mtx_;
 
  bool isFillResumed_;
 
};  // EpetraCrsMatrixT class
 
#endif  //#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
 
#ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
template <>
class EpetraCrsMatrixT<long long, EpetraNode>
  : public CrsMatrix<double, int, long long, EpetraNode> {
  typedef long long GlobalOrdinal;
  typedef EpetraNode Node;
  typedef typename CrsMatrix<double, int, GlobalOrdinal, Node>::scalar_type Scalar;
  typedef typename CrsMatrix<double, int, GlobalOrdinal, Node>::local_ordinal_type LocalOrdinal;
 
  // The following typedefs are used by the Kokkos interface
#ifdef HAVE_XPETRA_TPETRA
 public:
  typedef typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_type 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;
  typedef typename Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::node_type node_type;
 
 private:
#endif
 
 public:
 
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), maxNumEntriesPerRow, true)))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : isFillResumed_(false) {
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), numEntriesPerRowToAlloc.getRawPtr(), true));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, size_t maxNumEntriesPerRow, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), maxNumEntriesPerRow, true)))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap, const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : isFillResumed_(false) {
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), numEntriesPerRowToAlloc.getRawPtr(), true));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > &graph, const Teuchos::RCP<Teuchos::ParameterList> & /* plist */ = Teuchos::null)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(graph))))
    , isFillResumed_(false) {}
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &importer,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tImporter.getEpetra_Import(), myDomainMap, myRangeMap, restrictComm));
    if (restrictComm && mtx_->NumMyRows() == 0)
      mtx_ = Teuchos::null;
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
                   const Teuchos::RCP<Teuchos::ParameterList> &params                           = Teuchos::null)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
 
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tExporter.getEpetra_Export(), myDomainMap, myRangeMap, restrictComm));
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Import<LocalOrdinal, GlobalOrdinal, Node> &RowImporter,
                   const Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > DomainImporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, RowImporter, tImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
    XPETRA_RCP_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, DomainImporter, tdImporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraImportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tImporter.getEpetra_Import(), tdImporter->getEpetra_Import().get(), myDomainMap, myRangeMap, restrictComm));
    if (restrictComm && mtx_->NumMyRows() == 0)
      mtx_ = Teuchos::null;
  }
 
  EpetraCrsMatrixT(const Teuchos::RCP<const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > &sourceMatrix,
                   const Export<LocalOrdinal, GlobalOrdinal, Node> &RowExporter,
                   const Teuchos::RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > DomainExporter,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                   const Teuchos::RCP<Teuchos::ParameterList> &params)
    : isFillResumed_(false) {
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, *sourceMatrix, tSourceMatrix, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraCrsMatrixT as an input argument.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, RowExporter, tExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
    XPETRA_RCP_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, DomainExporter, tdExporter, "Xpetra::EpetraCrsMatrixT constructor only accepts Xpetra::EpetraExportT as an input argument.");
 
    const Epetra_Map *myDomainMap = (domainMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(domainMap) : 0;
    const Epetra_Map *myRangeMap  = (rangeMap != Teuchos::null) ? &toEpetra<GlobalOrdinal, Node>(rangeMap) : 0;
 
    // Follows the Tpetra parameters
    bool restrictComm = false;
    if (!params.is_null()) restrictComm = params->get("Restrict Communicator", restrictComm);
 
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(*tSourceMatrix.getEpetra_CrsMatrix(), *tExporter.getEpetra_Export(), tdExporter->getEpetra_Export().get(), myDomainMap, myRangeMap, restrictComm));
  }
 
#ifdef HAVE_XPETRA_TPETRA
  EpetraCrsMatrixT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
                   const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
                   const local_matrix_type &lclMatrix,
                   const Teuchos::RCP<Teuchos::ParameterList> &params = null)
    : EpetraCrsMatrixT(lclMatrix, rowMap, colMap, Teuchos::null, Teuchos::null, params) {}
 
  EpetraCrsMatrixT(
      const local_matrix_type &lclMatrix,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rowMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &colMap,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap = Teuchos::null,
      const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap  = Teuchos::null,
      const Teuchos::RCP<Teuchos::ParameterList> &params                           = null) {
    // local typedefs from local_matrix_type
    // typedef typename local_matrix_type::size_type size_type;
    typedef typename local_matrix_type::value_type value_type;
    typedef typename local_matrix_type::ordinal_type ordinal_type;
 
    // The number of rows in the sparse matrix.
    ordinal_type lclNumRows = lclMatrix.numRows();
    ordinal_type lclNumCols = lclMatrix.numCols();  // do we need this?
 
    // plausibility checks
    TEUCHOS_TEST_FOR_EXCEPTION(lclNumRows != Teuchos::as<ordinal_type>(rowMap->getLocalNumElements()), Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrixT: number of rows in local matrix and number of local entries in row map do not match!");
    TEUCHOS_TEST_FOR_EXCEPTION(lclNumCols != Teuchos::as<ordinal_type>(colMap->getLocalNumElements()), Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraCrsMatrixT: number of columns in local matrix and number of local entries in column map do not match!");
 
    Teuchos::ArrayRCP<size_t> NumEntriesPerRowToAlloc(lclNumRows);
    for (ordinal_type r = 0; r < lclNumRows; ++r) {
      // extract data from current row r
      auto rowview               = lclMatrix.row(r);
      NumEntriesPerRowToAlloc[r] = rowview.length;
    }
 
    // setup matrix
    isFillResumed_ = false;
    Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end());  // convert array of "size_t" to array of "int"
    mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra<GlobalOrdinal, Node>(rowMap), toEpetra<GlobalOrdinal, Node>(colMap), numEntriesPerRowToAlloc.getRawPtr(), true));
 
    // loop over all rows and colums of local matrix and fill matrix
    for (ordinal_type r = 0; r < lclNumRows; ++r) {
      // extract data from current row r
      auto rowview = lclMatrix.row(r);
 
      // arrays for current row data
      Teuchos::ArrayRCP<ordinal_type> indout(rowview.length, Teuchos::ScalarTraits<ordinal_type>::zero());
      Teuchos::ArrayRCP<value_type> valout(rowview.length, Teuchos::ScalarTraits<value_type>::zero());
 
      for (ordinal_type c = 0; c < rowview.length; c++) {
        value_type value    = rowview.value(c);
        ordinal_type colidx = rowview.colidx(c);
 
        TEUCHOS_TEST_FOR_EXCEPTION(colMap->isNodeLocalElement(colidx) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::EpetraCrsMatrixT: local matrix contains column elements which are not in the provided column map!");
 
        indout[c] = colidx;
        valout[c] = value;
      }
      insertLocalValues(r, indout.view(0, indout.size()), valout.view(0, valout.size()));
    }
 
    // call fill complete
    if (!domainMap.is_null() && !rangeMap.is_null())
      this->fillComplete(domainMap, rowMap, params);
    else
      this->fillComplete(rowMap, rowMap, params);
  }
#endif
 
  virtual ~EpetraCrsMatrixT() {}
 
 
 
 
  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
    XPETRA_MONITOR("EpetraCrsMatrixT::insertGlobalValues");
    XPETRA_ERR_CHECK(mtx_->InsertGlobalValues(globalRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
  }
 
  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals) {
    XPETRA_MONITOR("EpetraCrsMatrixT::insertLocalValues");
    XPETRA_ERR_CHECK(mtx_->InsertMyValues(localRow, vals.size(), vals.getRawPtr(), cols.getRawPtr()));
  }
 
  void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &indices, const ArrayView<const Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceGlobalValues");
 
    {
      const std::string tfecfFuncName("replaceGlobalValues");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!isFillActive(), std::runtime_error,
                                            ": Fill must be active in order to call this method.  If you have already "
                                            "called fillComplete(), you need to call resumeFill() before you can "
                                            "replace values.");
 
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
                                            std::runtime_error, ": values.size() must equal indices.size().");
    }
 
    XPETRA_ERR_CHECK(mtx_->ReplaceGlobalValues(globalRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
  }
 
  void replaceLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &indices, const ArrayView<const Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceLocalValues");
 
    {
      const std::string tfecfFuncName("replaceLocalValues");
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!isFillActive(), std::runtime_error,
                                            ": Fill must be active in order to call this method.  If you have already "
                                            "called fillComplete(), you need to call resumeFill() before you can "
                                            "replace values.");
 
      TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
                                            std::runtime_error, ": values.size() must equal indices.size().");
    }
 
    XPETRA_ERR_CHECK(mtx_->ReplaceMyValues(localRow, indices.size(), values.getRawPtr(), indices.getRawPtr()));
  }
 
  void setAllToScalar(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraCrsMatrixT::setAllToScalar");
    mtx_->PutScalar(alpha);
  }
 
  void scale(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraCrsMatrixT::scale");
    mtx_->Scale(alpha);
  }
 
  //** \warning This is an expert-only routine and should not be called from user code. */
  void allocateAllValues(size_t numNonZeros, ArrayRCP<size_t> &rowptr, ArrayRCP<LocalOrdinal> &colind, ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::allocateAllValues");
 
    // Row offsets
    // Unfortunately, we cannot do this in the same manner as column indices
    // and values (see below).  The problem is that Tpetra insists on using
    // size_t, and Epetra uses int internally.  So we only resize here, and
    // will need to copy in setAllValues
    rowptr.resize(getLocalNumRows() + 1);
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    // Column indices
    // Extract, resize, set colind
    Epetra_IntSerialDenseVector &myColind = mtx_->ExpertExtractIndices();
    myColind.Resize(numNonZeros);
    colind = Teuchos::arcp(myColind.Values(), lowerOffset, numNonZeros, ownMemory);
 
    // Values
    // Extract, reallocate, set values
    double *&myValues = mtx_->ExpertExtractValues();
    delete[] myValues;
    myValues = new double[numNonZeros];
    values   = Teuchos::arcp(myValues, lowerOffset, numNonZeros, ownMemory);
  }
 
  void setAllValues(const ArrayRCP<size_t> &rowptr, const ArrayRCP<LocalOrdinal> &colind, const ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::setAllValues");
 
    // Check sizes
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(rowptr.size()) != getLocalNumRows() + 1, Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that the size of your rowptr array is incorrect.");
    TEUCHOS_TEST_FOR_EXCEPTION(values.size() != colind.size(), Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that you mismatched your pointers.");
 
    // Check pointers
    if (values.size() > 0) {
      TEUCHOS_TEST_FOR_EXCEPTION(colind.getRawPtr() != mtx_->ExpertExtractIndices().Values(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
      TEUCHOS_TEST_FOR_EXCEPTION(values.getRawPtr() != mtx_->ExpertExtractValues(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
    }
 
    // We have to make a copy here, it is unavoidable
    // See comments in allocateAllValues
    const size_t N = getLocalNumRows();
 
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    myRowptr.Resize(N + 1);
    for (size_t i = 0; i < N + 1; i++)
      myRowptr[i] = Teuchos::as<int>(rowptr[i]);
  }
 
  void getAllValues(ArrayRCP<const size_t> &rowptr, ArrayRCP<const LocalOrdinal> &colind, ArrayRCP<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getAllValues");
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    const size_t n   = getLocalNumRows();
    const size_t nnz = getLocalNumEntries();
 
    // Row offsets
    // We have to make a copy here, it is unavoidable (see comments in allocateAllValues)
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    rowptr.resize(n + 1);
    for (size_t i = 0; i < n + 1; i++)
      (*const_cast<size_t *>(&rowptr[i])) = Teuchos::as<size_t>(myRowptr[i]);
 
    // Column indices
    colind = Teuchos::arcp(mtx_->ExpertExtractIndices().Values(), lowerOffset, nnz, ownMemory);
 
    // Values
    values = Teuchos::arcp(mtx_->ExpertExtractValues(), lowerOffset, nnz, ownMemory);
  }
 
  void getAllValues(ArrayRCP<Scalar> &values) {
    XPETRA_MONITOR("EpetraCrsMatrixT::getAllValues");
 
    int lowerOffset = 0;
    bool ownMemory  = false;
 
    const size_t nnz = getLocalNumEntries();
    // Values
    values = Teuchos::arcp(mtx_->ExpertExtractValues(), lowerOffset, nnz, ownMemory);
  }
 
  // Epetra always has global constants
  bool haveGlobalConstants() const { return true; }
 
  //** \warning This is an expert-only routine and should not be called from user code. */
  void expertStaticFillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap,
                                const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap,
                                const RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &importer = Teuchos::null,
                                const RCP<const Export<LocalOrdinal, GlobalOrdinal, Node> > &exporter = Teuchos::null,
                                const RCP<ParameterList> & /* params */                               = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::expertStaticFillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation.
    isFillResumed_ = false;
 
    int rv                        = 0;
    const Epetra_Import *myimport = 0;
    const Epetra_Export *myexport = 0;
 
    if (!importer.is_null()) {
      XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, *importer, eImporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
      myimport = eImporter.getEpetra_Import().getRawPtr();
    }
    if (!exporter.is_null()) {
      XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, *exporter, eExporter, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete only accepts Xpetra::EpetraImportT.");
      myexport = eExporter.getEpetra_Export().getRawPtr();
    }
 
    rv = mtx_->ExpertStaticFillComplete(toEpetra<GlobalOrdinal, Node>(domainMap), toEpetra<GlobalOrdinal, Node>(rangeMap), myimport, myexport);
 
    TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::expertStaticFillComplete FAILED!");
  }
 
 
 
  void resumeFill(const RCP<ParameterList> & /* params */ = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::resumeFill");
 
    // According to Tpetra documentation, resumeFill() may be called repeatedly.
    isFillResumed_ = true;
  }
 
  void fillComplete(const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &domainMap, const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &rangeMap, const RCP<ParameterList> &params = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation. There is no need for a fillComplete after some resumeFill() operations.
    if (isFillResumed_ == true) {
      isFillResumed_ = false;
      return;
    }
 
    bool doOptimizeStorage = true;
    if (params != null && params->get("Optimize Storage", true) == false) doOptimizeStorage = false;
    mtx_->FillComplete(toEpetra<GlobalOrdinal, Node>(domainMap), toEpetra<GlobalOrdinal, Node>(rangeMap), doOptimizeStorage);
  }
 
  void fillComplete(const RCP<ParameterList> &params = Teuchos::null) {
    XPETRA_MONITOR("EpetraCrsMatrixT::fillComplete");
 
    // For Epetra matrices, resumeFill() is a fictive operation. There is no need for a fillComplete after some resumeFill() operations.
    if (isFillResumed_ == true) {
      isFillResumed_ = false;
      return;
    }
 
    bool doOptimizeStorage = true;
    if (params != null && params->get("Optimize Storage", true) == false) doOptimizeStorage = false;
    mtx_->FillComplete(doOptimizeStorage);
  }
 
  void replaceDomainMapAndImporter(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newDomainMap, Teuchos::RCP<const Import<LocalOrdinal, GlobalOrdinal, Node> > &newImporter) {
    XPETRA_MONITOR("EpetraCrsMatrixT::replaceDomainMapAndImporter");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, *newImporter, eImporter, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter only accepts Xpetra::EpetraImportT.");
 
    const RCP<const Epetra_Import> &myImport = eImporter.getEpetra_Import();
    int rv                                   = 0;
    if (myImport == Teuchos::null)
      rv = mtx_->ReplaceDomainMapAndImporter(toEpetra<GlobalOrdinal, Node>(newDomainMap), 0);
    else
      rv = mtx_->ReplaceDomainMapAndImporter(toEpetra<GlobalOrdinal, Node>(newDomainMap), &*myImport);
    TEUCHOS_TEST_FOR_EXCEPTION(rv != 0, std::runtime_error, "Xpetra::EpetraCrsMatrixT::replaceDomainMapAndImporter FAILED!");
  }
 
 
 
 
  const RCP<const Comm<int> > getComm() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getComm");
    return toXpetra(mtx_->Comm());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRowMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getRowMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->RowMap());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getColMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getColMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->ColMap());
  }
 
  RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> > getCrsGraph() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getCrsGraph");
    return toXpetra<GlobalOrdinal, Node>(mtx_->Graph());
  }
 
  global_size_t getGlobalNumRows() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumRows");
    return mtx_->NumGlobalRows64();
  }
 
  global_size_t getGlobalNumCols() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumCols");
    return mtx_->NumGlobalCols64();
  }
 
  size_t getLocalNumRows() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumRows");
    return mtx_->NumMyRows();
  }
 
  size_t getLocalNumCols() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumCols");
    return mtx_->NumMyCols();
  }
 
  global_size_t getGlobalNumEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalNumEntries");
    return mtx_->NumGlobalNonzeros64();
  }
 
  size_t getLocalNumEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalNumEntries");
    return mtx_->NumMyNonzeros();
  }
 
  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getNumEntriesInLocalRow");
    return mtx_->NumMyEntries(localRow);
  }
 
  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getNumEntriesInGlobalRow");
    return mtx_->NumGlobalEntries(globalRow);
  }
 
  size_t getGlobalMaxNumRowEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalMaxNumRowEntries");
    return mtx_->GlobalMaxNumEntries();
  }
 
  size_t getLocalMaxNumRowEntries() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalMaxNumRowEntries");
    return mtx_->MaxNumEntries();
  }
 
  bool isLocallyIndexed() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isLocallyIndexed");
    return mtx_->IndicesAreLocal();
  }
 
  bool isGloballyIndexed() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isGloballyIndexed");
    return mtx_->IndicesAreGlobal();
  }
 
  bool isFillComplete() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isFillComplete");
    if (isFillResumed_)
      return false;
    else
      return mtx_->Filled();
  }
 
  bool isFillActive() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::isFillActive");
    return !isFillComplete();
  }
 
  typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getFrobeniusNorm");
    return mtx_->NormFrobenius();
  }
 
  bool supportsRowViews() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::supportsRowViews");
    return true;
  }
 
  void getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView<LocalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowCopy");
 
    int numEntries = -1;
    XPETRA_ERR_CHECK(mtx_->ExtractMyRowCopy(LocalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
    NumEntries = numEntries;
  }
 
  void getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView<GlobalOrdinal> &Indices, const ArrayView<Scalar> &Values, size_t &NumEntries) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowCopy");
 
    int numEntries = -1;
    XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowCopy(GlobalRow, Indices.size(), numEntries, Values.getRawPtr(), Indices.getRawPtr()));
    NumEntries = numEntries;
  }
 
  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getGlobalRowView");
 
    int numEntries;
    double *eValues;
    GlobalOrdinal *eIndices;
 
    XPETRA_ERR_CHECK(mtx_->ExtractGlobalRowView(GlobalRow, numEntries, eValues, eIndices));
    if (numEntries == 0) {
      eValues  = NULL;
      eIndices = NULL;
    }  // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
 
    indices = ArrayView<const GlobalOrdinal>(eIndices, numEntries);
    values  = ArrayView<const double>(eValues, numEntries);
  }
 
  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalRowView");
 
    int numEntries;
    double *eValues;
    int *eIndices;
 
    XPETRA_ERR_CHECK(mtx_->ExtractMyRowView(LocalRow, numEntries, eValues, eIndices));
    if (numEntries == 0) {
      eValues  = NULL;
      eIndices = NULL;
    }  // Cf. TEUCHOS_TEST_FOR_EXCEPT( p == 0 && size_in != 0 ) in Teuchos ArrayView constructor.
 
    indices = ArrayView<const int>(eIndices, numEntries);
    values  = ArrayView<const double>(eValues, numEntries);
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getLocalDiagCopy");
    mtx_->ExtractDiagonalCopy(toEpetra<GlobalOrdinal, Node>(diag));
  }
 
  void getLocalDiagOffsets(Teuchos::ArrayRCP<size_t> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagOffsets() is not implemented or supported.");
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & /* diag */, const Teuchos::ArrayView<const size_t> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagCopy using offsets is not implemented or supported.");
  }
 
  void getLocalDiagCopy(Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & /* diag */, const Kokkos::View<const size_t *, typename Node::device_type, Kokkos::MemoryUnmanaged> & /* offsets */) const {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.getLocalDiagCopy using offsets is not implemented or supported.");
  }
 
  void replaceDiag(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &diag) {
    mtx_->ReplaceDiagonalValues(toEpetra<GlobalOrdinal, Node>(diag));
  }
 
  void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x) {
    XPETRA_ERR_CHECK(mtx_->LeftScale(toEpetra<GlobalOrdinal, Node>(x)));
  };
  void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &x) {
    XPETRA_ERR_CHECK(mtx_->RightScale(toEpetra<GlobalOrdinal, Node>(x)));
  };
 
 
 
 
  void apply(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X, MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &Y, Teuchos::ETransp mode = Teuchos::NO_TRANS, Scalar alpha = ScalarTraits<Scalar>::one(), Scalar beta = ScalarTraits<Scalar>::zero()) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::apply");
 
    // TEUCHOS_TEST_FOR_EXCEPTION((alpha != 1) || (beta != 0), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.multiply() only accept alpha==1 and beta==0");
 
    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, X, eX, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, Y, eY, "Xpetra::EpetraCrsMatrixT->apply() only accept Xpetra::EpetraMultiVectorT as input arguments.");
 
    TEUCHOS_TEST_FOR_EXCEPTION((mode != Teuchos::NO_TRANS) && (mode != Teuchos::TRANS), Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT->apply() only accept mode == NO_TRANS or mode == TRANS");
    bool eTrans = toEpetra(mode);
 
    // /!\ UseTranspose value
    TEUCHOS_TEST_FOR_EXCEPTION(mtx_->UseTranspose(), Xpetra::Exceptions::NotImplemented, "An exception is throw to let you know that Xpetra::EpetraCrsMatrixT->apply() do not take into account the UseTranspose() parameter of Epetra_CrsMatrix.");
 
    RCP<Epetra_MultiVector> epY = eY.getEpetra_MultiVector();
 
    // helper vector: tmp = A*x
    RCP<Epetra_MultiVector> tmp = Teuchos::rcp(new Epetra_MultiVector(*epY));
    tmp->PutScalar(0.0);
    XPETRA_ERR_CHECK(mtx_->Multiply(eTrans, *eX.getEpetra_MultiVector(), *tmp));
 
    // calculate alpha * A * x + beta * y
    XPETRA_ERR_CHECK(eY.getEpetra_MultiVector()->Update(alpha, *tmp, beta));
  }
 
  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 {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "Xpetra::EpetraCrsMatrixT.apply() is not implemented for region matricies.");
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getDomainMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->DomainMap());
  }
 
  const RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getRangeMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->RangeMap());
  }
 
 
 
 
  std::string description() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::description");
 
    // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
    std::ostringstream oss;
    // TODO: oss << DistObject<char, LocalOrdinal,GlobalOrdinal>::description();
    if (isFillComplete()) {
      oss << "{status = fill complete"
          << ", global rows = " << getGlobalNumRows()
          << ", global cols = " << getGlobalNumCols()
          << ", global num entries = " << getGlobalNumEntries()
          << "}";
    } else {
      oss << "{status = fill not complete"
          << ", global rows = " << getGlobalNumRows()
          << "}";
    }
    return oss.str();
  }
 
  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
    XPETRA_MONITOR("EpetraCrsMatrixT::describe");
 
    // This implementation come from Tpetra_CrsMatrix_def.hpp (without modification)
    using std::endl;
    using std::setw;
    using Teuchos::VERB_DEFAULT;
    using Teuchos::VERB_EXTREME;
    using Teuchos::VERB_HIGH;
    using Teuchos::VERB_LOW;
    using Teuchos::VERB_MEDIUM;
    using Teuchos::VERB_NONE;
    Teuchos::EVerbosityLevel vl = verbLevel;
    if (vl == VERB_DEFAULT) vl = VERB_LOW;
    RCP<const Comm<int> > comm = this->getComm();
    const int myImageID        = comm->getRank(),
              numImages        = comm->getSize();
    size_t width               = 1;
    for (size_t dec = 10; dec < getGlobalNumRows(); dec *= 10) {
      ++width;
    }
    width = std::max<size_t>(width, 11) + 2;
    Teuchos::OSTab tab(out);
    //    none: print nothing
    //     low: print O(1) info from node 0
    //  medium: print O(P) info, num entries per node
    //    high: print O(N) info, num entries per row
    // extreme: print O(NNZ) info: print indices and values
    //
    // for medium and higher, print constituent objects at specified verbLevel
    if (vl != VERB_NONE) {
      if (myImageID == 0) out << this->description() << std::endl;
      // O(1) globals, minus what was already printed by description()
      if (isFillComplete() && myImageID == 0) {
        out << "Global max number of entries = " << getGlobalMaxNumRowEntries() << std::endl;
      }
      // constituent objects
      if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
        if (myImageID == 0) out << "\nRow map: " << std::endl;
        getRowMap()->describe(out, vl);
        //
        if (getColMap() != null) {
          if (getColMap() == getRowMap()) {
            if (myImageID == 0) out << "\nColumn map is row map.";
          } else {
            if (myImageID == 0) out << "\nColumn map: " << std::endl;
            getColMap()->describe(out, vl);
          }
        }
        if (getDomainMap() != null) {
          if (getDomainMap() == getRowMap()) {
            if (myImageID == 0) out << "\nDomain map is row map.";
          } else if (getDomainMap() == getColMap()) {
            if (myImageID == 0) out << "\nDomain map is row map.";
          } else {
            if (myImageID == 0) out << "\nDomain map: " << std::endl;
            getDomainMap()->describe(out, vl);
          }
        }
        if (getRangeMap() != null) {
          if (getRangeMap() == getDomainMap()) {
            if (myImageID == 0) out << "\nRange map is domain map." << std::endl;
          } else if (getRangeMap() == getRowMap()) {
            if (myImageID == 0) out << "\nRange map is row map." << std::endl;
          } else {
            if (myImageID == 0) out << "\nRange map: " << std::endl;
            getRangeMap()->describe(out, vl);
          }
        }
        if (myImageID == 0) out << std::endl;
      }
      // O(P) data
      if (vl == VERB_MEDIUM || vl == VERB_HIGH || vl == VERB_EXTREME) {
        for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
          if (myImageID == imageCtr) {
            out << "Node ID = " << imageCtr << std::endl;
            // TODO: need a graph
            //               if (staticGraph_->indicesAreAllocated() == false) {
            //                 out << "Node not allocated" << std::endl;
            //               }
            //               else {
            //                 out << "Node number of allocated entries = " << staticGraph_->getLocalAllocationSize() << std::endl;
            //               }
 
            // TMP:
            //            const Epetra_CrsGraph & staticGraph_ = mtx_->Graph();
            // End of TMP
 
            out << "Node number of entries = " << getLocalNumEntries() << std::endl;
            out << "Node max number of entries = " << getLocalMaxNumRowEntries() << std::endl;
          }
          comm->barrier();
          comm->barrier();
          comm->barrier();
        }
      }
      // O(N) and O(NNZ) data
      if (vl == VERB_HIGH || vl == VERB_EXTREME) {
        for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
          if (myImageID == imageCtr) {
            out << std::setw(width) << "Node ID"
                << std::setw(width) << "Global Row"
                << std::setw(width) << "Num Entries";
            if (vl == VERB_EXTREME) {
              out << std::setw(width) << "(Index,Value)";
            }
            out << std::endl;
            for (size_t r = 0; r < getLocalNumRows(); ++r) {
              const size_t nE   = getNumEntriesInLocalRow(r);
              GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
              out << std::setw(width) << myImageID
                  << std::setw(width) << gid
                  << std::setw(width) << nE;
              if (vl == VERB_EXTREME) {
                if (isGloballyIndexed()) {
                  ArrayView<const GlobalOrdinal> rowinds;
                  ArrayView<const Scalar> rowvals;
                  getGlobalRowView(gid, rowinds, rowvals);
                  for (size_t j = 0; j < nE; ++j) {
                    out << " (" << rowinds[j]
                        << ", " << rowvals[j]
                        << ") ";
                  }
                } else if (isLocallyIndexed()) {
                  ArrayView<const LocalOrdinal> rowinds;
                  ArrayView<const Scalar> rowvals;
                  getLocalRowView(r, rowinds, rowvals);
                  for (size_t j = 0; j < nE; ++j) {
                    out << " (" << getColMap()->getGlobalElement(rowinds[j])
                        << ", " << rowvals[j]
                        << ") ";
                  }
                }
              }
              out << std::endl;
            }
          }
          comm->barrier();
          comm->barrier();
          comm->barrier();
        }
      }
    }
  }
 
  void setObjectLabel(const std::string &objectLabel) {
    Teuchos::LabeledObject::setObjectLabel(objectLabel);
    mtx_->SetLabel(objectLabel.c_str());
  }
 
  EpetraCrsMatrixT(const EpetraCrsMatrixT &matrix)
    : mtx_(Teuchos::rcp(new Epetra_CrsMatrix(*(matrix.mtx_))))
    , isFillResumed_(false) {}
 
  //{@
 
  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const {
    XPETRA_MONITOR("EpetraCrsMatrixT::getMap");
    return toXpetra<GlobalOrdinal, Node>(mtx_->Map());
  }
 
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
                const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
    int err                       = mtx_->Import(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
                const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
    int err                       = mtx_->Export(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doImport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &source,
                const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doImport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tSource.getEpetra_CrsMatrix();
    int err                       = mtx_->Import(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void doExport(const DistObject<char, LocalOrdinal, GlobalOrdinal, Node> &dest,
                const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraCrsMatrixT::doExport");
 
    XPETRA_DYNAMIC_CAST(const EpetraCrsMatrixT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraCrsMatrixT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraCrsMatrixT::doImport only accept Xpetra::EpetraImportT as input arguments.");
 
    RCP<const Epetra_CrsMatrix> v = tDest.getEpetra_CrsMatrix();
    int err                       = mtx_->Export(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }
 
  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &newMap) {
    XPETRA_MONITOR("EpetraCrsMatrixT::removeEmptyProcessesInPlace");
    const Epetra_Map *newMapEpetra = (!newMap.is_null()) ? &toEpetra<GlobalOrdinal, Node>(newMap) : 0;
    int err                        = mtx_->RemoveEmptyProcessesInPlace(newMapEpetra);
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Caught nonzero error code " + std::to_string(err) + " returned by Epetra.");
  }
 
 
 
 
  bool hasMatrix() const { return !mtx_.is_null(); }
 
  EpetraCrsMatrixT(const Teuchos::RCP<Epetra_CrsMatrix> &mtx)
    : mtx_(mtx)
    , isFillResumed_(false) {}
 
  RCP<const Epetra_CrsMatrix> getEpetra_CrsMatrix() const { return mtx_; }
 
  RCP<Epetra_CrsMatrix> getEpetra_CrsMatrixNonConst() const { return mtx_; }  // TODO: remove
 
#ifdef HAVE_XPETRA_TPETRA
  local_matrix_type getLocalMatrix() const {
    RCP<Epetra_CrsMatrix> matrix = getEpetra_CrsMatrixNonConst();
 
    const int numRows = matrix->NumMyRows();
    const int numCols = matrix->NumMyCols();
    const int nnz     = matrix->NumMyNonzeros();
 
    int *rowptr;
    int *colind;
    double *vals;
    int rv = matrix->ExtractCrsDataPointers(rowptr, colind, vals);
    TEUCHOS_TEST_FOR_EXCEPTION(rv, std::runtime_error, "Xpetra::CrsMatrix<>::getLocalMatrix: failed in ExtractCrsDataPointers");
 
    // Transform int* rowptr array to size_type* array
    typename local_matrix_type::row_map_type::non_const_type kokkosRowPtr(Kokkos::ViewAllocateWithoutInitializing("local row map"), numRows + 1);
    for (size_t i = 0; i < kokkosRowPtr.size(); i++)
      kokkosRowPtr(i) = Teuchos::asSafe<typename local_matrix_type::row_map_type::value_type>(rowptr[i]);
 
    // create Kokkos::Views
    typename local_matrix_type::index_type kokkosColind(colind, nnz);
    typename local_matrix_type::values_type kokkosVals(vals, nnz);
 
    local_matrix_type localMatrix = local_matrix_type("LocalMatrix", numRows, numCols, nnz, kokkosVals, kokkosRowPtr, kokkosColind);
 
    return localMatrix;
  }
 
  void setAllValues(const typename local_matrix_type::row_map_type &ptr,
                    const typename local_matrix_type::StaticCrsGraphType::entries_type::non_const_type &ind,
                    const typename local_matrix_type::values_type &val) {
    // Check sizes
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(ptr.size()) != getLocalNumRows() + 1, Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that the size of your rowptr array is incorrect.");
    TEUCHOS_TEST_FOR_EXCEPTION(val.size() != ind.size(), Xpetra::Exceptions::RuntimeError,
                               "An exception is thrown to let you know that you mismatched your pointers.");
 
    // Check pointers
    if (val.size() > 0) {
      TEUCHOS_TEST_FOR_EXCEPTION(ind.data() != mtx_->ExpertExtractIndices().Values(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
      TEUCHOS_TEST_FOR_EXCEPTION(val.data() != mtx_->ExpertExtractValues(), Xpetra::Exceptions::RuntimeError,
                                 "An exception is thrown to let you know that you mismatched your pointers.");
    }
 
    // We have to make a copy here, it is unavoidable
    // See comments in allocateAllValues
    const size_t N = getLocalNumRows();
 
    Epetra_IntSerialDenseVector &myRowptr = mtx_->ExpertExtractIndexOffset();
    myRowptr.Resize(N + 1);
    for (size_t i = 0; i < N + 1; i++)
      myRowptr[i] = Teuchos::as<int>(ptr(i));
  }
 
  LocalOrdinal GetStorageBlockSize() const { return 1; }
 
 private:
#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
 
 
  void residual(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &X,
                const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B,
                MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &R) const {
    Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one;
    apply(X, R);
    R.update(one, B, negone);
  }
 
 private:
  RCP<Epetra_CrsMatrix> mtx_;
 
  bool isFillResumed_;  //< For Epetra, fillResume() is a fictive operation but we need to keep track of it. This boolean is true only is resumeFill() have been called and fillComplete() have not been called afterward.
 
};  // EpetraCrsMatrixT class
 
#endif  // #ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
 
}  // namespace Xpetra
 
#define XPETRA_EPETRACRSMATRIX_SHORT
#endif  // XPETRA_EPETRACRSMATRIX_HPP