Xpetra_EpetraMultiVector.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//             Xpetra: A linear algebra interface package
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef XPETRA_EPETRAMULTIVECTOR_HPP
#define XPETRA_EPETRAMULTIVECTOR_HPP

/* this file is automatically generated - do not edit (see script/epetra.py) */

#include <Kokkos_Core.hpp>
#include <Kokkos_DualView.hpp>

#include "Xpetra_EpetraConfigDefs.hpp"

#include "Xpetra_MultiVector.hpp"
#include "Xpetra_Vector.hpp"

#include "Xpetra_EpetraMap.hpp"
#include "Xpetra_EpetraExport.hpp"
#include "Xpetra_Utils.hpp"
#include "Xpetra_EpetraUtils.hpp"
#include "Xpetra_EpetraImport.hpp"
#include "Xpetra_Exceptions.hpp"
#include "Epetra_SerialComm.h"

#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>

#if defined(XPETRA_ENABLE_DEPRECATED_CODE)
#ifdef __GNUC__
#warning "The header file Trilinos/packages/xpetra/src/MultiVector/Xpetra_EpetraMultiVector.hpp is deprecated."
#endif
#else
#error "The header file Trilinos/packages/xpetra/src/MultiVector/Xpetra_EpetraMultiVector.hpp is deprecated."
#endif

namespace Xpetra {

// TODO: move that elsewhere
template <class GlobalOrdinal, class Node>
XPETRA_DEPRECATED const Epetra_MultiVector &toEpetra(const MultiVector<double, int, GlobalOrdinal, Node> &);
template <class GlobalOrdinal, class Node>
XPETRA_DEPRECATED Epetra_MultiVector &toEpetra(MultiVector<double, int, GlobalOrdinal, Node> &);
template <class GlobalOrdinal, class Node>
XPETRA_DEPRECATED RCP<MultiVector<double, int, GlobalOrdinal, Node> > toXpetra(RCP<Epetra_MultiVector> vec);

// we need this forward declaration
#ifndef DOXYGEN_SHOULD_SKIP_THIS
template <class GlobalOrdinal, class Node>
class EpetraVectorT;
#endif

template <class EpetraGlobalOrdinal, class Node>
class XPETRA_DEPRECATED EpetraMultiVectorT
  : public virtual MultiVector<double, int, EpetraGlobalOrdinal, Node> {
  typedef double Scalar;
  typedef int LocalOrdinal;
  typedef EpetraGlobalOrdinal GlobalOrdinal;

 public:


  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, size_t NumVectors, bool zeroOut = true) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraMultiVector only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }

  EpetraMultiVectorT(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Teuchos::DataAccess copyOrView = Teuchos::Copy) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraMultiVector only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }

  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> > &ArrayOfPtrs, size_t NumVectors) {
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraMultiVector only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }

  virtual ~EpetraMultiVectorT() {}




  void replaceGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {}

  void sumIntoGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {}

  void replaceLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {}

  void sumIntoLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {}

  void putScalar(const Scalar &value) {}




  Teuchos::RCP<const Vector<double, int, GlobalOrdinal, Node> > getVector(size_t j) const {
    return Teuchos::null;
  }

  Teuchos::RCP<Vector<double, int, GlobalOrdinal, Node> > getVectorNonConst(size_t j) {
    return Teuchos::null;
  }

  Teuchos::ArrayRCP<const Scalar> getData(size_t j) const {
    return ArrayRCP<const Scalar>();
  }

  Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j) {
    return ArrayRCP<Scalar>();
  }




  void dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Teuchos::ArrayView<Scalar> &dots) const {}

  void abs(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {}

  void reciprocal(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {}

  void scale(const Scalar &alpha) {}

  void scale(Teuchos::ArrayView<const Scalar> alpha) {}

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta) {}

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &gamma) {}

  void norm1(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {}

  void norm2(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {}

  void normInf(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {}

  void meanValue(const Teuchos::ArrayView<Scalar> &means) const {}

  void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &beta) {}

  void elementWiseMultiply(Scalar scalarAB, const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, Scalar scalarThis) {}




  size_t getNumVectors() const { return 0; }

  size_t getLocalLength() const { return 0; }

  global_size_t getGlobalLength() const { return 0; }

  // \brief Checks to see if the local length, number of vectors and size of Scalar type match
  bool isSameSize(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &vec) const { return false; }




  std::string description() const { return std::string(""); }

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


  void randomize(bool bUseXpetraImplementation = false) {}

  void randomize(const Scalar &minVal, const Scalar &maxVal, bool bUseXpetraImplementation = false) {}

  //{@

  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const { return Teuchos::null; }

  void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {}

  void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {}

  void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {}

  void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {}

  void replaceMap(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map) {}




  EpetraMultiVectorT(const RCP<Epetra_MultiVector> &vec) {  // TODO removed const
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::RuntimeError,
                               "Xpetra::EpetraMultiVector only available for GO=int or GO=long long with EpetraNode (Serial or OpenMP depending on configuration)");
  }

  RCP<Epetra_MultiVector> getEpetra_MultiVector() const { return Teuchos::null; }

  void setSeed(unsigned int seed) {}


 protected:
  virtual void
  assign(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &rhs) {}

};  // EpetraMultiVectorT class

