Domi_MDVector.hpp

// @HEADER
// ***********************************************************************
//
//     Domi: Multi-dimensional Distributed Linear Algebra Services
//                 Copyright (2014) 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 William F. Spotz (wfspotz@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef DOMI_MDVECTOR_HPP
#define DOMI_MDVECTOR_HPP

// #define DOMI_MDVECTOR_VERBOSE

// Standard includes
#include <ctime>

// Domi includes
#include "Domi_ConfigDefs.hpp"
#include "Domi_MDMap.hpp"
#include "Domi_MDArrayRCP.hpp"

// Teuchos includes
#include "Teuchos_DataAccess.hpp"
#include "Teuchos_Describable.hpp"
#include "Teuchos_ScalarTraitsDecl.hpp"
#include "Teuchos_Comm.hpp"
#include "Teuchos_CommHelpers.hpp"

#ifdef HAVE_EPETRA
#include "Epetra_IntVector.h"
#include "Epetra_Vector.h"
#endif

#ifdef HAVE_TPETRA
#include "Tpetra_Vector.hpp"
#endif


#ifdef OPEN_MPI
// I provide dummy definitions for MPI structs so that
// typeid(MPI_Datatype) and typeid(MPI_Request) will compile.
// Defining empty structs like this should be fine since only the guts
// of OpenMPI will ever see the real definitions.  This is a silly game
// we are playing with the C++ type system here but it should work
// just fine.
struct ompi_datatype_t {};
//struct ompi_request_t {};
#endif

#include <iostream>
using std::cout;
using std::endl;


namespace Domi
{

template< class Scalar >
class MDVector : public Teuchos::Describable
{
public:


  MDVector(const Teuchos::RCP< const MDMap > & mdMap,
           bool zeroOut = true);

  MDVector(const Teuchos::RCP< const MDMap > & mdMap,
           const dim_type leadingDim,
           const dim_type trailingDim = 0,
           bool zeroOut = true);

  MDVector(const Teuchos::RCP< const MDMap > & mdMap,
           const MDArrayView< Scalar > & source);

  MDVector(const Teuchos::RCP< const MDMap > & mdMap,
           const MDArrayRCP< Scalar > & source);

  MDVector(const MDVector< Scalar > & source,
           Teuchos::DataAccess access = Teuchos::View);

  MDVector(const Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm,
           Teuchos::ParameterList & plist);

  MDVector(const Teuchos::RCP< const MDComm > mdComm,
           Teuchos::ParameterList & plist);

  MDVector(const MDVector< Scalar > & parent,
           int axis,
           dim_type index);

  MDVector(const MDVector< Scalar > & parent,
           int                              axis,
           const Slice &                    slice,
           int                              bndryPad = 0);

  MDVector(const MDVector< Scalar > &    parent,
           const Teuchos::ArrayView< Slice > & slices,
           const Teuchos::ArrayView< int > &   bndryPad =
             Teuchos::ArrayView< int >());

  MDVector< Scalar > &
  operator=(const MDVector< Scalar > & source);

  virtual ~MDVector();



  inline const Teuchos::RCP< const MDMap >
  getMDMap() const;

  inline bool onSubcommunicator() const;

  inline Teuchos::RCP< const Teuchos::Comm< int > > getTeuchosComm() const;

  inline int numDims() const;

  inline int getCommDim(int axis) const;

  inline bool isPeriodic(int axis) const;

  inline int getCommIndex(int axis) const;

  inline int getLowerNeighbor(int axis) const;

  inline int getUpperNeighbor(int axis) const;

  inline dim_type getGlobalDim(int axis, bool withBndryPad=false) const;

  inline Slice getGlobalBounds(int axis, bool withBndryPadding=false) const;

  inline Slice getGlobalRankBounds(int axis, bool withBndryPad=false) const;

  inline dim_type getLocalDim(int axis, bool withCommPad=false) const;

  inline Teuchos::ArrayView< const Slice > getLocalBounds() const;

  inline Slice getLocalBounds(int axis, bool withPad=false) const;

  inline Slice getLocalInteriorBounds(int axis) const;

  inline bool hasPadding() const;

  inline int getLowerPadSize(int axis) const;

  inline int getUpperPadSize(int axis) const;

  inline int getCommPadSize(int axis) const;

  inline int getLowerBndryPad(int axis) const;

  inline int getUpperBndryPad(int axis) const;

  inline int getBndryPadSize(int axis) const;

  void setLowerPad(int axis, const Scalar value);

  void setUpperPad(int axis, const Scalar value);

  inline bool isReplicatedBoundary(int axis) const;

  inline Layout getLayout() const;

  inline bool isContiguous() const;



#ifdef HAVE_EPETRA

  Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorView() const;

  Teuchos::RCP< Epetra_Vector > getEpetraVectorView() const;

  Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorView() const;

  Teuchos::RCP< Epetra_IntVector > getEpetraIntVectorCopy() const;

  Teuchos::RCP< Epetra_Vector > getEpetraVectorCopy() const;

  Teuchos::RCP< Epetra_MultiVector > getEpetraMultiVectorCopy() const;

#endif

#ifdef HAVE_TPETRA

  template< class LocalOrdinal,
            class GlobalOrdinal = TpetraGOType>
  Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal > >
  getTpetraVectorView() const;

  template < class LocalOrdinal,
             class GlobalOrdinal = TpetraGOType>
  Teuchos::RCP< Tpetra::MultiVector< Scalar,
                                     LocalOrdinal,
                                     GlobalOrdinal > >
  getTpetraMultiVectorView() const;

  template< class LocalOrdinal,
            class GlobalOrdinal = TpetraGOType>
  Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal > >
  getTpetraVectorCopy() const;

  template < class LocalOrdinal,
             class GlobalOrdinal = TpetraGOType>
  Teuchos::RCP< Tpetra::MultiVector< Scalar,
                                     LocalOrdinal,
                                     GlobalOrdinal > >
  getTpetraMultiVectorCopy() const;

#endif



  MDArrayView< Scalar > getDataNonConst(bool includePadding = true);

  MDArrayView< const Scalar > getData(bool includePadding = true) const;

  MDArrayView< Scalar > getLowerPadDataNonConst(int axis);

  MDArrayView< const Scalar > getLowerPadData(int axis) const;

  MDArrayView< Scalar > getUpperPadDataNonConst(int axis);

  MDArrayView< const Scalar > getUpperPadData(int axis) const;



  Scalar
  dot(const MDVector< Scalar > & a) const;

  typename Teuchos::ScalarTraits< Scalar >::magnitudeType norm1() const;

  typename Teuchos::ScalarTraits< Scalar >::magnitudeType norm2() const;

  typename Teuchos::ScalarTraits< Scalar >::magnitudeType normInf() const;

  typename Teuchos::ScalarTraits< Scalar >::magnitudeType
  normWeighted(const MDVector< Scalar > & weights) const;

  Scalar meanValue() const;



  virtual std::string description() const;

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



  void putScalar(const Scalar & value,
                 bool includePadding = true);

  void randomize();



  // /** \brief Sum values of a locally replicated multivector across all
  //  *         processes.
  //  */
  // void reduce();

  void updateCommPad();

  /* \brief Update the data in the communication padding along the
   *        specified axis
   *
   * \param axis [in] the axis along which communication will be
   *        performed
   */
  void updateCommPad(int axis);

  void startUpdateCommPad(int axis);

  void endUpdateCommPad(int axis);



  MDVector< Scalar >
  operator[](dim_type index) const;

  MDVector< Scalar >
  operator[](Slice slice) const;



  void writeBinary(const std::string & filename,
                   bool includeBndryPad = false) const;

  void readBinary(const std::string & filename,
                  bool includeBndryPad = false);


private:

  // The Teuchos communicator.  Note that this is always a reference
  // to the communicator of the _mdMap, and is stored only for
  // convenience
  Teuchos::RCP< const Teuchos::Comm< int > > _teuchosComm;

  // The MDMap that describes the domain decomposition of this
  // MDVector
  Teuchos::RCP< const MDMap > _mdMap;

  // The MDArrayRCP that stores the data of this MDVector
  MDArrayRCP< Scalar > _mdArrayRcp;

  // The MDArrayView of the (possibly whole) sub-view into this
  // MDVector's MDArrayRCP
  MDArrayView< Scalar > _mdArrayView;

  // The operator[](int) and operator[](Slice) methods are applied to
  // a specific axis, namely this internally stored and updated next
  // axis
  int _nextAxis;

  // *** Communication Support *** //

#ifdef HAVE_MPI
  // An array of MPI_Request objects for supporting non-blocking sends
  // and receives
  Teuchos::Array< MPI_Request > _requests;
#endif

  // Define a struct for storing all the information needed for a
  // single message: a pointer to the buffer, a reference to the
  // message's MPI data type, the rank of the communication partner,
  // and the axis alng which we are communicating
  struct MessageInfo
  {
    // Pointer to first element of data buffer
    void *buffer;
#ifdef HAVE_MPI
    // MPI data type (strided vector)
    Teuchos::RCP< MPI_Datatype > datatype;
#else
    // Teuchos ArrayView for periodic domains
    MDArrayView< Scalar > dataview;
#endif
    // Processor rank for communication partner
    int proc;
    // Communication is along this axis
    int axis;
  };

  // An array of MessageInfo objects representing all active send
  // buffers.  The outer array represents the axis and the inner
  // 2-Tuple represents the lower and upper boundaries.
  Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _sendMessages;

  // An array of MessageInfo objects representing all active receive
  // buffers.  The outer array represents the axis and the inner
  // 2-Tuple represents the lower and upper boundaries.
  Teuchos::Array< Teuchos::Tuple< MessageInfo, 2 > > _recvMessages;

