libtrilinos-xpetra-dev 12.12.1-5 / usr / include / trilinos / Xpetra_BlockedMultiVector.hpp

// @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)
//                    Tobias Wiesner    (tawiesn@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef XPETRA_BLOCKEDMULTIVECTOR_HPP
#define XPETRA_BLOCKEDMULTIVECTOR_HPP

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

#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Map.hpp"
#include "Xpetra_MultiVector.hpp"

#include "Xpetra_BlockedMap.hpp"

namespace Xpetra {

#ifndef DOXYGEN_SHOULD_SKIP_THIS
  // forward declaration of Vector, needed to prevent circular inclusions
  template<class S, class LO, class GO, class N> class Vector;
  template<class S, class LO, class GO, class N> class BlockedVector;
  template<class S, class LO, class GO, class N> class MapExtractor;
  template<class S, class LO, class GO, class N> class MultiVectorFactory;
#endif

  template <class Scalar = double,
            class LocalOrdinal = Map<>::local_ordinal_type,
            class GlobalOrdinal = typename Map<LocalOrdinal>::global_ordinal_type,
            class Node = typename Map<LocalOrdinal, GlobalOrdinal>::node_type>
  class BlockedMultiVector
    : public MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >
  {
  public:
    typedef Scalar scalar_type;
    typedef LocalOrdinal local_ordinal_type;
    typedef GlobalOrdinal global_ordinal_type;
    typedef Node node_type;

  private:
#undef XPETRA_BLOCKEDMULTIVECTOR_SHORT
#include "Xpetra_UseShortNames.hpp"

  public:
    //! @name Constructor/Destructor Methods
    //@{

    //! Constructor

    /*!
     * Const version of constructor which accepts a const version
     * of a blocked map
     *
     * \param map BlockedMap defining the block structure of the multi vector
     * \param NumVectors Number of vector columns in multi vector
     * \param zeroOut If true initialize multivector with zeros
     */
    BlockedMultiVector(const Teuchos::RCP< const BlockedMap > &map,
        size_t NumVectors,
        bool zeroOut=true) :
          map_(map) {

      numVectors_ = NumVectors;

      vv_.reserve(map->getNumMaps());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < map->getNumMaps(); ++r)
        vv_.push_back(Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(map->getMap(r, map_->getThyraMode()), NumVectors, zeroOut));
    };

    /*!
     * Const version of constructor which accepts a const version
     * of the multi-vector
     *
     * \note If you change the information in input vector v the data in the
     *       blocked multi-vector are not affected (and vice versa). Consider
     *       the blocked multivector to be a copy of the input multivector (not a view)
     *
     * \param bmap BlockedMap object containing information about the block splitting
     * \param v MultiVector that is to be splitted into a blocked multi vector
     */
    BlockedMultiVector(Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > bmap,
                     Teuchos::RCP<const MultiVector> v) {
      XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " << bmap->getMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements() << ".");
      XPETRA_TEST_FOR_EXCEPTION(bmap->getMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has " << bmap->getMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements() << ".");
      //TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap. The BlockedMap has " << bmap->getFullMap()->getNodeNumElements() << " local elements. The vector has " << v->getMap()->getNodeNumElements() << ".");
      //TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap. The BlockedMap has " << bmap->getFullMap()->getGlobalNumElements() << " local elements. The vector has " << v->getMap()->getGlobalNumElements() << ".");

      numVectors_ = v->getNumVectors();

      map_ = bmap;

      // Extract vector
      vv_.reserve(bmap->getNumMaps());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < bmap->getNumMaps(); ++r)
        vv_.push_back(this->ExtractVector(v, r, bmap->getThyraMode()));
    }

    /*!
     * NonConst version of constructor which accepts a const version
     * of the multi-vector
     *
     * \note If you change the information in input vector v the data in the
     *       blocked multi-vector are not affected (and vice versa). Consider
     *       the blocked multivector to be a copy of the input multivector (not a view)
     *
     * \param bmap BlockedMap object containing information about the block splitting
     * \param v MultiVector that is to be splitted into a blocked multi vector
     */
    BlockedMultiVector(Teuchos::RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > bmap,
                     Teuchos::RCP<MultiVector> v) {
      TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getNodeNumElements() != v->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of local elements of MultiVector and BlockedMap");
      TEUCHOS_TEST_FOR_EXCEPTION(bmap->getFullMap()->getGlobalNumElements() != v->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector: inconsistent number of global elements of MultiVector and BlockedMap");

      numVectors_ = v->getNumVectors();

      map_ = bmap;

      // Extract vector
      vv_.reserve(bmap->getNumMaps());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < bmap->getNumMaps(); ++r)
        vv_.push_back(this->ExtractVector(v, r, bmap->getThyraMode()));
    }

    /*!
     * Const version of constructor which accepts a const version
     * of the multi-vector
     *
     * \note If you change the information in input vector v the data in the
     *       blocked multi-vector are not affected (and vice versa). Consider
     *       the blocked multivector to be a copy of the input multivector (not a view)
     *
     * \param mapExtractor MapExtractor object containing information about the block splitting
     * \param v MultiVector that is to be splitted into a blocked multi vector
     */
    BlockedMultiVector(Teuchos::RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mapExtractor,
                     Teuchos::RCP<const MultiVector> v) {
      numVectors_ = v->getNumVectors();

      // create blocked map out of MapExtractor
      std::vector<RCP<const Map> > maps;
      maps.reserve(mapExtractor->NumMaps());
      for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
        maps.push_back(mapExtractor->getMap(r,mapExtractor->getThyraMode()));
      map_ = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(mapExtractor->getFullMap(),maps,mapExtractor->getThyraMode()));