// specialization on GO=int and Node=EpetraNode
#ifndef XPETRA_EPETRA_NO_32BIT_GLOBAL_INDICES
template <>
class EpetraMultiVectorT<int, EpetraNode>
  : public virtual MultiVector<double, int, int, EpetraNode> {
  typedef double Scalar;
  typedef int LocalOrdinal;
  typedef int GlobalOrdinal;
  typedef EpetraNode Node;

 public:


  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, size_t NumVectors, bool zeroOut = true)
    : vec_(Teuchos::rcp(new Epetra_MultiVector(toEpetra<GlobalOrdinal, Node>(map), Teuchos::as<int>(NumVectors), zeroOut))) {}

  EpetraMultiVectorT(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Teuchos::DataAccess copyOrView = Teuchos::Copy) {
    if (copyOrView == Teuchos::Copy)
      vec_ = Teuchos::rcp(new Epetra_MultiVector(toEpetra<GlobalOrdinal, Node>(source)));
    else {
      int *indices = new int[source.getNumVectors()];
      for (size_t i = 0; i < source.getNumVectors(); i++)
        indices[i] = i;
      vec_ = Teuchos::rcp(new Epetra_MultiVector(View, toEpetra<GlobalOrdinal, Node>(source), indices, source.getNumVectors()));
      delete[] indices;
    }
  }

  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> > &ArrayOfPtrs, size_t NumVectors) {
    // TODO: input argument 'NumVectors' is not necessary in both Xpetra and Tpetra interface. Should it be removed?

    const std::string tfecfFuncName("MultiVector(ArrayOfPtrs)");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(NumVectors < 1 || NumVectors != Teuchos::as<size_t>(ArrayOfPtrs.size()), std::runtime_error,
                                          ": ArrayOfPtrs.size() must be strictly positive and as large as ArrayOfPtrs.");

#ifdef HAVE_XPETRA_DEBUG
    // This cannot be tested by Epetra itself
    {
      size_t localLength = map->getLocalNumElements();
      for (int j = 0; j < ArrayOfPtrs.size(); j++) {
        TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != localLength, std::runtime_error,
                                              ": ArrayOfPtrs[" << j << "].size() (== " << ArrayOfPtrs[j].size() << ") is not equal to getLocalLength() (== " << localLength);
      }
    }