  // A private method to initialize the _sendMessages and
  // _recvMessages arrays.
  void initializeMessages();

  // *** Input/Output Support *** //

  // Define a struct for storing all of the information needed to
  // write or read the MDVector to a file: arrays that store the file
  // shape, the local buffer shape, the local data shape, the file
  // starting coordinates, and the data starting coordinates.
  struct FileInfo
  {
    Teuchos::Array< dim_type > fileShape;
    Teuchos::Array< dim_type > bufferShape;
    Teuchos::Array< dim_type > dataShape;
    Teuchos::Array< dim_type > fileStart;
    Teuchos::Array< dim_type > dataStart;
#ifdef HAVE_MPI
    Teuchos::RCP< MPI_Datatype > filetype;
    Teuchos::RCP< MPI_Datatype > datatype;
#endif
  };

  // FileInfo struct for an input or output file that does not store
  // boundary padding.  This is mutable because it does not get set
  // until the first time the MDVector is read or written to a file,
  // and the writeBinary() method should logically be const.
  mutable Teuchos::RCP< FileInfo > _fileInfo;

  // FileInfo struct for an input or output file that does store
  // boundary padding.  This is mutable because it does not get set
  // until the first time the MDVector is read or written to a file,
  // and the writeBinary() method should logically be const.
  mutable Teuchos::RCP< FileInfo > _fileInfoWithBndry;

  // Compute either the _fileInfo or _fileInfoWithBndry data members.
  // This private method gets called by the writeBinary() and
  // readBinary() methods, and sets the requested fileInfo object,
  // unless it has already been set.  This is const so that it can be
  // called by writeBinary(), but its whole purpose is to change
  // mutable data members.
  Teuchos::RCP< FileInfo > & computeFileInfo(bool includeBndryPad) const;

};

// *** Implementations *** //


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const MDMap > & mdMap,
         bool zeroOut) :
  _teuchosComm(mdMap->getTeuchosComm()),
  _mdMap(mdMap),
  _mdArrayRcp(),
  _mdArrayView(),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  // Obtain the array of dimensions
  int numDims = _mdMap->numDims();
  Teuchos::Array< dim_type > dims(numDims);
  for (int axis = 0; axis < numDims; ++axis)
    dims[axis] = _mdMap->getLocalDim(axis,true);

  // Resize the MDArrayRCP and set the MDArrayView
  _mdArrayRcp.resize(dims);
  _mdArrayView = _mdArrayRcp();
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const MDMap > & mdMap,
         const dim_type leadingDim,
         const dim_type trailingDim,
         bool zeroOut) :
  _teuchosComm(mdMap->getTeuchosComm()),
  _mdMap(mdMap->getAugmentedMDMap(leadingDim, trailingDim)),
  _mdArrayRcp(),
  _mdArrayView(),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  // Obtain the array of dimensions
  int numDims = _mdMap->numDims();
  Teuchos::Array< dim_type > dims(numDims);
  for (int axis = 0; axis < numDims; ++axis)
    dims[axis] = _mdMap->getLocalDim(axis,true);

  // Resize the MDArrayRCP and set the MDArrayView
  _mdArrayRcp.resize(dims);
  _mdArrayView = _mdArrayRcp();
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const MDMap > & mdMap,
         const MDArrayView< Scalar > & source) :
  _mdMap(mdMap),
  _mdArrayRcp(source),
  _mdArrayView(_mdArrayRcp()),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");
  int numDims = _mdMap->numDims();
  TEUCHOS_TEST_FOR_EXCEPTION(
    numDims != _mdArrayRcp.numDims(),
    InvalidArgument,
    "MDMap and source array do not have the same number of dimensions");

  for (int axis = 0; axis < numDims; ++axis)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
      _mdMap->getLocalDim(axis) != _mdArrayRcp.dimension(axis),
      InvalidArgument,
      "Axis " << axis << ": MDMap dimension = " << _mdMap->getLocalDim(axis)
      << ", MDArray dimension = " << _mdArrayRcp.dimension(axis));
  }
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const MDMap > & mdMap,
         const MDArrayRCP< Scalar > & mdArrayRcp) :
  _mdMap(mdMap),
  _mdArrayRcp(mdArrayRcp),
  _mdArrayView(_mdArrayRcp()),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
#ifdef DOMI_MDVECTOR_VERBOSE
  cout << "_mdArrayRcp  = " << _mdArrayRcp  << endl;
  cout << "_mdArrayView.getRawPtr() = " << _mdArrayView.getRawPtr()
            << " (constructor)" << endl;
  cout << "_mdArrayView = " << _mdArrayView << endl;
#endif
  setObjectLabel("Domi::MDVector");
  int numDims = _mdMap->numDims();
  TEUCHOS_TEST_FOR_EXCEPTION(
    numDims != _mdArrayRcp.numDims(),
    InvalidArgument,
    "MDMap and source array do not have the same number of dimensions");

  for (int axis = 0; axis < numDims; ++axis)
  {
    TEUCHOS_TEST_FOR_EXCEPTION(
      _mdMap->getLocalDim(axis) != _mdArrayRcp.dimension(axis),
      InvalidArgument,
      "Axis " << axis << ": MDMap dimension = " << _mdMap->getLocalDim(axis)
      << ", MDArray dimension = " << _mdArrayRcp.dimension(axis));
  }
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const MDVector< Scalar > & source,
         Teuchos::DataAccess access) :
  _teuchosComm(source.getMDMap()->getTeuchosComm()),
  _mdMap(source.getMDMap()),
  _mdArrayRcp(source._mdArrayRcp),
  _mdArrayView(source._mdArrayView),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  if (access == Teuchos::Copy)
  {
#ifdef DOMI_MDVECTOR_VERBOSE
    cout << "Inside MDVector copy constructor with copy access" << endl;
#endif
    // Obtain the array of dimensions
    int numDims = _mdMap->numDims();
    Teuchos::Array< dim_type > dims(numDims);
    for (int axis = 0; axis < numDims; ++axis)
      dims[axis] = _mdMap->getLocalDim(axis,true);

    // Reset the MDArrayRCP and set the MDArrayView
    _mdArrayRcp  = MDArrayRCP< Scalar >(dims, 0, source.getLayout());
    _mdArrayView = _mdArrayRcp();

    // Copy the source data to the new MDVector
    typedef typename MDArrayView< Scalar >::iterator iterator;
    typedef typename MDArrayView< Scalar >::const_iterator const_iterator;
    const_iterator src = source.getData().begin();
    for (iterator trg = _mdArrayView.begin(); trg != _mdArrayView.end(); ++trg)
    {
      *trg = *src;
      ++src;
    }
  }
#ifdef DOMI_MDVECTOR_VERBOSE
  else
  {
    cout << "Inside MDVector copy constructor with view access"
              << endl;
  }
#endif
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const Teuchos::Comm< int > > teuchosComm,
         Teuchos::ParameterList & plist) :
  _teuchosComm(teuchosComm),
  _mdMap(),
  _mdArrayRcp(),
  _mdArrayView(),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  // Compute the MDComm and MDMap
  MDMap * myMdMap = new MDMap(teuchosComm, plist);
  dim_type leadingDim  = plist.get("leading dimension" , 0);
  dim_type trailingDim = plist.get("trailing dimension", 0);
  if (leadingDim + trailingDim > 0)
  {
    _mdMap = myMdMap->getAugmentedMDMap(leadingDim, trailingDim);
    delete myMdMap;
  }
  else
    _mdMap = Teuchos::rcp(myMdMap);

  // Obtain the array of dimensions
  int numDims = _mdMap->numDims();
  Teuchos::Array< dim_type > dims(numDims);
  for (int axis = 0; axis < numDims; ++axis)
    dims[axis] = _mdMap->getLocalDim(axis,true);

  // Resize the MDArrayRCP and set the MDArrayView
  _mdArrayRcp.resize(dims);
  _mdArrayView = _mdArrayRcp();
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const Teuchos::RCP< const MDComm > mdComm,
         Teuchos::ParameterList & plist) :
  _teuchosComm(mdComm->getTeuchosComm()),
  _mdMap(Teuchos::rcp(new MDMap(mdComm, plist))),
  _mdArrayRcp(),
  _mdArrayView(),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  // Compute the MDMap
  MDMap * myMdMap = new MDMap(mdComm, plist);
  dim_type leadingDim  = plist.get("leading dimension" , 0);
  dim_type trailingDim = plist.get("trailing dimension", 0);
  if (leadingDim + trailingDim > 0)
  {
    _mdMap = myMdMap->getAugmentedMDMap(leadingDim, trailingDim);
    delete myMdMap;
  }
  else
    _mdMap = Teuchos::rcp(myMdMap);

  // Obtain the array of dimensions
  int numDims = _mdMap->numDims();
  Teuchos::Array< dim_type > dims(numDims);
  for (int axis = 0; axis < numDims; ++axis)
    dims[axis] = _mdMap->getLocalDim(axis,true);

  // Resize the MDArrayRCP and set the MDArrayView
  _mdArrayRcp.resize(dims);
  _mdArrayView = _mdArrayRcp();
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const MDVector< Scalar > & parent,
         int axis,
         dim_type globalIndex) :
  _teuchosComm(parent._teuchosComm),
  _mdMap(),
  _mdArrayRcp(parent._mdArrayRcp),
  _mdArrayView(parent._mdArrayView),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
  setObjectLabel("Domi::MDVector");

  // Obtain the parent MDMap
  Teuchos::RCP< const MDMap > parentMdMap = parent.getMDMap();

  // Obtain the new, sliced MDMap
  _mdMap = Teuchos::rcp(new MDMap(*parentMdMap,
                                          axis,
                                          globalIndex));

  // Check that we are on the new sub-communicator
  if (_mdMap->onSubcommunicator())
  {
    // Convert the index from global to local.  We start by
    // determining the starting global index on this processor along
    // the given axis, ignoring the boundary padding.  We then
    // subtract the lower padding, whether it is communication padding
    // or boundary padding.
    dim_type origin = parentMdMap->getGlobalRankBounds(axis,false).start() -
                      parentMdMap->getLowerPadSize(axis);

    // The local index along the given axis is the global axis minus
    // the starting index.  Since we are on the subcommunicator, this
    // should be valid.
    dim_type localIndex = globalIndex - origin;

    // Obtain the new MDArrayView using the local index
    MDArrayView< Scalar > newView(_mdArrayView, axis, localIndex);
    _mdArrayView = newView;
  }
  else
  {
    // We are not on the sub-communicator, so clear out the data
    // buffer and view
    _mdArrayRcp.clear();
    _mdArrayView = MDArrayView< Scalar >();
  }
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const MDVector< Scalar > & parent,
         int axis,
         const Slice & slice,
         int bndryPad) :
  _teuchosComm(),
  _mdMap(),
  _mdArrayRcp(parent._mdArrayRcp),
  _mdArrayView(parent._mdArrayView),
  _nextAxis(0),