      // Extract vector
      vv_.reserve(mapExtractor->NumMaps());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
        vv_.push_back(this->ExtractVector(v, r, mapExtractor->getThyraMode()));
    }

    /*!
     * Non-const version of constructor which accepts a non-const version
     * of the multi-vector
     *
     * \note If you change the information in input vector v the data in the
     *       blocked multi-vector are not affected (and vice versa). Consider
     *       the blocked multivector to be a copy of the input multivector (not a view)
     *
     * \param mapExtractor MapExtractor object containing information about the block splitting
     * \param v MultiVector that is to be splitted into a blocked multi vector
     */
    BlockedMultiVector(Teuchos::RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mapExtractor,
                     Teuchos::RCP<MultiVector> v) {
      numVectors_ = v->getNumVectors();

      // create blocked map out of MapExtractor
      std::vector<RCP<const Map> > maps;
      maps.reserve(mapExtractor->NumMaps());
      for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
        maps.push_back(mapExtractor->getMap(r,mapExtractor->getThyraMode()));
      map_ = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(mapExtractor->getFullMap(),maps,mapExtractor->getThyraMode()));

      vv_.reserve(mapExtractor->NumMaps());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < mapExtractor->NumMaps(); ++r)
        vv_.push_back(this->ExtractVector(v, r, mapExtractor->getThyraMode()));
    }

    //! Destructor.
    virtual ~BlockedMultiVector() {
      for (size_t r = 0; r < vv_.size(); ++r)
        vv_[r] = Teuchos::null;
      map_ = Teuchos::null;
      numVectors_ = 0;
    }

    /// \brief Assignment operator: Does a deep copy.
    ///
    /// The assignment operator does a deep copy, just like
    /// subclasses' copy constructors.
    ///
    /// \note This currently only works if both <tt>*this</tt> and the
    ///   input argument are instances of the same subclass.
    BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&
    operator= (const MultiVector& rhs) {
      assign (rhs); // dispatch to protected virtual method
      return *this;
    }

    //@}
    //! @name Post-construction modification routines
    //@{

    //! Replace value, using global (row) index.
    virtual void replaceGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceGlobalValue: Not (yet) supported by BlockedMultiVector.");
    }

    //! Add value to existing value, using global (row) index.
    virtual void sumIntoGlobalValue(GlobalOrdinal globalRow, size_t vectorIndex, const Scalar &value) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoGlobalValue: Not (yet) supported by BlockedMultiVector.");
    }

    //! Replace value, using local (row) index.
    virtual void replaceLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceLocalValue: Not supported by BlockedMultiVector.");
    }

    //! Add value to existing value, using local (row) index.
    virtual void sumIntoLocalValue(LocalOrdinal myRow, size_t vectorIndex, const Scalar &value) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::sumIntoLocalValue:Not (yet) supported by BlockedMultiVector.");
    }

    //! Set all values in the multivector with the given value.
    virtual void putScalar(const Scalar &value) {
      XPETRA_MONITOR("BlockedMultiVector::putScalar");
      for(size_t r = 0; r < map_->getNumMaps(); r++) {
        getMultiVector(r)->putScalar(value);
      }
    }

    //@}

    //! @name Data Copy and View get methods
    //@{

    //! Return a Vector which is a const view of column j.
    virtual Teuchos::RCP< const Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > getVector(size_t j) const {
      RCP<Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ret =
          Teuchos::rcp(new Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(this->getBlockedMap(),false));

      for(size_t r=0; r<getBlockedMap()->getNumMaps(); r++) {
        RCP<const Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > subvec =
            this->getMultiVector(r,this->getBlockedMap()->getThyraMode())->getVector(j);
        ret->setMultiVector(r,subvec,this->getBlockedMap()->getThyraMode());
      }
      return ret;
    }

    //! Return a Vector which is a nonconst view of column j.
    virtual Teuchos::RCP< Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > getVectorNonConst(size_t j) {
      RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ret =
          Teuchos::rcp_const_cast<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(getVector(j));
      return ret;
    }

    //! Const view of the local values in a particular vector of this multivector.
    virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const {
      if(map_->getNumMaps() == 1) {
        return vv_[0]->getData(j);
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getData: Not (yet) supported by BlockedMultiVector.");
      return Teuchos::null;
    }

    //! View of the local values in a particular vector of this multivector.
    virtual Teuchos::ArrayRCP< Scalar > getDataNonConst(size_t j) {
      if(map_->getNumMaps() == 1) {
        return vv_[0]->getDataNonConst(j);
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getDataNonConst: Not (yet) supported by BlockedMultiVector.");
      return Teuchos::null;
    }

    //@}

    //! @name Mathematical methods
    //@{

    //! Compute dot product of each corresponding pair of vectors, dots[i] = this[i].dot(A[i]).
    virtual void dot(const MultiVector&A, const Teuchos::ArrayView< Scalar > &dots) const {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::dot: Not (yet) supported by BlockedMultiVector.");
    }

    //! Put element-wise absolute values of input Multi-vector in target: A = abs(this).
    virtual void abs(const MultiVector&A) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::abs: Not (yet) supported by BlockedMultiVector.");
    }

    //! Put element-wise reciprocal values of input Multi-vector in target, this(i,j) = 1/A(i,j).
    virtual void reciprocal(const MultiVector&A) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::reciprocal: Not (yet) supported by BlockedMultiVector.");
    }

    //! Scale the current values of a multi-vector, this = alpha*this.
    virtual void scale(const Scalar &alpha) {
      XPETRA_MONITOR("BlockedMultiVector::scale (Scalar)");
      for (size_t r = 0; r < map_->getNumMaps(); ++r) {
        if(getMultiVector(r)!=Teuchos::null) {
          getMultiVector(r)->scale(alpha);
        }
      }
    }

    //! Scale the current values of a multi-vector, this[j] = alpha[j]*this[j].
    virtual void scale (Teuchos::ArrayView< const Scalar > alpha) {
      XPETRA_MONITOR("BlockedMultiVector::scale (Array)");
      for (size_t r = 0; r < map_->getNumMaps(); ++r) {
        if(getMultiVector(r)!=Teuchos::null) {
          getMultiVector(r)->scale(alpha);
        }
      }
    }

    //! Update multi-vector values with scaled values of A, this = beta*this + alpha*A.
    virtual void update(const Scalar &alpha, const MultiVector&A, const Scalar &beta) {
      XPETRA_MONITOR("BlockedMultiVector::update");
      Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
      Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
      TEUCHOS_TEST_FOR_EXCEPTION(numVectors_ != rcpA->getNumVectors(),Xpetra::Exceptions::RuntimeError,"BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector).");
      if(bA != Teuchos::null) {
        // A is a BlockedMultiVector (and compatible with this)
        // Call update recursively on all sub vectors
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different thyra mode).");
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors).");
        for(size_t r = 0; r < map_->getNumMaps(); r++) {

          XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != bA->getMultiVector(r)->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: in subvector " << r << ": Cannot add a vector of (local) length " << bA->getMultiVector(r)->getMap()->getNodeNumElements() << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
          XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() != bA->getMultiVector(r)->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: in subvector " << r << ": Cannot add a vector of length " << bA->getMultiVector(r)->getMap()->getGlobalNumElements() << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
          // TAW: 12/6 We basically want to do something like:
          //   getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta);
          // But if the left hand side is a standard MultiVector and bA->getMultiVector(r) is
          // a blocked MultiVector this will not work, as the update implementation of the standard
          // multivector cannot deal with blocked multivectors.
          // The only use case where this could happen is, if the rhs vector is a single block multivector
          Teuchos::RCP<MultiVector> lmv = getMultiVector(r);
          Teuchos::RCP<MultiVector> rmv = bA->getMultiVector(r);

          // check whether lmv/rmv are blocked or not
          Teuchos::RCP<BlockedMultiVector> blmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(lmv);
          Teuchos::RCP<BlockedMultiVector> brmv = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rmv);

          if(blmv.is_null() == true && brmv.is_null() == false) {
            // special case: lmv is standard MultiVector but rmv is BlockedMultiVector
            TEUCHOS_TEST_FOR_EXCEPTION(brmv->getBlockedMap()->getNumMaps() > 1, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: Standard MultiVector object does not accept BlockedMultVector object as parameter in update call.");
            lmv->update(alpha, *(brmv->getMultiVector(0)), beta);
          } else
            lmv->update(alpha, *rmv, beta);
        }
      } else {
        // A is a MultiVector
        // If this is a BlockedMultiVector with only one sub-vector of same length we can just update
        // Otherwise, A is not compatible with this as BlockedMultiVector and we have to extract the vector data from A

        if(getBlockedMap()->getNumMaps() == 1) {
          // Actually call "update" on the underlying MultiVector (Epetra or Tpetra).
          // The maps have to be compatible.
          XPETRA_TEST_FOR_EXCEPTION(getMultiVector(0)->getMap()->isSameAs(*(rcpA->getMap()))==false, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor).");
          getMultiVector(0)->update(alpha,*rcpA,beta);
        } else {
          // general case: A has to be splitted and subparts have to be extracted and stored in this BlockedMultiVector
          //XPETRA_TEST_FOR_EXCEPTION(map_->getFullMap()->isSameAs(*(rcpA->getMap()))==false, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (maps of full vector do not match with map in MapExtractor). - Note: This test might be too strict and can probably be relaxed!");
          for(size_t r = 0; r < map_->getNumMaps(); r++) {
            // Call "update" on the subvector. Note, that getMultiVector(r) could return another BlockedMultiVector.
            // That is, in Thyra mode the maps could differ (local Xpetra versus Thyra style gids)
            Teuchos::RCP<const MultiVector> part = this->ExtractVector(rcpA, r, map_->getThyraMode());
            XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getNodeNumElements() != part->getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: in subvector " << r << ": Cannot add a vector of (local) length " << part->getMap()->getNodeNumElements() << " to the existing vector with " << getMultiVector(r)->getMap()->getNodeNumElements() << " entries.");
            XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->getGlobalNumElements() !=part->getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: in subvector " << r << ": Cannot add a vector of length " << part->getMap()->getGlobalNumElements() << " to the existing vector with " << getMultiVector(r)->getMap()->getGlobalNumElements() << " entries.");
            getMultiVector(r)->update(alpha, *part, beta);
          }
        }
      }
    }

    //! Update multi-vector with scaled values of A and B, this = gamma*this + alpha*A + beta*B.
    virtual void update(const Scalar &alpha, const MultiVector&A, const Scalar &beta, const MultiVector&B, const Scalar &gamma) {
      XPETRA_MONITOR("BlockedMultiVector::update2");
      Teuchos::RCP<const MultiVector> rcpA = Teuchos::rcpFromRef(A);
      Teuchos::RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
      Teuchos::RCP<const BlockedMultiVector> bA = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpA);
      Teuchos::RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
      if(bA != Teuchos::null && bB != Teuchos::null) {
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bA->getBlockedMap()->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector A).");
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bA->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector A).");
        TEUCHOS_TEST_FOR_EXCEPTION(numVectors_ != bA->getNumVectors(),Xpetra::Exceptions::RuntimeError,"BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector A).");
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() != bB->getBlockedMap()->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different thyra mode in vector B).");
        TEUCHOS_TEST_FOR_EXCEPTION(map_->getNumMaps() != bB->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different number of partial vectors in vector B).");
        TEUCHOS_TEST_FOR_EXCEPTION(numVectors_ != bB->getNumVectors(),Xpetra::Exceptions::RuntimeError,"BlockedMultiVector::update: update with incompatible vector (different number of vectors in multivector in vector B).");

        for(size_t r = 0; r < map_->getNumMaps(); r++) {
          XPETRA_TEST_FOR_EXCEPTION(getMultiVector(r)->getMap()->isSameAs(*(bA->getMultiVector(r)->getMap()))==false, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::update: update with incompatible vector (different maps in partial vector " << r << ").");
          getMultiVector(r)->update(alpha, *(bA->getMultiVector(r)), beta, *(bB->getMultiVector(r)), gamma);
        }
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::update: only supports update with other BlockedMultiVector.");
    }

    //! Compute 1-norm of each vector in multi-vector.
    virtual void norm1(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const {
      XPETRA_MONITOR("BlockedMultiVector::norm1");
      typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
      Array<Magnitude> temp_norms(getNumVectors());
      std::fill(norms.begin(),norms.end(),ScalarTraits<Magnitude>::zero());
      std::fill(temp_norms.begin(),temp_norms.end(),ScalarTraits<Magnitude>::zero());
      for (size_t r = 0; r < map_->getNumMaps(); ++r) {
        if(getMultiVector(r)!=Teuchos::null) {
          getMultiVector(r)->norm1(temp_norms);
          for (size_t c = 0; c < getNumVectors(); ++c)
            norms[c] += temp_norms[c];
        }
      }
    }

    //!
    virtual void norm2(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const {
      XPETRA_MONITOR("BlockedMultiVector::norm2");
      typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
      Array<Magnitude> results(getNumVectors());
      Array<Magnitude> temp_norms(getNumVectors());
      std::fill(norms.begin(),norms.end(),ScalarTraits<Magnitude>::zero());
      std::fill(results.begin(),results.end(),ScalarTraits<Magnitude>::zero());
      std::fill(temp_norms.begin(),temp_norms.end(),ScalarTraits<Magnitude>::zero());
      for (size_t r = 0; r < map_->getNumMaps(); ++r) {
        if(getMultiVector(r)!=Teuchos::null) {
          getMultiVector(r)->norm2(temp_norms);
          for (size_t c = 0; c < getNumVectors(); ++c)
            results[c] += temp_norms[c] * temp_norms[c];
        }
      }
      for (size_t c = 0; c < getNumVectors(); ++c)
        norms[c] = Teuchos::ScalarTraits<Magnitude >::squareroot(results[c]);
    }

    //! Compute Inf-norm of each vector in multi-vector.
    virtual void normInf(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const {
      XPETRA_MONITOR("BlockedMultiVector::normInf");
      typedef typename ScalarTraits<Scalar>::magnitudeType Magnitude;
      Array<Magnitude> temp_norms(getNumVectors());
      std::fill(norms.begin(),norms.end(),ScalarTraits<Magnitude>::zero());
      std::fill(temp_norms.begin(),temp_norms.end(),ScalarTraits<Magnitude>::zero());
      for (size_t r = 0; r < map_->getNumMaps(); ++r) {
        if(getMultiVector(r)!=Teuchos::null) {
          getMultiVector(r)->normInf(temp_norms);
          for (size_t c = 0; c < getNumVectors(); ++c)
            norms[c] = std::max(norms[c],temp_norms[c]);
        }
      }
    }

    //! Compute mean (average) value of each vector in multi-vector. The outcome of this routine is undefined for non-floating point scalar types (e.g., int).
    virtual void meanValue(const Teuchos::ArrayView< Scalar > &means) const {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::meanValue: Not (yet) supported by BlockedMultiVector.");
    }

    //! Matrix-matrix multiplication: this = beta*this + alpha*op(A)*op(B).
    virtual void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const Scalar &alpha, const MultiVector&A, const MultiVector&B, const Scalar &beta) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::multiply: Not (yet) supported by BlockedMultiVector.");
    }

    //! Element-wise multiply of a Vector A with a MultiVector B.
    virtual void elementWiseMultiply(Scalar scalarAB, const Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&A, const MultiVector&B, Scalar scalarThis) {
      XPETRA_MONITOR("BlockedMultiVector::elementWiseMultiply");
      XPETRA_TEST_FOR_EXCEPTION(B.getMap()->isSameAs(*(this->getMap()))==false, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: B must have same blocked map than this.");
      //XPETRA_TEST_FOR_EXCEPTION(A.getMap()->isSameAs(*(this->getMap()))==false, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: A must have same blocked map than this.");
      TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getNodeNumElements() != B.getMap()->getNodeNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getNodeNumElements() << " elements, B has " << B.getMap()->getNodeNumElements() << ".");
      TEUCHOS_TEST_FOR_EXCEPTION(A.getMap()->getGlobalNumElements() != B.getMap()->getGlobalNumElements(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: A has " << A.getMap()->getGlobalNumElements() << " elements, B has " << B.getMap()->getGlobalNumElements() << ".");

      RCP<const BlockedMap> bmap = getBlockedMap();
      RCP<const Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rcpA = Teuchos::rcpFromRef(A);
      RCP<const Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
          Teuchos::rcp_dynamic_cast<const Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);

      RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
      RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
      TEUCHOS_TEST_FOR_EXCEPTION(bB.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::elementWiseMultipy: B must be a BlockedMultiVector.");

      //RCP<Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > me = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(bmap));

      for(size_t m = 0; m < bmap->getNumMaps(); m++) {
        RCP<const Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > partA =
            bA->getMultiVector(m, bmap->getThyraMode())->getVector(0);
        RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > partB =
            bB->getMultiVector(m, bmap->getThyraMode());
        RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > thisPart =
            this->getMultiVector(m, bmap->getThyraMode());

        thisPart->elementWiseMultiply(scalarAB,*partA,*partB,scalarThis);
      }
    }

    //@}

    //! @name Attribute access functions
    //@{

    //! Number of columns in the multivector.
    virtual size_t getNumVectors() const {
      return numVectors_;
    }

    //! Local number of rows on the calling process.
    virtual size_t getLocalLength() const {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::getLocalLength: routine not implemented. It has no value as one must iterate on the partial vectors.");
      return 0;
    }

    //! Global number of rows in the multivector.
    virtual global_size_t getGlobalLength() const {
      XPETRA_MONITOR("BlockedMultiVector::getGlobalLength()");
      return map_->getFullMap()->getGlobalNumElements();
    }

    //@}

    //! @name Overridden from Teuchos::Describable
    //@{

    //! A simple one-line description of this object.
    virtual std::string description() const {
      return std::string("BlockedMultiVector");
    }

    //! Print the object with the given verbosity level to a FancyOStream.
    virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const {
      out << description() << std::endl;
      for(size_t r = 0; r < map_->getNumMaps(); r++)
        getMultiVector(r)->describe(out, verbLevel);
    }

    virtual void replaceMap(const RCP<const Map>& map) {
      XPETRA_MONITOR("BlockedMultiVector::replaceMap");
      RCP<const BlockedMap> bmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map);
      if (bmap.is_null() == true) {
        // if this has more than 1 sub blocks but "map" is not a blocked map, they are very likely not compatible
        if (this->getBlockedMap()->getNumMaps() > 1) {
          throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::replaceMap: map is not of type BlockedMap. General implementation not available, yet.");
          return;
        }
        // special case: this is a blocked map with only one block
        // TODO add more debug check (especially for Thyra mode)
        std::vector<Teuchos::RCP<const Map> > subMaps (1, map);
        map_ = Teuchos::rcp(new BlockedMap(map, subMaps,this->getBlockedMap()->getThyraMode()));
        this->getMultiVector(0)->replaceMap(map);
        return;
      }
      RCP<const BlockedMap> mybmap = Teuchos::rcp_dynamic_cast<const BlockedMap>(map_);

      XPETRA_TEST_FOR_EXCEPTION(mybmap->getThyraMode() != bmap->getThyraMode(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::replaceMap: inconsistent Thyra mode");
      map_ = bmap;
      for(size_t r = 0; r < map_->getNumMaps(); r++)
        getMultiVector(r)->replaceMap(bmap->getMap(r,map_->getThyraMode()));
    }

    //! Import.
    virtual void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Import& importer, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
    }

    //! Export.
    virtual void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Import& importer, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
    }

    //! Import (using an Exporter).
    virtual void doImport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &source, const Export& exporter, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doImport: Not supported by BlockedMultiVector.");
    }

    //! Export (using an Importer).
    virtual void doExport(const DistObject<Scalar, LocalOrdinal, GlobalOrdinal, Node> &dest, const Export& exporter, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::doExport: Not supported by BlockedMultiVector.");
    }

    //@}

    //! @name Xpetra specific
    //@{

    //! Set seed for Random function.
    virtual void setSeed(unsigned int seed) {
      for(size_t r = 0; r < map_->getNumMaps(); ++r) {
        getMultiVector(r)->setSeed(seed);
      }
    }


    virtual void randomize(bool bUseXpetraImplementation = false) {
      for(size_t r = 0; r < map_->getNumMaps(); ++r) {
        getMultiVector(r)->randomize(bUseXpetraImplementation);
      }
    }

    //! Set multi-vector values to random numbers. XPetra implementation
    virtual void Xpetra_randomize()
    {
      Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Xpetra_randomize();
    }

#ifdef HAVE_XPETRA_KOKKOS_REFACTOR
    typedef typename Kokkos::Details::ArithTraits<Scalar>::val_type impl_scalar_type;
    typedef Kokkos::DualView<impl_scalar_type**, Kokkos::LayoutStride,
        typename node_type::execution_space,
        Kokkos::MemoryUnmanaged> dual_view_type;
    typedef typename dual_view_type::host_mirror_space host_execution_space;
    typedef typename dual_view_type::t_dev::execution_space dev_execution_space;

    /// \brief Return an unmanaged non-const view of the local data on a specific device.
    /// \tparam TargetDeviceType The Kokkos Device type whose data to return.
    ///
    /// \warning Be aware that the view on the multivector data is non-persisting, i.e.
    ///          only valid as long as the multivector does not run of scope!
    template<class TargetDeviceType>
    typename Kokkos::Impl::if_c<
      Kokkos::Impl::is_same<
        typename dev_execution_space::memory_space,
        typename TargetDeviceType::memory_space>::value,
        typename dual_view_type::t_dev_um,
        typename dual_view_type::t_host_um>::type
    getLocalView () const {
      if(Kokkos::Impl::is_same<
                   typename host_execution_space::memory_space,
                   typename TargetDeviceType::memory_space
         >::value) {
        return getHostLocalView();
      } else {
        return getDeviceLocalView();
      }
    }

    virtual typename dual_view_type::t_host_um getHostLocalView ()  const {
      typename dual_view_type::t_host_um test;
      return test;
    }
    virtual typename dual_view_type::t_dev_um  getDeviceLocalView() const {
      typename dual_view_type::t_dev_um test;
      return test;
    }

#endif

    //@}

    //! Access function for the underlying Map this DistObject was constructed with.
    Teuchos::RCP< const Map> getMap() const { XPETRA_MONITOR("BlockedMultiVector::getMap"); return map_; }

    //! Access function for the underlying Map this DistObject was constructed with.
    Teuchos::RCP< const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > getBlockedMap() const { XPETRA_MONITOR("BlockedMultiVector::getBlockedMap"); return map_; }

    /// return partial multivector associated with block row r
    Teuchos::RCP<MultiVector> getMultiVector(size_t r) const {
      XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r)");
      TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(), std::out_of_range, "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
      return vv_[r];
    }

    /// return partial multivector associated with block row r
    Teuchos::RCP<MultiVector> getMultiVector(size_t r, bool bThyraMode) const {
      XPETRA_MONITOR("BlockedMultiVector::getMultiVector(r,bThyraMode)");
      TEUCHOS_TEST_FOR_EXCEPTION(r > map_->getNumMaps(), std::out_of_range, "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::getMultiVector: inconsistent Thyra mode");
      return vv_[r];
    }

    /// set partial multivector associated with block row r
    void setMultiVector(size_t r, Teuchos::RCP<const MultiVector> v, bool bThyraMode) {
      // The map of the MultiVector should be the same as the stored submap
      // In thyra mode the vectors should live on the thyra maps
      // in xpetra mode the should live in the xpetra maps
      // that should be also ok in the nested case for thyra (if the vectors are distributed accordingly)
      XPETRA_MONITOR("BlockedMultiVector::setMultiVector");
      XPETRA_TEST_FOR_EXCEPTION(r >= map_->getNumMaps(), std::out_of_range, "Error, r = " << r << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(numVectors_ != v->getNumVectors(),Xpetra::Exceptions::RuntimeError,"The BlockedMultiVectors expects " << getNumVectors() << " vectors. The provided partial multivector has " << v->getNumVectors() << " vectors.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() != bThyraMode, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::setMultiVector: inconsistent Thyra mode");
      Teuchos::RCP<MultiVector> vv = Teuchos::rcp_const_cast<MultiVector>(v);
      TEUCHOS_TEST_FOR_EXCEPTION(vv==Teuchos::null, Xpetra::Exceptions::RuntimeError, "Partial vector must not be Teuchos::null");
      vv_[r] = vv;
    }

    /// merge BlockedMultiVector blocks to a single MultiVector
    Teuchos::RCP<MultiVector> Merge() const {
      XPETRA_MONITOR("BlockedMultiVector::Merge");
      using Teuchos::RCP;

      RCP<MultiVector> v = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(map_->getFullMap(), getNumVectors());
      for(size_t r = 0; r < map_->getNumMaps(); ++r) {
        RCP<MultiVector> vi = getMultiVector(r);
        RCP<BlockedMultiVector> bvi = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vi);
        if(bvi.is_null() == true) {
          this->InsertVector(vi,r,v,map_->getThyraMode());
        } else {
          RCP<MultiVector> mvi = bvi->Merge();
          this->InsertVector(mvi,r,v,map_->getThyraMode());
        }
      }

      // TODO plausibility checks

      return v;
    }


  protected:
    /// \brief Implementation of the assignment operator (operator=);
    ///   does a deep copy.
    ///
    /// Each subclass must implement this.  This includes
    /// Xpetra::EpetraMultiVector and Xpetra::TpetraMultiVector as well as
    /// Xpetra::BockedMultiVector
    virtual void assign (const MultiVector& rhs) {
      Teuchos::RCP<const MultiVector> rcpRhs = Teuchos::rcpFromRef(rhs);
      Teuchos::RCP<const BlockedMultiVector> bRhs = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpRhs);
      if(bRhs == Teuchos::null)
        throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: argument is not a blocked multi vector.");

      map_ = Teuchos::rcp(new BlockedMap(*(bRhs->getBlockedMap())));
      vv_ = std::vector<Teuchos::RCP<MultiVector> >(map_->getNumMaps());
      for(size_t r = 0; r < map_->getNumMaps(); ++r) {
        // extract source vector (is of type MultiVector or BlockedMultiVector)
        RCP<MultiVector> src = bRhs->getMultiVector(r,map_->getThyraMode());

        // create new (empty) multivector (is of type MultiVector or BlockedMultiVector)
        RCP<MultiVector> vv = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(map_->getMap(r,bRhs->getBlockedMap()->getThyraMode()),rcpRhs->getNumVectors(),true);

        // check type of source and target vector
        RCP<BlockedMultiVector> bsrc = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(src);
        RCP<BlockedMultiVector> bvv  = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vv);
        if(bsrc.is_null() == true && bvv.is_null() == true) {
          *vv = *src; // deep copy
        } else if (bsrc.is_null() == true && bvv.is_null() == false) {
          // special case: source vector is a merged MultiVector but target vector is blocked
          *vv = *src; // deep copy (is this a problem???)
        } else if (bsrc.is_null() == false && bvv.is_null() == true) {
          // special case: source vector is blocked but target vector is a merged MultiVector
          // This is a problem and only works if bsrc has only one block
          if(bsrc->getBlockedMap()->getNumMaps() > 1) {
            throw Xpetra::Exceptions::RuntimeError("BlockedMultiVector::assign: source vector is of type BlockedMultiVector (with more than 1 blocks) and target is a MultiVector.");
            return;
          }
          RCP<MultiVector> ssrc = bsrc->getMultiVector(0,map_->getThyraMode());
          XPETRA_TEST_FOR_EXCEPTION(ssrc.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: cannot extract vector");
          XPETRA_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<BlockedMultiVector>(ssrc) != Teuchos::null, Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: sub block must not be of type BlockedMultiVector.");
          *vv = *ssrc;
        } else {
          // this should work (no exception necessary, i guess)
          XPETRA_TEST_FOR_EXCEPTION(bsrc->getBlockedMap()->getNumMaps() != bvv->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError, "BlockedMultiVector::assign: source and target are BlockedMultiVectors with a different number of submaps.");
          *vv = *src; // deep copy
        }
        vv_[r] = vv;
      }
      numVectors_ = rcpRhs->getNumVectors();
    }

  private:
    // helper routines for interaction of MultiVector and BlockedMultiVectors

    void ExtractVector(RCP<const MultiVector>& full, size_t block, RCP<MultiVector>& partial) const { ExtractVector(*full, block, *partial); }
    void ExtractVector(RCP<      MultiVector>& full, size_t block, RCP<MultiVector>& partial) const { ExtractVector(*full, block, *partial); }

    RCP<MultiVector> ExtractVector(RCP<const MultiVector>& full, size_t block, bool bThyraMode = false) const {
      XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(), std::out_of_range, "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false) == null, Xpetra::Exceptions::RuntimeError,
            "ExtractVector: map_->getmap(" << block << ",false) is null");
      RCP<const BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(full);
      if(bfull.is_null() == true) {
        // standard case: full is not of type BlockedMultiVector
        // first extract partial vector from full vector (using xpetra style GIDs)
        const RCP<MultiVector> vv = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(map_->getMap(block,false), full->getNumVectors(), false);
        //if(bThyraMode == false) {
        //  ExtractVector(*full, block, *vv);
        //  return vv;
        //} else {
          RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
          RCP<MultiVector> rcpNonConstFull = Teuchos::rcp_const_cast<MultiVector>(full);
          rcpNonConstFull->replaceMap(map_->getImporter(block)->getSourceMap());
          ExtractVector(*rcpNonConstFull, block, *vv);
          TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true, Xpetra::Exceptions::RuntimeError,
                     "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been created using Thyra-style numbered submaps.");
          if(bThyraMode == true)
            vv->replaceMap(map_->getMap(block,true)); // switch to Thyra-style map
          rcpNonConstFull->replaceMap(oldThyMapFull);
          return vv;
        //}
      } else {
        // special case: full is of type BlockedMultiVector
        XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError,
              "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be " << bfull->getBlockedMap()->getNumMaps() << " (number of blocks in BlockedMultiVector)");
        return bfull->getMultiVector(block,bThyraMode);
      }
    }
    RCP<MultiVector> ExtractVector(RCP<      MultiVector>& full, size_t block, bool bThyraMode = false) const {
      XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(), std::out_of_range, "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false) == null, Xpetra::Exceptions::RuntimeError,
            "ExtractVector: map_->getmap(" << block << ",false) is null");
      RCP<BlockedMultiVector> bfull = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(full);
      if(bfull.is_null() == true) {
        // standard case: full is not of type BlockedMultiVector
        // first extract partial vector from full vector (using xpetra style GIDs)
        const RCP<MultiVector> vv = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(map_->getMap(block,false), full->getNumVectors(), false);
        //if(bThyraMode == false) {
        //  ExtractVector(*full, block, *vv);
        //  return vv;
        //} else {
          RCP<const Map> oldThyMapFull = full->getMap(); // temporarely store map of full
          full->replaceMap(map_->getImporter(block)->getSourceMap());
          ExtractVector(*full, block, *vv);
          TEUCHOS_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true, Xpetra::Exceptions::RuntimeError,
                     "MapExtractor::ExtractVector: ExtractVector in Thyra-style numbering only possible if MapExtractor has been created using Thyra-style numbered submaps.");
          if(bThyraMode == true)
            vv->replaceMap(map_->getMap(block,true)); // switch to Thyra-style map
          full->replaceMap(oldThyMapFull);
          return vv;
        //}
      } else {
        // special case: full is of type BlockedMultiVector
        XPETRA_TEST_FOR_EXCEPTION(map_->getNumMaps() != bfull->getBlockedMap()->getNumMaps(), Xpetra::Exceptions::RuntimeError,
              "ExtractVector: Number of blocks in map extractor is " << map_->getNumMaps() << " but should be " << bfull->getBlockedMap()->getNumMaps() << " (number of blocks in BlockedMultiVector)");
        return bfull->getMultiVector(block,bThyraMode);
      }
    }
    void ExtractVector(const MultiVector& full, size_t block, MultiVector& partial) const {
      XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(), std::out_of_range, "ExtractVector: Error, block = " << block << " is too big. The BlockedMultiVector only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError,
            "ExtractVector: maps_[" << block << "] is null");
      partial.doImport(full, *(map_->getImporter(block)), Xpetra::INSERT);
    }

    void InsertVector(const MultiVector& partial, size_t block, MultiVector& full, bool bThyraMode = false) const {
      XPETRA_TEST_FOR_EXCEPTION(block >= map_->getNumMaps(), std::out_of_range, "ExtractVector: Error, block = " << block << " is too big. The MapExtractor only contains " << map_->getNumMaps() << " partial blocks.");
      XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false) == null, Xpetra::Exceptions::RuntimeError,
            "ExtractVector: map_->getmap(" << block << ",false) is null");
      XPETRA_TEST_FOR_EXCEPTION(map_->getThyraMode() == false && bThyraMode == true, Xpetra::Exceptions::RuntimeError,
                 "MapExtractor::InsertVector: InsertVector in Thyra-style numbering only possible if MapExtractor has been created using Thyra-style numbered submaps.");
      if(bThyraMode) {
        // NOTE: the importer objects in the BlockedMap are always using Xpetra GIDs (or Thyra style Xpetra GIDs)
        // The source map corresponds to the full map (in Xpetra GIDs) starting with GIDs from zero. The GIDs are consecutive in Thyra mode
        // The target map is the partial map (in the corresponding Xpetra GIDs)

        // TODO can we skip the Export call in special cases (i.e. Src = Target map, same length, etc...)

        // store original GIDs (could be Thyra GIDs)
        RCP<const MultiVector> rcpPartial = Teuchos::rcpFromRef(partial);
        RCP<MultiVector> rcpNonConstPartial = Teuchos::rcp_const_cast<MultiVector>(rcpPartial);
        RCP<const Map> oldThyMapPartial = rcpNonConstPartial->getMap(); // temporarely store map of partial
        RCP<const Map> oldThyMapFull = full.getMap(); // temporarely store map of full

        // check whether getMap(block,false) is identical to target map of importer
        //XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block,false)->isSameAs(*(map_->getImporter(block)->getTargetMap()))==false, Xpetra::Exceptions::RuntimeError,
        //           "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of partial vector are not identical to target Map of Importer. This should not be.");

        //XPETRA_TEST_FOR_EXCEPTION(full.getMap()->isSameAs(*(map_->getImporter(block)->getSourceMap()))==false, Xpetra::Exceptions::RuntimeError,
        //           "MapExtractor::InsertVector: InsertVector in Thyra-style mode: Xpetra GIDs of full vector are not identical to source Map of Importer. This should not be.");

        rcpNonConstPartial->replaceMap(map_->getMap(block,false)); // temporarely switch to xpetra-style map
        full.replaceMap(map_->getImporter(block)->getSourceMap());     // temporarely switch to Xpetra GIDs

        // do the Export
        full.doExport(*rcpNonConstPartial, *(map_->getImporter(block)), Xpetra::INSERT);

        // switch back to original maps
        full.replaceMap(oldThyMapFull); // reset original map (Thyra GIDs)
        rcpNonConstPartial->replaceMap(oldThyMapPartial); // change map back to original map
      } else {
        // Xpetra style numbering
        full.doExport(partial, *(map_->getImporter(block)), Xpetra::INSERT);
      }
    }
    void InsertVector(RCP<const MultiVector> partial, size_t block, RCP<MultiVector> full, bool bThyraMode = false) const {
      RCP<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bfull = Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(full);
      if(bfull.is_null() == true)
        InsertVector(*partial, block, *full, bThyraMode);
      else {
        XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError,
              "InsertVector: maps_[" << block << "] is null");
        full->setMultiVector(block, partial, bThyraMode);
      }
    }
    void InsertVector(RCP<      MultiVector> partial, size_t block, RCP<MultiVector> full, bool bThyraMode = false) const {
      RCP<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bfull = Teuchos::rcp_dynamic_cast<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(full);
      if(bfull.is_null() == true)
        InsertVector(*partial, block, *full, bThyraMode);
      else {
        XPETRA_TEST_FOR_EXCEPTION(map_->getMap(block) == null, Xpetra::Exceptions::RuntimeError,
              "InsertVector: maps_[" << block << "] is null");
        bfull->setMultiVector(block, partial, bThyraMode);
      }
    }

  private:
    Teuchos::RCP<const BlockedMap>          map_;         ///< blocked map containing the sub block maps (either thyra or xpetra mode)
    std::vector<Teuchos::RCP<MultiVector> > vv_;          ///< array containing RCPs of the partial vectors
    size_t                                  numVectors_;  ///< number of vectors (columns in multi vector)
  }; // BlockedMultiVector class

} // Xpetra namespace

#define XPETRA_BLOCKEDMULTIVECTOR_SHORT
#endif // XPETRA_BLOCKEDMULTIVECTOR_HPP