#endif

    // Convert Teuchos::ArrayView< const Teuchos::ArrayView< const Scalar > > to double**
    Array<const double *> arrayOfRawPtrs(ArrayOfPtrs.size());
    for (int i = 0; i < ArrayOfPtrs.size(); i++) {
      arrayOfRawPtrs[i] = ArrayOfPtrs[i].getRawPtr();
    }
    double **rawArrayOfRawPtrs = const_cast<double **>(arrayOfRawPtrs.getRawPtr());  // This const_cast should be fine, because Epetra_DataAccess=Copy.

    vec_ = Teuchos::rcp(new Epetra_MultiVector(Copy, toEpetra<GlobalOrdinal, Node>(map), rawArrayOfRawPtrs, static_cast<int>(NumVectors)));
  }

  virtual ~EpetraMultiVectorT() {}




  void replaceGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceGlobalValue");
    vec_->ReplaceGlobalValue(globalRow, Teuchos::as<int>(vectorIndex), value);
  }

  void sumIntoGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::sumIntoGlobalValue");
    vec_->SumIntoGlobalValue(globalRow, Teuchos::as<int>(vectorIndex), value);
  }

  void replaceLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceLocalValue");
    vec_->ReplaceMyValue(myRow, Teuchos::as<int>(vectorIndex), value);
  }

  void sumIntoLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::sumIntoLocalValue");
    vec_->SumIntoMyValue(myRow, Teuchos::as<int>(vectorIndex), value);
  }

  void putScalar(const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::putScalar");
    vec_->PutScalar(value);
  }




  Teuchos::RCP<const Vector<double, int, int, EpetraNode> > getVector(size_t j) const;

  Teuchos::RCP<Vector<double, int, int, EpetraNode> > getVectorNonConst(size_t j);

  Teuchos::ArrayRCP<const Scalar> getData(size_t j) const {
    XPETRA_MONITOR("EpetraMultiVectorT::getData");

    double **arrayOfPointers;

    vec_->ExtractView(&arrayOfPointers);

    double *data    = arrayOfPointers[j];
    int localLength = vec_->MyLength();

    return ArrayRCP<double>(data, 0, localLength, false);  // no ownership
  }

  Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j) {
    XPETRA_MONITOR("EpetraMultiVectorT::getDataNonConst");

    double **arrayOfPointers;

    vec_->ExtractView(&arrayOfPointers);

    double *data    = arrayOfPointers[j];
    int localLength = vec_->MyLength();

    return ArrayRCP<double>(data, 0, localLength, false);  // no ownership
  }




  void dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Teuchos::ArrayView<Scalar> &dots) const {
    XPETRA_MONITOR("EpetraMultiVectorT::dot");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, A, eA, "This Xpetra::EpetraMultiVectorT method only accept Xpetra::EpetraMultiVectorT as input arguments.");
    vec_->Dot(*eA.getEpetra_MultiVector(), dots.getRawPtr());
  }

  void abs(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
    XPETRA_MONITOR("EpetraMultiVectorT::abs");
    vec_->Abs(toEpetra<GlobalOrdinal, Node>(A));
  }

  void reciprocal(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
    XPETRA_MONITOR("EpetraMultiVectorT::reciprocal");
    vec_->Reciprocal(toEpetra<GlobalOrdinal, Node>(A));
  }

  void scale(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraMultiVectorT::scale");
    vec_->Scale(alpha);
  }

  void scale(Teuchos::ArrayView<const Scalar> alpha) {
    XPETRA_MONITOR("EpetraMultiVectorT::scale");
    // Epetra, unlike Tpetra, doesn't implement this version of
    // scale().  Deal with this by scaling one column at a time.
    const size_t numVecs = this->getNumVectors();
    for (size_t j = 0; j < numVecs; ++j) {
      Epetra_Vector *v = (*vec_)(j);
      v->Scale(alpha[j]);
    }
  }

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta) {
    XPETRA_MONITOR("EpetraMultiVectorT::update");
    vec_->Update(alpha, toEpetra<GlobalOrdinal, Node>(A), beta);
  }

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &gamma) {
    XPETRA_MONITOR("EpetraMultiVectorT::update");
    vec_->Update(alpha, toEpetra<GlobalOrdinal, Node>(A), beta, toEpetra<GlobalOrdinal, Node>(B), gamma);
  }

  void norm1(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::norm1");
    vec_->Norm1(norms.getRawPtr());
  }

  void norm2(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::norm2");
    vec_->Norm2(norms.getRawPtr());
  }

  void normInf(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::normInf");
    vec_->NormInf(norms.getRawPtr());
  }

  void meanValue(const Teuchos::ArrayView<Scalar> &means) const {
    XPETRA_MONITOR("EpetraMultiVectorT::meanValue");
    vec_->MeanValue(means.getRawPtr());
  }  // TODO: modify ArrayView size ??

  void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &beta) {
    XPETRA_MONITOR("EpetraMultiVectorT::multiply");
    vec_->Multiply(toEpetra(transA), toEpetra(transB), alpha, toEpetra(A), toEpetra(B), beta);
  }

  void elementWiseMultiply(Scalar scalarAB, const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, Scalar scalarThis) {
    XPETRA_MONITOR("EpetraMultiVectorT::elementWiseMultiply");
    vec_->Multiply(scalarAB, toEpetra<GlobalOrdinal, Node>(A), toEpetra<GlobalOrdinal, Node>(B), scalarThis);
  }




  size_t getNumVectors() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getNumVectors");
    return vec_->NumVectors();
  }

  size_t getLocalLength() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getLocalLength");
    return vec_->MyLength();
  }

  global_size_t getGlobalLength() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getGlobalLength");
    return vec_->GlobalLength64();
  }

  bool isSameSize(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &vec) const {
    XPETRA_MONITOR("EpetraMultiVectorT::isSameSize");
    auto vv = toEpetra<GlobalOrdinal, Node>(vec);
    return ((vec_->MyLength() == vv.MyLength()) && (vec_->NumVectors() == vv.NumVectors()));
  }




  std::string description() const {
    XPETRA_MONITOR("EpetraMultiVectorT::description");
    TEUCHOS_TEST_FOR_EXCEPTION(1, Xpetra::Exceptions::NotImplemented, "TODO");
    TEUCHOS_UNREACHABLE_RETURN("TODO");
  }

  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel /* verbLevel */ = Teuchos::Describable::verbLevel_default) const {
    XPETRA_MONITOR("EpetraMultiVectorT::describe");
    vec_->Print(out);
  }


  void randomize(bool bUseXpetraImplementation = false) {
    XPETRA_MONITOR("EpetraMultiVectorT::randomize");

    if (bUseXpetraImplementation)
      Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Xpetra_randomize();
    else
      vec_->Random();
  }

  void randomize(const Scalar &minVal, const Scalar &maxVal, bool bUseXpetraImplementation = false) {
    XPETRA_MONITOR("EpetraMultiVectorT::randomize");

    if (bUseXpetraImplementation)
      Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Xpetra_randomize(minVal, maxVal);
    else {
      vec_->Random();
      const size_t numVectors = getNumVectors();
      for (size_t i = 0; i < numVectors; i++) {
        Teuchos::ArrayRCP<Scalar> datai = getDataNonConst(i);

        const size_t myLength = getLocalLength();
        for (size_t j = 0; j < myLength; j++) {
          datai[j] = 0.5 * (maxVal - minVal) * datai[j] + 0.5 * (maxVal + minVal);
        }
      }
    }
  }

  //{@

  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getMap");
    return toXpetra<GlobalOrdinal, Node>(vec_->Map());
  }

  void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doImport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tSource.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->Import(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra is " << err);
  }

  void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doExport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tDest.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->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<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doImport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tSource.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->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<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doExport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tDest.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->Export(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }

  void replaceMap(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceMap");
    int err = 0;
    if (!map.is_null()) {
      err = this->getEpetra_MultiVector()->ReplaceMap(toEpetra<GlobalOrdinal, Node>(map));

    } else {
      // Replace map with a dummy map to avoid potential hangs later
      Epetra_SerialComm SComm;
      Epetra_Map NewMap((GlobalOrdinal)vec_->MyLength(), (GlobalOrdinal)vec_->Map().IndexBase64(), SComm);
      err = this->getEpetra_MultiVector()->ReplaceMap(NewMap);
    }
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }




  EpetraMultiVectorT(const RCP<Epetra_MultiVector> &vec)
    : vec_(vec) {}  // TODO removed const

  RCP<Epetra_MultiVector> getEpetra_MultiVector() const { return vec_; }

  void setSeed(unsigned int seed) {
    XPETRA_MONITOR("EpetraMultiVectorT::seedrandom");

    Teuchos::ScalarTraits<Scalar>::seedrandom(seed);
    vec_->SetSeed(seed);
  }

  typedef typename Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::dual_view_type dual_view_type;

  typename dual_view_type::t_host_const_um getHostLocalView(Access::ReadOnlyStruct) const override { return getHostLocalView(Access::ReadWrite); }

  typename dual_view_type::t_dev_const_um getDeviceLocalView(Access::ReadOnlyStruct) const override { return getDeviceLocalView(Access::ReadWrite); }

  typename dual_view_type::t_host_um getHostLocalView(Access::OverwriteAllStruct) const override { return getHostLocalView(Access::ReadWrite); }

  typename dual_view_type::t_dev_um getDeviceLocalView(Access::OverwriteAllStruct) const override { return getDeviceLocalView(Access::ReadWrite); }

  typename dual_view_type::t_host_um getHostLocalView(Access::ReadWriteStruct) const override {
    typedef Kokkos::View<typename dual_view_type::t_host::data_type,
                         Kokkos::LayoutLeft,
                         typename dual_view_type::t_host::device_type,
                         Kokkos::MemoryUnmanaged>
        epetra_view_type;

    // access Epetra multivector data
    double *data = NULL;
    int myLDA;
    vec_->ExtractView(&data, &myLDA);
    int localLength = vec_->MyLength();
    int numVectors  = getNumVectors();

    // create view
    epetra_view_type test                  = epetra_view_type(data, localLength, numVectors);
    typename dual_view_type::t_host_um ret = subview(test, Kokkos::ALL(), Kokkos::ALL());

    return ret;
  }

  typename dual_view_type::t_dev_um getDeviceLocalView(Access::ReadWriteStruct) const override { return getHostLocalView(Access::ReadWrite); }


 protected:
  virtual void
  assign(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &rhs) {
    typedef EpetraMultiVectorT this_type;
    const this_type *rhsPtr = dynamic_cast<const this_type *>(&rhs);
    TEUCHOS_TEST_FOR_EXCEPTION(
        rhsPtr == NULL, std::invalid_argument,
        "Xpetra::MultiVector::operator=: "
        "The left-hand side (LHS) of the assignment has a different type than "
        "the right-hand side (RHS).  The LHS has type Xpetra::EpetraMultiVectorT "
        "(which means it wraps an Epetra_MultiVector), but the RHS has some "
        "other type.  This probably means that the RHS wraps a Tpetra::Multi"
        "Vector.  Xpetra::MultiVector does not currently implement assignment "
        "from a Tpetra object to an Epetra object, though this could be added "
        "with sufficient interest.");

    RCP<const Epetra_MultiVector> rhsImpl = rhsPtr->getEpetra_MultiVector();
    RCP<Epetra_MultiVector> lhsImpl       = this->getEpetra_MultiVector();

    TEUCHOS_TEST_FOR_EXCEPTION(
        rhsImpl.is_null(), std::logic_error,
        "Xpetra::MultiVector::operator= "
        "(in Xpetra::EpetraMultiVectorT::assign): *this (the right-hand side of "
        "the assignment) has a null RCP<Epetra_MultiVector> inside.  Please "
        "report this bug to the Xpetra developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        lhsImpl.is_null(), std::logic_error,
        "Xpetra::MultiVector::operator= "
        "(in Xpetra::EpetraMultiVectorT::assign): The left-hand side of the "
        "assignment has a null RCP<Epetra_MultiVector> inside.  Please report "
        "this bug to the Xpetra developers.");

    // Epetra_MultiVector's assignment operator does a deep copy.
    *lhsImpl = *rhsImpl;
  }

 private:
  RCP<Epetra_MultiVector> vec_;

};  // EpetraMultiVectorT class (specialization on GO=int, NO=EpetraNode
#endif