#ifdef HAVE_MPI
  _requests(),
#endif
  _sendMessages(),
  _recvMessages()
{
#ifdef DOMI_MDVECTOR_VERBOSE
  cout << "slice axis " << axis << endl;
  cout << "  _mdArrayRcp  @ " << _mdArrayRcp.getRawPtr()  << endl;
  cout << "  _mdArrayView @ " << _mdArrayView.getRawPtr() << endl;
#endif

  setObjectLabel("Domi::MDVector");

  // Obtain the parent MDMap
  Teuchos::RCP< const MDMap > parentMdMap = parent.getMDMap();

  // Obtain the new, sliced MDMap
  _mdMap = Teuchos::rcp(new MDMap(*parentMdMap,
                                          axis,
                                          slice,
                                          bndryPad));
  _teuchosComm = _mdMap->getTeuchosComm();

  // Check that we are on the new sub-communicator
  if (_mdMap->onSubcommunicator())
  {
    // Get the concrete bounds
    Slice bounds = slice.bounds(parentMdMap->getGlobalDim(axis,true));

    // Convert the given Slice start index from global to local.  We
    // start by determining the starting global index on this
    // processor along the given axis, ignoring the boundary padding.
    // We then subtract the lower padding, whether it is communication
    // padding or boundary padding.
    dim_type origin = parentMdMap->getGlobalRankBounds(axis,false).start() -
                      parentMdMap->getLowerPadSize(axis);

    // Determine the starting index of our local slice.  This will be
    // the start of the given slice minus the starting global index on
    // this processor minus the given boundary pad.  If this is less
    // than zero, then the start is on a lower processor, so set the
    // local start to zero.
    dim_type start = std::max(0, bounds.start() - origin - bndryPad);

    // Now get the stop index of the local slice.  This will be the
    // stop of the given slice minus the starting global index on this
    // processor plus the given boundary pad.  If this is larger than
    // the local dimension, then set the local stop to the local
    // dimension.
    dim_type stop = std::min(bounds.stop() - origin + bndryPad,
                             parentMdMap->getLocalDim(axis,true));

    // Obtain the new MDArrayView using the local slice
    MDArrayView< Scalar > newView(_mdArrayView, axis, Slice(start,stop));
    _mdArrayView = newView;
  }
  else
  {
    // We are not on the sub-communicator, so clear out the data
    // buffer and view
    _mdArrayRcp.clear();
    _mdArrayView = MDArrayView< Scalar >();
  }
#ifdef DOMI_MDVECTOR_VERBOSE
  cout << "  _mdArrayView @ " << _mdArrayView.getRawPtr() << endl;
  cout << "  offset = " << int(_mdArrayView.getRawPtr() -
                                    _mdArrayRcp.getRawPtr()) << endl;
#endif
}


template< class Scalar >
MDVector< Scalar >::
MDVector(const MDVector< Scalar > & parent,
         const Teuchos::ArrayView< Slice > & slices,
         const Teuchos::ArrayView< int > & bndryPad)
{
  setObjectLabel("Domi::MDVector");

  // Temporarily store the number of dimensions
  int numDims = parent.numDims();

  // Sanity check on dimensions
  TEUCHOS_TEST_FOR_EXCEPTION(
    (slices.size() != numDims),
    InvalidArgument,
    "number of slices = " << slices.size() << " != parent MDVector number of "
    "dimensions = " << numDims);

  // Apply the single-Slice constructor to each axis in succession
  MDVector< Scalar > tempMDVector1(parent);
  for (int axis = 0; axis < numDims; ++axis)
  {
    int bndryPadding = (axis < bndryPad.size()) ? bndryPad[axis] : 0;
    MDVector< Scalar > tempMDVector2(tempMDVector1,
                                           axis,
                                           slices[axis],
                                           bndryPadding);
    tempMDVector1 = tempMDVector2;
  }
  *this = tempMDVector1;
}


template< class Scalar >
MDVector< Scalar > &
MDVector< Scalar >::
operator=(const MDVector< Scalar > & source)
{
  _teuchosComm  = source._teuchosComm;
  _mdMap        = source._mdMap;
  _mdArrayRcp   = source._mdArrayRcp;
  _mdArrayView  = source._mdArrayView;
  _nextAxis     = source._nextAxis;
#ifdef HAVE_MPI
  _requests     = source._requests;
#endif
  _sendMessages = source._sendMessages;
  _recvMessages = source._recvMessages;
  return *this;
}


template< class Scalar >
MDVector< Scalar >::~MDVector()
{
}


template< class Scalar >
const Teuchos::RCP< const MDMap >
MDVector< Scalar >::
getMDMap() const
{
  return _mdMap;
}


template< class Scalar >
bool
MDVector< Scalar >::
onSubcommunicator() const
{
  return _mdMap->onSubcommunicator();
}


template< class Scalar >
Teuchos::RCP< const Teuchos::Comm< int > >
MDVector< Scalar >::
getTeuchosComm() const
{
  return _mdMap->getTeuchosComm();
}


template< class Scalar >
int
MDVector< Scalar >::
numDims() const
{
  return _mdMap->numDims();
}


template< class Scalar >
int
MDVector< Scalar >::
getCommDim(int axis) const
{
  return _mdMap->getCommDim(axis);
}