// specialization on GO=long long and EpetraNode
#ifndef XPETRA_EPETRA_NO_64BIT_GLOBAL_INDICES
template <>
class EpetraMultiVectorT<long long, EpetraNode>
  : public virtual MultiVector<double, int, long long, EpetraNode> {
  typedef double Scalar;
  typedef int LocalOrdinal;
  typedef long long GlobalOrdinal;
  typedef EpetraNode Node;

 public:


  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, size_t NumVectors, bool zeroOut = true)
    : vec_(Teuchos::rcp(new Epetra_MultiVector(toEpetra<GlobalOrdinal, Node>(map), Teuchos::as<int>(NumVectors), zeroOut))) {}

  EpetraMultiVectorT(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source)
    : vec_(Teuchos::rcp(new Epetra_MultiVector(toEpetra<GlobalOrdinal, Node>(source)))) {}

  EpetraMultiVectorT(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map, const Teuchos::ArrayView<const Teuchos::ArrayView<const Scalar> > &ArrayOfPtrs, size_t NumVectors) {
    // TODO: input argument 'NumVectors' is not necessary in both Xpetra and Tpetra interface. Should it be removed?

    const std::string tfecfFuncName("MultiVector(ArrayOfPtrs)");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(NumVectors < 1 || NumVectors != Teuchos::as<size_t>(ArrayOfPtrs.size()), std::runtime_error,
                                          ": ArrayOfPtrs.size() must be strictly positive and as large as ArrayOfPtrs.");

#ifdef HAVE_XPETRA_DEBUG
    // This cannot be tested by Epetra itself
    {
      size_t localLength = map->getLocalNumElements();
      for (int j = 0; j < ArrayOfPtrs.size(); j++) {
        TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(Teuchos::as<size_t>(ArrayOfPtrs[j].size()) != localLength, std::runtime_error,
                                              ": ArrayOfPtrs[" << j << "].size() (== " << ArrayOfPtrs[j].size() << ") is not equal to getLocalLength() (== " << localLength);
      }
    }
#endif

    // Convert Teuchos::ArrayView< const Teuchos::ArrayView< const Scalar > > to double**
    Array<const double *> arrayOfRawPtrs(ArrayOfPtrs.size());
    for (int i = 0; i < ArrayOfPtrs.size(); i++) {
      arrayOfRawPtrs[i] = ArrayOfPtrs[i].getRawPtr();
    }
    double **rawArrayOfRawPtrs = const_cast<double **>(arrayOfRawPtrs.getRawPtr());  // This const_cast should be fine, because Epetra_DataAccess=Copy.

    vec_ = Teuchos::rcp(new Epetra_MultiVector(Copy, toEpetra<GlobalOrdinal, Node>(map), rawArrayOfRawPtrs, NumVectors));
  }

  virtual ~EpetraMultiVectorT() {}




  void replaceGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceGlobalValue");
    vec_->ReplaceGlobalValue(globalRow, Teuchos::as<int>(vectorIndex), value);
  }

  void sumIntoGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::sumIntoGlobalValue");
    vec_->SumIntoGlobalValue(globalRow, Teuchos::as<int>(vectorIndex), value);
  }

  void replaceLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceLocalValue");
    vec_->ReplaceMyValue(myRow, Teuchos::as<int>(vectorIndex), value);
  }

  void sumIntoLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::sumIntoLocalValue");
    vec_->SumIntoMyValue(myRow, Teuchos::as<int>(vectorIndex), value);
  }

  void putScalar(const Scalar &value) {
    XPETRA_MONITOR("EpetraMultiVectorT::putScalar");
    vec_->PutScalar(value);
  }




  Teuchos::RCP<const Vector<double, int, long long, EpetraNode> > getVector(size_t j) const;

  Teuchos::RCP<Vector<double, int, long long, EpetraNode> > getVectorNonConst(size_t j);

  Teuchos::ArrayRCP<const Scalar> getData(size_t j) const {
    XPETRA_MONITOR("EpetraMultiVectorT::getData");

    double **arrayOfPointers;

    vec_->ExtractView(&arrayOfPointers);

    double *data    = arrayOfPointers[j];
    int localLength = vec_->MyLength();

    return ArrayRCP<double>(data, 0, localLength, false);  // no ownership
  }

  Teuchos::ArrayRCP<Scalar> getDataNonConst(size_t j) {
    XPETRA_MONITOR("EpetraMultiVectorT::getDataNonConst");

    double **arrayOfPointers;

    vec_->ExtractView(&arrayOfPointers);

    double *data    = arrayOfPointers[j];
    int localLength = vec_->MyLength();

    return ArrayRCP<double>(data, 0, localLength, false);  // no ownership
  }




  void dot(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Teuchos::ArrayView<Scalar> &dots) const {
    XPETRA_MONITOR("EpetraMultiVectorT::dot");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, A, eA, "This Xpetra::EpetraMultiVectorT method only accept Xpetra::EpetraMultiVectorT as input arguments.");
    vec_->Dot(*eA.getEpetra_MultiVector(), dots.getRawPtr());
  }

  void abs(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
    XPETRA_MONITOR("EpetraMultiVectorT::abs");
    vec_->Abs(toEpetra<GlobalOrdinal, Node>(A));
  }

  void reciprocal(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A) {
    XPETRA_MONITOR("EpetraMultiVectorT::reciprocal");
    vec_->Reciprocal(toEpetra<GlobalOrdinal, Node>(A));
  }

  void scale(const Scalar &alpha) {
    XPETRA_MONITOR("EpetraMultiVectorT::scale");
    vec_->Scale(alpha);
  }

  void scale(Teuchos::ArrayView<const Scalar> alpha) {
    XPETRA_MONITOR("EpetraMultiVectorT::scale");
    // Epetra, unlike Tpetra, doesn't implement this version of
    // scale().  Deal with this by scaling one column at a time.
    const size_t numVecs = this->getNumVectors();
    for (size_t j = 0; j < numVecs; ++j) {
      Epetra_Vector *v = (*vec_)(j);
      v->Scale(alpha[j]);
    }
  }

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta) {
    XPETRA_MONITOR("EpetraMultiVectorT::update");
    vec_->Update(alpha, toEpetra<GlobalOrdinal, Node>(A), beta);
  }

  void update(const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const Scalar &beta, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &gamma) {
    XPETRA_MONITOR("EpetraMultiVectorT::update");
    vec_->Update(alpha, toEpetra<GlobalOrdinal, Node>(A), beta, toEpetra<GlobalOrdinal, Node>(B), gamma);
  }

  void norm1(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::norm1");
    vec_->Norm1(norms.getRawPtr());
  }

  void norm2(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::norm2");
    vec_->Norm2(norms.getRawPtr());
  }

  void normInf(const Teuchos::ArrayView<Teuchos::ScalarTraits<Scalar>::magnitudeType> &norms) const {
    XPETRA_MONITOR("EpetraMultiVectorT::normInf");
    vec_->NormInf(norms.getRawPtr());
  }

  void meanValue(const Teuchos::ArrayView<Scalar> &means) const {
    XPETRA_MONITOR("EpetraMultiVectorT::meanValue");
    vec_->MeanValue(means.getRawPtr());
  }  // TODO: modify ArrayView size ??

  void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const Scalar &alpha, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, const Scalar &beta) {
    XPETRA_MONITOR("EpetraMultiVectorT::multiply");
    vec_->Multiply(toEpetra(transA), toEpetra(transB), alpha, toEpetra(A), toEpetra(B), beta);
  }

  void elementWiseMultiply(Scalar scalarAB, const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &A, const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &B, Scalar scalarThis) {
    XPETRA_MONITOR("EpetraMultiVectorT::elementWiseMultiply");
    vec_->Multiply(scalarAB, toEpetra<GlobalOrdinal, Node>(A), toEpetra<GlobalOrdinal, Node>(B), scalarThis);
  }




  size_t getNumVectors() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getNumVectors");
    return vec_->NumVectors();
  }

  size_t getLocalLength() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getLocalLength");
    return vec_->MyLength();
  }

  global_size_t getGlobalLength() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getGlobalLength");
    return vec_->GlobalLength64();
  }

  bool isSameSize(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &vec) const {
    XPETRA_MONITOR("EpetraMultiVectorT::isSameSize");
    auto vv = toEpetra<GlobalOrdinal, Node>(vec);
    return ((vec_->MyLength() == vv.MyLength()) && (vec_->NumVectors() == vv.NumVectors()));
  }




  std::string description() const {
    XPETRA_MONITOR("EpetraMultiVectorT::description");
    TEUCHOS_TEST_FOR_EXCEPTION(true, Xpetra::Exceptions::NotImplemented, "TODO");
    // return "TODO"; // unreachable
  }

  void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel /* verbLevel */ = Teuchos::Describable::verbLevel_default) const {
    XPETRA_MONITOR("EpetraMultiVectorT::describe");
    vec_->Print(out);
  }


  void randomize(bool bUseXpetraImplementation = false) {
    XPETRA_MONITOR("EpetraMultiVectorT::randomize");

    if (bUseXpetraImplementation)
      Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Xpetra_randomize();
    else
      vec_->Random();
  }

  void randomize(const Scalar &minVal, const Scalar &maxVal, bool bUseXpetraImplementation = false) {
    XPETRA_MONITOR("EpetraMultiVectorT::randomize");

    if (bUseXpetraImplementation)
      Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Xpetra_randomize(minVal, maxVal);
    else {
      vec_->Random();
      const size_t numVectors = getNumVectors();
      for (size_t i = 0; i < numVectors; i++) {
        Teuchos::ArrayRCP<Scalar> datai = getDataNonConst(i);

        const size_t myLength = getLocalLength();
        for (size_t j = 0; j < myLength; j++) {
          datai[j] = 0.5 * (maxVal - minVal) * datai[j] + 0.5 * (maxVal + minVal);
        }
      }
    }
  }

  //{@

  Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > getMap() const {
    XPETRA_MONITOR("EpetraMultiVectorT::getMap");
    return toXpetra<GlobalOrdinal, Node>(vec_->Map());
  }

  void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doImport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tSource.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->Import(*v, *tImporter.getEpetra_Import(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra is " << err);
  }

  void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Import<LocalOrdinal, GlobalOrdinal, Node> &importer, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doExport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraImportT<GlobalOrdinal XPETRA_COMMA Node>, importer, tImporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tDest.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->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<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doImport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, source, tSource, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tSource.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->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<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Export<LocalOrdinal, GlobalOrdinal, Node> &exporter, CombineMode CM) {
    XPETRA_MONITOR("EpetraMultiVectorT::doExport");

    XPETRA_DYNAMIC_CAST(const EpetraMultiVectorT<GlobalOrdinal XPETRA_COMMA Node>, dest, tDest, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraMultiVectorT as input arguments.");
    XPETRA_DYNAMIC_CAST(const EpetraExportT<GlobalOrdinal XPETRA_COMMA Node>, exporter, tExporter, "Xpetra::EpetraMultiVectorT::doImport only accept Xpetra::EpetraImportT as input arguments.");

    RCP<Epetra_MultiVector> v = tDest.getEpetra_MultiVector();
    int err                   = this->getEpetra_MultiVector()->Export(*v, *tExporter.getEpetra_Export(), toEpetra(CM));
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }

  void replaceMap(const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> > &map) {
    XPETRA_MONITOR("EpetraMultiVectorT::replaceMap");
    int err = 0;
    if (!map.is_null()) {
      err = this->getEpetra_MultiVector()->ReplaceMap(toEpetra<GlobalOrdinal, Node>(map));

    } else {
      // Replace map with a dummy map to avoid potential hangs later
      Epetra_SerialComm SComm;
      Epetra_Map NewMap((GlobalOrdinal)vec_->MyLength(), (GlobalOrdinal)vec_->Map().IndexBase64(), SComm);
      err = this->getEpetra_MultiVector()->ReplaceMap(NewMap);
    }
    TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, "Catch error code returned by Epetra.");
  }




  EpetraMultiVectorT(const RCP<Epetra_MultiVector> &vec)
    : vec_(vec) {}  // TODO removed const

  RCP<Epetra_MultiVector> getEpetra_MultiVector() const { return vec_; }

  void setSeed(unsigned int seed) {
    XPETRA_MONITOR("EpetraMultiVectorT::seedrandom");

    Teuchos::ScalarTraits<Scalar>::seedrandom(seed);
    vec_->SetSeed(seed);
  }


 protected:
  virtual void
  assign(const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &rhs) {
    typedef EpetraMultiVectorT this_type;
    const this_type *rhsPtr = dynamic_cast<const this_type *>(&rhs);
    TEUCHOS_TEST_FOR_EXCEPTION(
        rhsPtr == NULL, std::invalid_argument,
        "Xpetra::MultiVector::operator=: "
        "The left-hand side (LHS) of the assignment has a different type than "
        "the right-hand side (RHS).  The LHS has type Xpetra::EpetraMultiVectorT "
        "(which means it wraps an Epetra_MultiVector), but the RHS has some "
        "other type.  This probably means that the RHS wraps a Tpetra::Multi"
        "Vector.  Xpetra::MultiVector does not currently implement assignment "
        "from a Tpetra object to an Epetra object, though this could be added "
        "with sufficient interest.");

    RCP<const Epetra_MultiVector> rhsImpl = rhsPtr->getEpetra_MultiVector();
    RCP<Epetra_MultiVector> lhsImpl       = this->getEpetra_MultiVector();

    TEUCHOS_TEST_FOR_EXCEPTION(
        rhsImpl.is_null(), std::logic_error,
        "Xpetra::MultiVector::operator= "
        "(in Xpetra::EpetraMultiVectorT::assign): *this (the right-hand side of "
        "the assignment) has a null RCP<Epetra_MultiVector> inside.  Please "
        "report this bug to the Xpetra developers.");
    TEUCHOS_TEST_FOR_EXCEPTION(
        lhsImpl.is_null(), std::logic_error,
        "Xpetra::MultiVector::operator= "
        "(in Xpetra::EpetraMultiVectorT::assign): The left-hand side of the "
        "assignment has a null RCP<Epetra_MultiVector> inside.  Please report "
        "this bug to the Xpetra developers.");

    // Epetra_MultiVector's assignment operator does a deep copy.
    *lhsImpl = *rhsImpl;
  }

 private:
  RCP<Epetra_MultiVector> vec_;

};  // EpetraMultiVectorT class (specialization on GO=long long, NO=EpetraNode
#endif

}  // namespace Xpetra

#include "Xpetra_EpetraVector.hpp"  // to avoid incomplete type instantiated above in out-of-body functions.

#endif  // XPETRA_EPETRAMULTIVECTOR_HPP