template< class Scalar >
bool
MDVector< Scalar >::
isPeriodic(int axis) const
{
  return _mdMap->isPeriodic(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getCommIndex(int axis) const
{
  return _mdMap->getCommIndex(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getLowerNeighbor(int axis) const
{
  return _mdMap->getLowerNeighbor(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getUpperNeighbor(int axis) const
{
  return _mdMap->getUpperNeighbor(axis);
}


template< class Scalar >
dim_type
MDVector< Scalar >::
getGlobalDim(int axis, bool withBndryPad) const
{
  return _mdMap->getGlobalDim(axis, withBndryPad);
}


template< class Scalar >
Slice
MDVector< Scalar >::
getGlobalBounds(int axis, bool withBndryPad) const
{
  return _mdMap->getGlobalBounds(axis, withBndryPad);
}


template< class Scalar >
Slice
MDVector< Scalar >::
getGlobalRankBounds(int axis, bool withBndryPad) const
{
  return _mdMap->getGlobalRankBounds(axis, withBndryPad);
}


template< class Scalar >
dim_type
MDVector< Scalar >::
getLocalDim(int axis, bool withCommPad) const
{
  return _mdMap->getLocalDim(axis, withCommPad);
}


template< class Scalar >
Teuchos::ArrayView< const Slice >
MDVector< Scalar >::
getLocalBounds() const
{
  return _mdMap->getLocalBounds();
}


template< class Scalar >
Slice
MDVector< Scalar >::
getLocalBounds(int axis, bool withCommPad) const
{
  return _mdMap->getLocalBounds(axis, withCommPad);
}


template< class Scalar >
Slice
MDVector< Scalar >::
getLocalInteriorBounds(int axis) const
{
  return _mdMap->getLocalInteriorBounds(axis);
}


template< class Scalar >
bool
MDVector< Scalar >::
hasPadding() const
{
  return _mdMap->hasPadding();
}


template< class Scalar >
int
MDVector< Scalar >::
getLowerPadSize(int axis) const
{
  return _mdMap->getLowerPadSize(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getUpperPadSize(int axis) const
{
  return _mdMap->getUpperPadSize(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getCommPadSize(int axis) const
{
  return _mdMap->getCommPadSize(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getLowerBndryPad(int axis) const
{
  return _mdMap->getLowerBndryPad(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getUpperBndryPad(int axis) const
{
  return _mdMap->getUpperBndryPad(axis);
}


template< class Scalar >
int
MDVector< Scalar >::
getBndryPadSize(int axis) const
{
  return _mdMap->getBndryPadSize(axis);
}


template< class Scalar >
void
MDVector< Scalar >::
setLowerPad(int axis,
            const Scalar value)
{
  MDArrayView< Scalar > lowerPad = getLowerPadDataNonConst(axis);
  lowerPad.assign(value);
}


template< class Scalar >
void
MDVector< Scalar >::
setUpperPad(int axis,
            const Scalar value)
{
  MDArrayView< Scalar > upperPad = getUpperPadDataNonConst(axis);
  upperPad.assign(value);
}


template< class Scalar >
bool
MDVector< Scalar >::
isReplicatedBoundary(int axis) const
{
  return _mdMap->isReplicatedBoundary(axis);
}


template< class Scalar >
Layout
MDVector< Scalar >::
getLayout() const
{
  return _mdMap->getLayout();
}


template< class Scalar >
bool
MDVector< Scalar >::
isContiguous() const
{
  return _mdMap->isContiguous();
}


#ifdef HAVE_EPETRA

// The getEpetraIntVectorView() method only makes sense for Scalar =
// int, because Epetra_IntVectors store data buffers of type int only.
// There is no convenient way to specialize just one (or a small
// handfull of) methods, instead you have to specialize the whole
// class.  So we allow getEpetraIntVectorView() to compile for any
// Scalar, but we will throw an exception if Scalar is not int.

template< class Scalar >
Teuchos::RCP< Epetra_IntVector >
MDVector< Scalar >::
getEpetraIntVectorView() const
{
  // Throw an exception if Scalar is not int
  TEUCHOS_TEST_FOR_EXCEPTION(
    typeid(Scalar) != typeid(int),
    TypeError,
    "MDVector is of scalar type '" << typeid(Scalar).name() << "', but "
    "Epetra_IntVector requires scalar type 'int'");

  // Throw an exception if this MDVector's MDMap is not contiguous
  TEUCHOS_TEST_FOR_EXCEPTION(
    !isContiguous(),
    MDMapNoncontiguousError,
    "This MDVector's MDMap is non-contiguous.  This can happen when you take "
    "a slice of a parent MDVector.");

  // Get the Epetra_Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(true);

  // Return the result.  We are changing a Scalar* to a double*, which
  // looks sketchy, but we have thrown an exception if Sca is not
  // double, so everything is kosher.
  void * buffer = (void*) _mdArrayView.getRawPtr();
  return Teuchos::rcp(new Epetra_IntVector(View,
                                           *epetraMap,
                                           (int*) buffer));
}


// The getEpetraVectorView() method only makes sense for Scalar =
// double, because Epetra_Vectors store data buffers of type double
// only.  There is no convenient way to specialize just one (or a
// small handfull of) methods, instead you have to specialize the
// whole class.  So we allow getEpetraVectorView() to compile for any
// Scalar, but we will throw an exception if Scalar is not double.

template< class Scalar >
Teuchos::RCP< Epetra_Vector >
MDVector< Scalar >::
getEpetraVectorView() const
{
  // Throw an exception if Scalar is not double
  TEUCHOS_TEST_FOR_EXCEPTION(
    typeid(Scalar) != typeid(double),
    TypeError,
    "MDVector is of scalar type '" << typeid(Scalar).name() << "', but "
    "Epetra_Vector requires scalar type 'double'");

  // Throw an exception if this MDVector's MDMap is not contiguous
  TEUCHOS_TEST_FOR_EXCEPTION(
    !isContiguous(),
    MDMapNoncontiguousError,
    "This MDVector's MDMap is non-contiguous.  This can happen when you take "
    "a slice of a parent MDVector.");

  // Get the Epetra_Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(true);

  // Return the result.  We are changing a Scalar* to a double*, which
  // looks sketchy, but we have thrown an exception if Sca is not
  // double, so everything is kosher.
  void * buffer = (void*) _mdArrayView.getRawPtr();
  return Teuchos::rcp(new Epetra_Vector(View,
                                        *epetraMap,
                                        (double*) buffer));
}


// The getEpetraMultiVectorView() method only makes sense for Scalar =
// double, because Epetra_MultiVectors store data buffers of type
// double only.  There is no convenient way to specialize just one (or
// a small handfull of) methods, instead you have to specialize the
// whole class.  So we allow getEpetraVectorView() to compile for any
// Scalar, but we will throw an exception if Scalar is not double.

template< class Scalar >
Teuchos::RCP< Epetra_MultiVector >
MDVector< Scalar >::
getEpetraMultiVectorView() const
{
  // Throw an exception if Scalar is not double
  TEUCHOS_TEST_FOR_EXCEPTION(
    typeid(Scalar) != typeid(double),
    TypeError,
    "MDVector is of scalar type '" << typeid(Scalar).name() << "', but "
    "Epetra_Vector requires scalar type 'double'");

  // Determine the vector axis and related info
  int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
  int padding    = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
  int commDim    = getCommDim(vectorAxis);
  int numVectors = getGlobalDim(vectorAxis);

  // Obtain the appropriate MDMap and check that it is contiguous
  Teuchos::RCP< const MDMap > newMdMap;
  if (padding == 0 && commDim == 1)
    newMdMap = Teuchos::rcp(new MDMap(*_mdMap, vectorAxis, 0));
  else
  {
    newMdMap = _mdMap;
    numVectors = 1;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    ! newMdMap->isContiguous(),
    MDMapNoncontiguousError,
    "This MDVector's MDMap is non-contiguous.  This can happen when you take "
    "a slice of a parent MDVector.");

  // Get the stride between vectors.  The MDMap strides are private,
  // but we know the new MDMap is contiguous, so we can calculate it
  // as the product of the new MDMap dimensions (including padding)
  size_type stride = newMdMap->getLocalDim(0,true);
  for (int axis = 1; axis < newMdMap->numDims(); ++axis)
    stride *= newMdMap->getLocalDim(axis,true);
  TEUCHOS_TEST_FOR_EXCEPTION(
    stride*numVectors > Teuchos::OrdinalTraits<int>::max(),
    MapOrdinalError,
    "Buffer size " << stride*numVectors << " is too large for Epetra int "
    "ordinals");
  int lda = (int)stride;

  // Get the Epetra_Map that is equivalent to newMdMap
  Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(true);

  // Return the result.  We are changing a Scalar* to a double*, which
  // looks sketchy, but we have thrown an exception if Sca is not
  // double, so everything is kosher.
  void * buffer = (void*) _mdArrayView.getRawPtr();
  return Teuchos::rcp(new Epetra_MultiVector(View,
                                             *epetraMap,
                                             (double*) buffer,
                                             lda,
                                             numVectors));
}


template< class Scalar >
Teuchos::RCP< Epetra_IntVector >
MDVector< Scalar >::
getEpetraIntVectorCopy() const
{
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  // Get the Epetra_Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(true);

  // Construct the result
  Teuchos::RCP< Epetra_IntVector > result =
    Teuchos::rcp(new Epetra_IntVector(*epetraMap));

  // Copy the MDVector data buffer to the Epetra_IntVector, even if the
  // MDVector is non-contiguous
  int ii = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    result->operator[](ii++) = (int) *it;

  // Return the result
  return result;
}


template< class Scalar >
Teuchos::RCP< Epetra_Vector >
MDVector< Scalar >::
getEpetraVectorCopy() const
{
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  // Get the Epetra_Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Epetra_Map > epetraMap = _mdMap->getEpetraMap(true);

  // Construct the result
  Teuchos::RCP< Epetra_Vector > result =
    Teuchos::rcp(new Epetra_Vector(*epetraMap));

  // Copy the MDVector data buffer to the Epetra_Vector, even if the
  // MDVector is non-contiguous
  int ii = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    result->operator[](ii++) = (double) *it;

  // Return the result
  return result;
}


template< class Scalar >
Teuchos::RCP< Epetra_MultiVector >
MDVector< Scalar >::
getEpetraMultiVectorCopy() const
{
  typedef typename MDArrayView< Scalar >::iterator iterator;
  typedef typename MDArrayView< const Scalar >::iterator citerator;

  // Determine the vector axis and related info
  int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
  int padding    = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
  int commDim    = getCommDim(vectorAxis);
  int numVectors = getGlobalDim(vectorAxis);

  // Obtain the appropriate MDMap
  Teuchos::RCP< const MDMap > newMdMap;
  if (padding == 0 && commDim == 1)
    newMdMap = Teuchos::rcp(new MDMap(*_mdMap, vectorAxis, 0));
  else
  {
    newMdMap = _mdMap;
    numVectors = 1;
  }

  // Get the Epetra_Map that is equivalent to newMdMap
  Teuchos::RCP< const Epetra_Map > epetraMap = newMdMap->getEpetraMap(true);

  // Construct the result
  Teuchos::RCP< Epetra_MultiVector > result =
    Teuchos::rcp(new Epetra_MultiVector(*epetraMap, numVectors));

  // Copy the MDVector data to the Epetra_MultiVector, even if the
  // MDVector is non-contiguous
  int ii = 0;
  if (numVectors == 1)
  {
    for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
      result->operator[](0)[ii++] = (double) *it;
  }
  else
  {
    for (int iv = 0; iv < numVectors; ++iv)
    {
      ii = 0;
      MDArrayView< Scalar > subArray(_mdArrayView, vectorAxis, iv);
      for (iterator it = subArray.begin(); it != subArray.end(); ++it)
        result->operator[](iv)[ii++] = (double) *it;
    }
  }

  // Return the result
  return result;
}

#endif


#ifdef HAVE_TPETRA

template< class Scalar >
template< class LocalOrdinal,
          class GlobalOrdinal >
Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal > >
MDVector< Scalar >::
getTpetraVectorView() const
{
  // Throw an exception if this MDVector's MDMap is not contiguous
  TEUCHOS_TEST_FOR_EXCEPTION(
    !isContiguous(),
    MDMapNoncontiguousError,
    "This MDVector's MDMap is non-contiguous.  This can happen when you take "
    "a slice of a parent MDVector.");

  // Get the Tpetra::Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Tpetra::Map< LocalOrdinal,
                                   GlobalOrdinal > > tpetraMap =
    _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal >(true);

  // Return the result
  return Teuchos::rcp(new Tpetra::Vector< Scalar,
                                          LocalOrdinal,
                                          GlobalOrdinal >(tpetraMap,
                                                          _mdArrayView.arrayView()));
}


template< class Scalar >
template< class LocalOrdinal,
          class GlobalOrdinal >
Teuchos::RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal > >
MDVector< Scalar >::
getTpetraVectorCopy() const
{
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  // Get the Tpetra::Map that is equivalent to this MDVector's MDMap
  Teuchos::RCP< const Tpetra::Map< LocalOrdinal,
                                   GlobalOrdinal > > tpetraMap =
    _mdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal >(true);

  // Construct the result
  Teuchos::RCP< Tpetra::Vector< Scalar,
                                LocalOrdinal,
                                GlobalOrdinal > > result =
    Teuchos::rcp(new Tpetra::Vector< Scalar,
                                     LocalOrdinal,
                                     GlobalOrdinal >(tpetraMap) );

  // Copy the MDVector data to the Tpetra::Vector, even if the
  // MDVector is non-contiguous
  Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst();
  int ii = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    tpetraVectorArray[ii++] = *it;

  // Return the result
  return result;
}


template< class Scalar >
template < class LocalOrdinal,
           class GlobalOrdinal >
Teuchos::RCP< Tpetra::MultiVector< Scalar,
                                   LocalOrdinal,
                                   GlobalOrdinal > >
MDVector< Scalar >::
getTpetraMultiVectorView() const
{
  // Determine the vector axis and related info
  int vectorAxis    = (getLayout() == C_ORDER) ? 0 : numDims()-1;
  int padding       = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
  int commDim       = getCommDim(vectorAxis);
  size_t numVectors = getGlobalDim(vectorAxis);

  // Obtain the appropriate MDMap and check that it is contiguous
  Teuchos::RCP< const MDMap > newMdMap;
  if (padding == 0 && commDim == 1)
    newMdMap = Teuchos::rcp(new MDMap(*_mdMap, vectorAxis, 0));
  else
  {
    newMdMap = _mdMap;
    numVectors = 1;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    ! newMdMap->isContiguous(),
    MDMapNoncontiguousError,
    "This MDVector's MDMap is non-contiguous.  This can happen when you take "
    "a slice of a parent MDVector.");

  // Get the stride between vectors.  The MDMap strides are private,
  // but we know the new MDMap is contiguous, so we can calculate it
  // as the product of the new MDMap dimensions (including padding)
  size_type stride = newMdMap->getLocalDim(0,true);
  for (int axis = 1; axis < newMdMap->numDims(); ++axis)
    stride *= newMdMap->getLocalDim(axis,true);
  TEUCHOS_TEST_FOR_EXCEPTION(
    (long long int)(stride*numVectors) > Teuchos::OrdinalTraits<GlobalOrdinal>::max(),
    MapOrdinalError,
    "Buffer size " << stride*numVectors << " is too large for Tpetra "
    "GlobalOrdinal = " << typeid(GlobalOrdinal).name() );
  size_t lda = (size_t)stride;

  // Get the Tpetra::Map that is equivalent to newMdMap
  Teuchos::RCP< const Tpetra::Map< LocalOrdinal,
                                   GlobalOrdinal > > tpetraMap =
    newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal >(true);

  // Return the result
  return Teuchos::rcp(new Tpetra::MultiVector< Scalar,
                                               LocalOrdinal,
                                               GlobalOrdinal >(tpetraMap,
                                                               _mdArrayView.arrayView(),
                                                               lda,
                                                               numVectors));
}


template< class Scalar >
template < class LocalOrdinal,
           class GlobalOrdinal >
Teuchos::RCP< Tpetra::MultiVector< Scalar,
                                   LocalOrdinal,
                                   GlobalOrdinal > >
MDVector< Scalar >::
getTpetraMultiVectorCopy() const
{
  typedef typename MDArrayView< Scalar >::iterator iterator;
  typedef typename MDArrayView< const Scalar >::iterator citerator;

  // Determine the vector axis and related info
  int vectorAxis = (getLayout() == C_ORDER) ? 0 : numDims()-1;
  int padding    = getLowerPadSize(vectorAxis) + getUpperPadSize(vectorAxis);
  int commDim    = getCommDim(vectorAxis);
  int numVectors = getGlobalDim(vectorAxis);

  // Obtain the appropriate MDMap
  Teuchos::RCP< const MDMap > newMdMap;
  if (padding == 0 && commDim == 1)
    newMdMap = Teuchos::rcp(new MDMap(*_mdMap, vectorAxis, 0));
  else
  {
    newMdMap = _mdMap;
    numVectors = 1;
  }

  // Get the Tpetra::Map that is equivalent to newMdMap
  Teuchos::RCP< const Tpetra::Map< LocalOrdinal,
                                   GlobalOrdinal > > tpetraMap =
    newMdMap->template getTpetraMap< LocalOrdinal, GlobalOrdinal >(true);

  // Construct the result
  Teuchos::RCP< Tpetra::MultiVector< Scalar,
                                     LocalOrdinal,
                                     GlobalOrdinal > > result =
    Teuchos::rcp(new Tpetra::MultiVector< Scalar,
                                          LocalOrdinal,
                                          GlobalOrdinal >(tpetraMap, numVectors));

  // Copy the MDVector to the Tpetra::MultiVector, even if the
  // MDVector is non-contiguous
  int ii = 0;
  if (numVectors == 1)
  {
    Teuchos::ArrayRCP< Scalar > tpetraVectorArray = result->getDataNonConst(0);
    for (citerator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
      tpetraVectorArray[ii++] = (double) *it;
  }
  else
  {
    for (int iv = 0; iv < numVectors; ++iv)
    {
      ii = 0;
      Teuchos::ArrayRCP< Scalar > tpetraVectorArray =
        result->getDataNonConst(iv);
      MDArrayView< Scalar > subArray(_mdArrayView, vectorAxis, iv);
      for (iterator it = subArray.begin(); it != subArray.end(); ++it)
        tpetraVectorArray[ii++] = (double) *it;
    }
  }

  // Return the result
  return result;
}

#endif


template< class Scalar >
MDArrayView< Scalar >
MDVector< Scalar >::
getDataNonConst(bool includePadding)
{
#ifdef DOMI_MDVECTOR_VERBOSE
  cout << "_mdArrayView.getRawPtr() = " << _mdArrayView.getRawPtr()
            << endl;
  cout << "_mdArrayView = " << _mdArrayView << endl;
#endif
  if (includePadding)
    return _mdArrayView;
  MDArrayView< Scalar > newArray(_mdArrayView);
  for (int axis = 0; axis < numDims(); ++axis)
  {
    int lo = getLowerPadSize(axis);
    int hi = getLocalDim(axis,true) - getUpperPadSize(axis);
    newArray = MDArrayView< Scalar >(newArray, axis, Slice(lo,hi));
  }
  return newArray;
}


template< class Scalar >
MDArrayView< const Scalar >
MDVector< Scalar >::
getData(bool includePadding) const
{
  if (includePadding)
    return _mdArrayView.getConst();
  MDArrayView< Scalar > newArray(_mdArrayView);
  for (int axis = 0; axis < numDims(); ++axis)
  {
    int lo = getLowerPadSize(axis);
    int hi = getLocalDim(axis,true) - getUpperPadSize(axis);
    newArray = MDArrayView< Scalar >(newArray, axis, Slice(lo,hi));
  }
  return newArray.getConst();
}


template< class Scalar >
MDArrayView< Scalar >
MDVector< Scalar >::
getLowerPadDataNonConst(int axis)
{
  MDArrayView< Scalar > newArrayView(_mdArrayView,
                                     axis,
                                     Slice(getLowerPadSize(axis)));
  return newArrayView;
}


template< class Scalar >
MDArrayView< const Scalar >
MDVector< Scalar >::
getLowerPadData(int axis) const
{
  MDArrayView< const Scalar > newArrayView(_mdArrayView.getConst(),
                                           axis,
                                           Slice(getLowerPadSize(axis)));
  return newArrayView;
}


template< class Scalar >
MDArrayView< Scalar >
MDVector< Scalar >::
getUpperPadDataNonConst(int axis)
{
  dim_type n  = getLocalDim(axis,true);
  int     pad = getUpperPadSize(axis);
  Slice   slice;
  if (pad) slice = Slice(n-pad,n);
  else     slice = Slice(n-1,n-1);
  MDArrayView< Scalar > newArrayView(_mdArrayView,
                                     axis,
                                     slice);
  return newArrayView;
}


template< class Scalar >
MDArrayView< const Scalar >
MDVector< Scalar >::
getUpperPadData(int axis) const
{
  dim_type n  = getLocalDim(axis,true);
  int     pad = getUpperPadSize(axis);
  Slice   slice;
  if (pad) slice = Slice(n-pad,n);
  else     slice = Slice(n-1,n-1);
  MDArrayView< const Scalar > newArrayView(_mdArrayView.getConst(),
                                           axis,
                                           slice);
  return newArrayView;
}


template< class Scalar >
Scalar
MDVector< Scalar >::
dot(const MDVector< Scalar > & a) const
{
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  TEUCHOS_TEST_FOR_EXCEPTION(
    ! _mdMap->isCompatible(*(a._mdMap)),
    MDMapError,
    "MDMap of calling MDVector and argument 'a' are incompatible");

  MDArrayView< const Scalar > aView = a.getData();
  Scalar local_dot = 0;
  iterator a_it = aView.begin();
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
       ++it, ++a_it)
    local_dot += *it * *a_it;
  Scalar global_dot = 0;
  Teuchos::reduceAll(*_teuchosComm,
                     Teuchos::REDUCE_SUM,
                     1,
                     &local_dot,
                     &global_dot);
  return global_dot;
}


template< class Scalar >
typename Teuchos::ScalarTraits< Scalar >::magnitudeType
MDVector< Scalar >::
norm1() const
{
  typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  mag local_norm1 = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    local_norm1 += std::abs(*it);
  mag global_norm1 = 0;
  Teuchos::reduceAll(*_teuchosComm,
                     Teuchos::REDUCE_SUM,
                     1,
                     &local_norm1,
                     &global_norm1);
  return global_norm1;
}


template< class Scalar >
typename Teuchos::ScalarTraits< Scalar >::magnitudeType
MDVector< Scalar >::
norm2() const
{
  typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
  mag norm2 = dot(*this);
  return Teuchos::ScalarTraits<mag>::squareroot(norm2);
}


template< class Scalar >
typename Teuchos::ScalarTraits< Scalar >::magnitudeType
MDVector< Scalar >::
normInf() const
{
  typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  mag local_normInf = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    local_normInf = std::max(local_normInf, std::abs(*it));
  mag global_normInf = 0;
  Teuchos::reduceAll(*_teuchosComm,
                     Teuchos::REDUCE_MAX,
                     1,
                     &local_normInf,
                     &global_normInf);
  return global_normInf;
}


template< class Scalar >
typename Teuchos::ScalarTraits< Scalar >::magnitudeType
MDVector< Scalar >::
normWeighted(const MDVector< Scalar > & weights) const
{
  typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  TEUCHOS_TEST_FOR_EXCEPTION(
    ! _mdMap->isCompatible(*(weights._mdMap)),
    MDMapError,
    "MDMap of calling MDVector and argument 'weights' are incompatible");

  MDArrayView< const Scalar > wView = weights.getData();
  mag local_wNorm = 0;
  iterator w_it = wView.begin();
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end();
       ++it, ++w_it)
    local_wNorm += *it * *it * *w_it;
  mag global_wNorm = 0;
  Teuchos::reduceAll(*_teuchosComm,
                     Teuchos::REDUCE_SUM,
                     1,
                     &local_wNorm,
                     &global_wNorm);
  Teuchos::Array< dim_type > dimensions(numDims());
  for (int i = 0; i < numDims(); ++i)
    dimensions[i] = _mdMap->getGlobalDim(i);
  size_type n = computeSize(dimensions);
  if (n == 0) return 0;
  return Teuchos::ScalarTraits<mag>::squareroot(global_wNorm / n);
}


template< class Scalar >
Scalar
MDVector< Scalar >::
meanValue() const
{
  typedef typename Teuchos::ScalarTraits< Scalar >::magnitudeType mag;
  typedef typename MDArrayView< const Scalar >::iterator iterator;

  mag local_sum = 0;
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    local_sum += *it;
  mag global_sum = 0;
  Teuchos::reduceAll(*_teuchosComm,
                     Teuchos::REDUCE_SUM,
                     1,
                     &local_sum,
                     &global_sum);
  Teuchos::Array< dim_type > dimensions(numDims());
  for (int i = 0; i < numDims(); ++i)
    dimensions[i] = _mdMap->getGlobalDim(i);
  size_type n = computeSize(dimensions);
  if (n == 0) return 0;
  return global_sum / n;
}


template< class Scalar >
std::string
MDVector< Scalar >::
description() const
{
  using Teuchos::TypeNameTraits;

  Teuchos::Array< dim_type > dims(numDims());
  for (int axis = 0; axis < numDims(); ++axis)
    dims[axis] = getGlobalDim(axis, true);

  std::ostringstream oss;
  oss << "\"Domi::MDVector\": {"
      << "Template parameters: {"
      << "Scalar: " << TypeNameTraits<Scalar>::name()
      << "}";
  if (this->getObjectLabel() != "")
    oss << ", Label: \"" << this->getObjectLabel () << "\", ";
  oss << "Global dimensions: " << dims << " }";
  return oss.str();
}


template< class Scalar >
void
MDVector< Scalar >::
describe(Teuchos::FancyOStream &out,
         const Teuchos::EVerbosityLevel verbLevel) const
{
  using std::setw;
  using Teuchos::Comm;
  using Teuchos::RCP;
  using Teuchos::TypeNameTraits;
  using Teuchos::VERB_DEFAULT;
  using Teuchos::VERB_NONE;
  using Teuchos::VERB_LOW;
  using Teuchos::VERB_MEDIUM;
  using Teuchos::VERB_HIGH;
  using Teuchos::VERB_EXTREME;

  const Teuchos::EVerbosityLevel vl =
    (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;

  const MDMap & mdMap = *(getMDMap());
  Teuchos::RCP< const Teuchos::Comm< int > > comm = mdMap.getTeuchosComm();
  const int myImageID = comm->getRank();
  const int numImages = comm->getSize();
  Teuchos::OSTab tab0(out);

  if (vl != VERB_NONE)
  {
    if (myImageID == 0)
    {
      out << "\"Domi::MDVector\":" << endl;
    }
    Teuchos::OSTab tab1(out);// applies to all processes
    if (myImageID == 0)
    {
      out << "Template parameters:";
      {
        Teuchos::OSTab tab2(out);
        out << "Scalar: " << TypeNameTraits<Scalar>::name() << endl;
      }
      out << endl;
      if (this->getObjectLabel() != "")
      {
        out << "Label: \"" << getObjectLabel() << "\"" << endl;
      }
      Teuchos::Array< dim_type > globalDims(numDims());
      for (int axis = 0; axis < numDims(); ++axis)
        globalDims[axis] = getGlobalDim(axis, true);
      out << "Global dimensions: " << globalDims << endl;
    }
    for (int imageCtr = 0; imageCtr < numImages; ++imageCtr)
    {
      if (myImageID == imageCtr)
      {
        if (vl != VERB_LOW)
        {
          // VERB_MEDIUM and higher prints getLocalLength()
          out << "Process: " << myImageID << endl;
          Teuchos::OSTab tab2(out);
          Teuchos::Array< dim_type > localDims(numDims());
          for (int axis = 0; axis < numDims(); ++axis)
            localDims[axis] = getLocalDim(axis, true);
          out << "Local dimensions: " << localDims << endl;
        }
        std::flush(out); // give output time to complete
      }
      comm->barrier(); // give output time to complete
      comm->barrier();
      comm->barrier();
    }
  }
}


template< class Scalar >
void
MDVector< Scalar >::
putScalar(const Scalar & value,
          bool includePadding)
{
  typedef typename MDArrayView< Scalar >::iterator iterator;

  MDArrayView< Scalar > newArray = getDataNonConst(includePadding);
  for (iterator it = newArray.begin(); it != newArray.end(); ++it)
    *it = value;
}


template< class Scalar >
void
MDVector< Scalar >::
randomize()
{
  typedef typename MDArrayView< Scalar >::iterator iterator;

  Teuchos::ScalarTraits< Scalar >::seedrandom(time(NULL));
  for (iterator it = _mdArrayView.begin(); it != _mdArrayView.end(); ++it)
    *it = Teuchos::ScalarTraits< Scalar >::random();
}


template< class Scalar >
void
MDVector< Scalar >::
updateCommPad()
{
  for (int axis = 0; axis < numDims(); ++axis)
  {
    updateCommPad(axis);
  }
}


template< class Scalar >
void
MDVector< Scalar >::
updateCommPad(int axis)
{
  startUpdateCommPad(axis);
  endUpdateCommPad(axis);
}


template< class Scalar >
void
MDVector< Scalar >::
startUpdateCommPad(int axis)
{
  // #define DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD

  // Initialize the _sendMessages and _recvMessages members on the
  // first call to startUpdateCommPad(int).
  if (_sendMessages.empty()) initializeMessages();

#ifdef HAVE_MPI
  int rank    = _teuchosComm->getRank();
  int numProc = _teuchosComm->getSize();
  int tag;
  // Since HAVE_MPI is defined, we know that _teuchosComm points to a
  // const Teuchos::MpiComm< int >.  We downcast, extract and
  // dereference so that we can get access to the MPI_Comm used to
  // construct it.
  Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
    Teuchos::rcp_dynamic_cast< const Teuchos::MpiComm< int > >(_teuchosComm);
  const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
    *(mpiComm->getRawMpiComm());

  // Post the non-blocking sends
  MPI_Request request;
  for (int boundary = 0; boundary < 2; ++boundary)
  {
    MessageInfo message = _sendMessages[axis][boundary];
    if (message.proc >= 0)
    {
      tag = 2 * (rank * numProc + message.proc) + boundary;

#ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD
      cout << rank << ": post send for axis " << axis << ", boundary "
           << boundary << ", buffer = " << message.buffer << ", proc = "
           << message.proc << ", tag = " << tag << endl;
#endif

      if (MPI_Isend(message.buffer,
                    1,
                    *(message.datatype),
                    message.proc,
                    tag,
                    communicator(),
                    &request))
        throw std::runtime_error("Domi::MDVector: Error in MPI_Isend");
      _requests.push_back(request);
    }
  }

  // Post the non-blocking receives
  for (int boundary = 0; boundary < 2; ++boundary)
  {
    MessageInfo message = _recvMessages[axis][boundary];
    if (message.proc >= 0)
    {
      tag = 2 * (message.proc * numProc + rank) + (1-boundary);

#ifdef DOMI_MDVECTOR_OUTPUT_UPDATECOMMPAD
      cout << rank << ": post recv for axis " << axis << ", boundary "
           << boundary << ", buffer = " << message.buffer << ", proc = "
           << message.proc << ", tag = " << tag << endl;
#endif

      if (MPI_Irecv(message.buffer,
                    1,
                    *(message.datatype),
                    message.proc,
                    tag,
                    communicator(),
                    &request))
        throw std::runtime_error("Domi::MDVector: Error in MPI_Irecv");
      _requests.push_back(request);
    }
  }
#else
  // HAVE_MPI is not defined, so we are on a single processor.
  // However, if the axis is periodic, we need to copy the appropriate
  // data to the communication padding.
  if (isPeriodic(axis))
  {
    for (int sendBndry = 0; sendBndry < 2; ++sendBndry)
    {
      int recvBndry = 1 - sendBndry;
      // Get the receive and send data views
      MDArrayView< Scalar > recvView = _recvMessages[axis][recvBndry].dataview;
      MDArrayView< Scalar > sendView = _sendMessages[axis][sendBndry].dataview;

      // Initialize the receive and send data view iterators
      typename MDArrayView< Scalar >::iterator it_recv = recvView.begin();
      typename MDArrayView< Scalar >::iterator it_send = sendView.begin();

      // Copy the send buffer to the receive buffer
      for ( ; it_recv != recvView.end(); ++it_recv, ++it_send)
        *it_recv = *it_send;
    }
  }
#endif
}


template< class Scalar >
void
MDVector< Scalar >::
endUpdateCommPad(int axis)
{
#ifdef HAVE_MPI
  if (_requests.size() > 0)
  {
    Teuchos::Array< MPI_Status > status(_requests.size());
    if (MPI_Waitall(_requests.size(),
                    &(_requests[0]),
                    &(status[0]) ) )
      throw std::runtime_error("Domi::MDVector: Error in MPI_Waitall");
    _requests.clear();
  }
#endif
}


template< class Scalar >
MDVector< Scalar >
MDVector< Scalar >::
operator[](dim_type index) const
{
  MDVector< Scalar > result(*this,
                                  _nextAxis,
                                  index);
  int newAxis = _nextAxis + 1;
  if (newAxis >= result.numDims()) newAxis = 0;
  result._nextAxis = newAxis;
  return result;
}


template< class Scalar >
MDVector< Scalar >
MDVector< Scalar >::
operator[](Slice slice) const
{
  MDVector< Scalar > result(*this,
                                  _nextAxis,
                                  slice     );
  int newAxis = _nextAxis + 1;
  if (newAxis >= result.numDims()) newAxis = 0;
  result._nextAxis = newAxis;
  return result;
}


template< class Scalar >
void
MDVector< Scalar >::
initializeMessages()
{
  // #define DOMI_MDVECTOR_MESSAGE_INITIALIZE

  int ndims = numDims();
  Teuchos::Array<int> sizes(ndims);
  Teuchos::Array<int> subsizes(ndims);
  Teuchos::Array<int> starts(ndims);
  MessageInfo messageInfo;

  _sendMessages.resize(ndims);
  _recvMessages.resize(ndims);

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
  std::stringstream msg;
  int rank = _teuchosComm->getRank();
#endif

#ifdef HAVE_MPI
  int order = mpiOrder(getLayout());
  MPI_Datatype datatype = mpiType< Scalar >();
#endif

  // Loop over the axes we are going to send messages along
  for (int msgAxis = 0; msgAxis < ndims; ++msgAxis)
  {
    // Set the initial values for sizes, subsizes and starts
    for (int axis = 0; axis < ndims; ++axis)
    {
      sizes[axis]    = _mdArrayRcp.dimension(axis);
      subsizes[axis] = _mdArrayView.dimension(axis);
      starts[axis]   = 0;
    }

    // Set the lower receive and send messages

    int proc = getLowerNeighbor(msgAxis);

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
    msg << endl << "P" << rank << ": axis " << msgAxis << ", lower neighbor = "
        << proc << endl;
#endif

    // Fix the subsize along this message axis
    subsizes[msgAxis] = getLowerPadSize(msgAxis);
    // Fix the proc if the subsize is zero
    if (subsizes[msgAxis] == 0) proc = -1;
    // Assign the non-MPI members of messageInfo
    messageInfo.buffer = (void*) getData().getRawPtr();
    messageInfo.proc   = proc;
    messageInfo.axis   = msgAxis;

    if (proc >= 0)
    {
      // Lower receive message

#ifdef HAVE_MPI

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
      msg << endl << "P" << rank << ": axis " << msgAxis
          << ", Lower receive message:" << endl << "  ndims    = " << ndims
          << endl << "  sizes    = " << sizes << endl << "  subsizes = "
          << subsizes << endl << "  starts   = " << starts << endl
          << "  order    = " << order << endl;
#endif
      Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(new MPI_Datatype);
      MPI_Type_create_subarray(ndims,
                               &sizes[0],
                               &subsizes[0],
                               &starts[0],
                               order,
                               datatype,
                               commPad.get());
      MPI_Type_commit(commPad.get());
      messageInfo.datatype = commPad;
#else
      messageInfo.dataview = _mdArrayView;
      for (int axis = 0; axis < numDims(); ++axis)
      {
        Slice slice(starts[axis], starts[axis] + subsizes[axis]);
        messageInfo.dataview = MDArrayView< Scalar >(messageInfo.dataview,
                                                     axis,
                                                     slice);
      }
#endif

    }
    _recvMessages[msgAxis][0] = messageInfo;

    // Lower send message

    starts[msgAxis] = getLowerPadSize(msgAxis);
    if (isReplicatedBoundary(msgAxis) && getCommIndex(msgAxis) == 0)
      starts[msgAxis] += 1;
    if (proc >= 0)
    {

#ifdef HAVE_MPI

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
      msg << endl << "P" << rank << ": axis " << msgAxis
          << ", Lower send message:" << endl << "  ndims    = " << ndims
          << endl << "  sizes    = " << sizes << endl << "  subsizes = "
          << subsizes << endl << "  starts   = " << starts << endl
          << "  order    = " << order << endl;
#endif
      Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(new MPI_Datatype);
      MPI_Type_create_subarray(ndims,
                               &sizes[0],
                               &subsizes[0],
                               &starts[0],
                               order,
                               datatype,
                               commPad.get());
      MPI_Type_commit(commPad.get());
      messageInfo.datatype = commPad;
#else
      messageInfo.dataview = _mdArrayView;
      for (int axis = 0; axis < numDims(); ++axis)
      {
        Slice slice(starts[axis], starts[axis] + subsizes[axis]);
        messageInfo.dataview = MDArrayView< Scalar >(messageInfo.dataview,
                                                     axis,
                                                     slice);
      }
#endif

    }
    _sendMessages[msgAxis][0] = messageInfo;

    // Set the upper receive and send messages

    proc = getUpperNeighbor(msgAxis);

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
    msg << endl << "P" << rank << ": axis " << msgAxis << ", upper neighbor = "
        << proc << endl;
#endif

    subsizes[msgAxis] = getUpperPadSize(msgAxis);
    starts[msgAxis]  = _mdArrayView.dimension(msgAxis) -
                       getUpperPadSize(msgAxis);
    if (subsizes[msgAxis] == 0) proc = -1;
    messageInfo.proc = proc;
    if (proc >= 0)
    {
      // Upper receive message

#ifdef HAVE_MPI

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
      msg << endl << "P" << rank << ": axis " << msgAxis
          << ", Upper receive message:" << endl << "  ndims    = " << ndims
          << endl << "  sizes    = " << sizes << endl << "  subsizes = "
          << subsizes << endl << "  starts   = " << starts << endl
          << "  order    = " << order << endl;
#endif
      Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(new MPI_Datatype);
      MPI_Type_create_subarray(ndims,
                               &sizes[0],
                               &subsizes[0],
                               &starts[0],
                               order,
                               datatype,
                               commPad.get());
      MPI_Type_commit(commPad.get());
      messageInfo.datatype = commPad;
#else
      messageInfo.dataview = _mdArrayView;
      for (int axis = 0; axis < numDims(); ++axis)
      {
        Slice slice(starts[axis], starts[axis] + subsizes[axis]);
        messageInfo.dataview = MDArrayView< Scalar >(messageInfo.dataview,
                                                     axis,
                                                     slice);
      }
#endif
    }
    _recvMessages[msgAxis][1] = messageInfo;

    // Upper send message

    starts[msgAxis] -= getUpperPadSize(msgAxis);
    if (isReplicatedBoundary(msgAxis) &&
        getCommIndex(msgAxis) == getCommDim(msgAxis)-1)
      starts[msgAxis] -= 1;
    if (proc >= 0)
    {

#ifdef HAVE_MPI

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
      msg << endl << "P" << rank << ": axis " << msgAxis
          << ", Upper send message:" << endl << "  ndims    = " << ndims
          << endl << "  sizes    = " << sizes << endl << "  subsizes = "
          << subsizes << endl << "  starts   = " << starts << endl
          << "  order    = " << order << endl;
#endif
      Teuchos::RCP< MPI_Datatype > commPad = Teuchos::rcp(new MPI_Datatype);
      MPI_Type_create_subarray(ndims,
                               &sizes[0],
                               &subsizes[0],
                               &starts[0],
                               order,
                               datatype,
                               commPad.get());
      MPI_Type_commit(commPad.get());
      messageInfo.datatype = commPad;
#else
      messageInfo.dataview = _mdArrayView;
      for (int axis = 0; axis < numDims(); ++axis)
      {
        Slice slice(starts[axis], starts[axis] + subsizes[axis]);
        messageInfo.dataview = MDArrayView< Scalar >(messageInfo.dataview,
                                                     axis,
                                                     slice);
      }
#endif

    }
    _sendMessages[msgAxis][1] = messageInfo;
  }

#ifdef DOMI_MDVECTOR_MESSAGE_INITIALIZE
  for (int proc = 0; proc < _teuchosComm->getSize(); ++proc)
  {
    if (proc == rank)
    {
      cout << msg.str();
      msg.flush();
    }
    _teuchosComm->barrier();
  }
#endif

}


template< class Scalar >
void
MDVector< Scalar >::
writeBinary(const std::string & filename,
            bool includeBndryPad) const
{
  FILE * datafile;
  // If we are using MPI and overwriting an existing file, and the new
  // file is shorter than the old file, it appears that the new file
  // will retain the old file's length and trailing data.  To prevent
  // this behavior, we open and close the file to give it zero length.
  int pid = _teuchosComm->getRank();
  if (pid == 0)
  {
    datafile = fopen(filename.c_str(), "w");
    fclose(datafile);
  }
  _teuchosComm->barrier();

  // Get the pointer to this MDVector's MDArray, including all padding
  const Scalar * buffer = getData(true).getRawPtr();

  // Compute either _fileInfo or _fileInfoWithBndry, whichever is
  // appropriate, and return a reference to that fileInfo object
  Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);

  // Parallel output
#ifdef HAVE_MPI

  // Since HAVE_MPI is defined, we know that _teuchosComm points to a
  // const Teuchos::MpiComm< int >.  We downcast, extract and
  // dereference so that we can get access to the MPI_Comm used to
  // construct it.
  Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
    Teuchos::rcp_dynamic_cast< const Teuchos::MpiComm< int > >(_teuchosComm);
  const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
    *(mpiComm->getRawMpiComm());

  // Compute the access mode
  int access = MPI_MODE_WRONLY | MPI_MODE_CREATE;

  // I copy the filename C string, because the c_str() method returns
  // a const char*, and the MPI_File_open() function requires
  // (incorrectly) a non-const char*.
  char * cstr = new char[filename.size()+1];
  std::strcpy(cstr, filename.c_str());

  // Use MPI I/O to write the binary file
  MPI_File   mpiFile;
  MPI_Status status;
  char       datarep[7] = "native";
  MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
  MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
                    *(fileInfo->filetype), datarep, MPI_INFO_NULL);
  MPI_File_write_all(mpiFile, (void*)buffer, 1, *(fileInfo->datatype),
                     &status);
  MPI_File_close(&mpiFile);

  // Delete the C string
  delete [] cstr;

  // Serial output
#else

  // Get the number of dimensions
  // int ndims = _mdMap->numDims();

  // Initialize the data file
  datafile = fopen(filename.c_str(), "w");

  // Obtain the data to write, including the boundary padding if requested
  MDArrayView< const Scalar > mdArrayView = getData(includeBndryPad);

  // Iterate over the data and write it to the data file
  typedef typename MDArrayView< Scalar >::const_iterator const_iterator;
  for (const_iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
  {
    fwrite((const void *) &(*it), sizeof(Scalar), 1, datafile);
  }

  // Close the data file
  fclose(datafile);

#endif

}


template< class Scalar >
void
MDVector< Scalar >::
readBinary(const std::string & filename,
           bool includeBndryPad)
{
  // Get the pointer to this MDVector's MDArray, including all padding
  const Scalar * buffer = getDataNonConst(true).getRawPtr();

  // Compute either _fileInfo or _fileInfoWithBndry, whichever is
  // appropriate, and return a reference to that fileInfo object
  Teuchos::RCP< FileInfo > & fileInfo = computeFileInfo(includeBndryPad);

  // Parallel input
#ifdef HAVE_MPI

  // Since HAVE_MPI is defined, we know that _teuchosComm points to a
  // const Teuchos::MpiComm< int >.  We downcast, extract and
  // dereference so that we can get access to the MPI_Comm used to
  // construct it.
  Teuchos::RCP< const Teuchos::MpiComm< int > > mpiComm =
    Teuchos::rcp_dynamic_cast< const Teuchos::MpiComm< int > >(_teuchosComm);
  const Teuchos::OpaqueWrapper< MPI_Comm > & communicator =
    *(mpiComm->getRawMpiComm());

  // Compute the access mode
  int access = MPI_MODE_RDONLY;

  // I copy the filename C string, because the c_str() method returns
  // a const char*, and the MPI_File_open() function requires
  // (incorrectly) a non-const char*.
  char * cstr = new char[filename.size()+1];
  std::strcpy(cstr, filename.c_str());

  // Use MPI I/O to read the binary file
  MPI_File   mpiFile;
  MPI_Status status;
  char       datarep[7] = "native";
  MPI_File_open(communicator(), cstr, access, MPI_INFO_NULL, &mpiFile);
  MPI_File_set_view(mpiFile, 0, mpiType< Scalar >(),
                    *(fileInfo->filetype), datarep, MPI_INFO_NULL);
  MPI_File_read_all(mpiFile, (void*)buffer, 1, *(fileInfo->datatype),
                    &status);
  MPI_File_close(&mpiFile);

  // Delete the C string
  delete [] cstr;

  // Serial output
#else

  // Get the number of dimensions
  int ndims = _mdMap->numDims();

  // Initialize the data file
  FILE * datafile;
  datafile = fopen(filename.c_str(), "r");

  // Obtain the MDArrayView to read into, including the boundary
  // padding if requested
  MDArrayView< Scalar > mdArrayView = getDataNonConst(includeBndryPad);

  // Initialize the set of indexes
  Teuchos::Array< Ordinal > index(3);
  for (int axis = 0; axis < ndims; ++axis)
    index[axis] = fileInfo->dataStart[axis];

  // Iterate over the data and read it from the data file
  typedef typename MDArrayView< Scalar >::iterator iterator;
  for (iterator it = mdArrayView.begin(); it != mdArrayView.end(); ++it)
  {
    fread(&(*it), sizeof(Scalar), 1, datafile);
  }

  // Close the data file
  fclose(datafile);

#endif

}


template< class Scalar >
Teuchos::RCP< typename MDVector< Scalar >::FileInfo > &
MDVector< Scalar >::
computeFileInfo(bool includeBndryPad) const
{
  // Work with the appropriate FileInfo object.  By using a reference
  // here, we are working directly with the member data.
  Teuchos::RCP< MDVector< Scalar >::FileInfo > & fileInfo =
    includeBndryPad ? _fileInfoWithBndry : _fileInfo;

  // If the fileInfo object already has been set, our work is done
  if (!fileInfo.is_null()) return fileInfo;

  // Initialize the new FileInfo object
  int ndims = _mdMap->numDims();
  fileInfo.reset(new MDVector< Scalar >::FileInfo);
  fileInfo->fileShape.resize(ndims);
  fileInfo->bufferShape.resize(ndims);
  fileInfo->dataShape.resize(ndims);
  fileInfo->fileStart.resize(ndims);
  fileInfo->dataStart.resize(ndims);

  // Initialize the shapes and starts.
  for (int axis = 0; axis < ndims; ++axis)
  {
    // Initialize the FileInfo arrays using includeBndryPad where
    // appropriate
    fileInfo->fileShape[axis]   = getGlobalDim(axis,includeBndryPad);
    fileInfo->bufferShape[axis] = getLocalDim(axis,true );
    fileInfo->dataShape[axis]   = getLocalDim(axis,false);
    fileInfo->fileStart[axis]   = getGlobalRankBounds(axis,includeBndryPad).start();
    fileInfo->dataStart[axis]   = getLocalBounds(axis).start();
    // Modify dataShape and dataStart if boundary padding is included
    if (includeBndryPad)
    {
      int commIndex = _mdMap->getCommIndex(axis);
      if (commIndex == 0)
      {
        int pad = getLowerBndryPad(axis);
        fileInfo->dataShape[axis] += pad;
        fileInfo->dataStart[axis] -= pad;
      }
      if (commIndex == _mdMap->getCommDim(axis)-1)
      {
        fileInfo->dataShape[axis] += getUpperBndryPad(axis);
      }
    }
  }

#ifdef DOMI_MDVECTOR_DEBUG_IO
  cout << pid << ": fileShape   = " << fileInfo->fileShape()   << endl;
  cout << pid << ": bufferShape = " << fileInfo->bufferShape() << endl;
  cout << pid << ": dataShape   = " << fileInfo->dataShape()   << endl;
  cout << pid << ": fileStart   = " << fileInfo->fileStart()   << endl;
  cout << pid << ": dataStart   = " << fileInfo->dataStart()   << endl;
#endif

#ifdef HAVE_MPI
  int mpi_order = getLayout() == C_ORDER ? MPI_ORDER_C : MPI_ORDER_FORTRAN;
  // Build the MPI_Datatype for the file
  fileInfo->filetype = Teuchos::rcp(new MPI_Datatype);
  MPI_Type_create_subarray(ndims,
                           fileInfo->fileShape.getRawPtr(),
                           fileInfo->dataShape.getRawPtr(),
                           fileInfo->fileStart.getRawPtr(),
                           mpi_order,
                           mpiType< Scalar >(),
                           fileInfo->filetype.get());
  MPI_Type_commit(fileInfo->filetype.get());

  // Build the MPI_Datatype for the data
  fileInfo->datatype = Teuchos::rcp(new MPI_Datatype);
  MPI_Type_create_subarray(ndims,
                           fileInfo->bufferShape.getRawPtr(),
                           fileInfo->dataShape.getRawPtr(),
                           fileInfo->dataStart.getRawPtr(),
                           mpi_order,
                           mpiType< Scalar >(),
                           fileInfo->datatype.get());
  MPI_Type_commit(fileInfo->datatype.get());
#endif  // DGM_PARALLEL

  return fileInfo;
}


}  // Namespace Domi

#endif