Xpetra_BlockedMap_def.hpp

Go to the documentation of this file.

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

#include "Xpetra_BlockedMap_decl.hpp"

#include "Xpetra_Exceptions.hpp"
#include "Xpetra_ImportFactory.hpp"
#include "Xpetra_MapFactory_decl.hpp"



namespace Xpetra {



template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
BlockedMap()
{
    bThyraMode_ = false;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
BlockedMap(const RCP<const Map>& fullmap, const std::vector<RCP<const Map>>& maps, bool bThyraMode)
{
    bThyraMode_ = bThyraMode;

    if(bThyraMode == false)
    {
        // use Xpetra-style numbering for sub-block maps
        // That is, all sub-block maps have unique GIDs which may not be contiguous and start with GIDs different than zero.

        // plausibility check
        size_t numAllElements = 0;
        for(size_t v = 0; v < maps.size(); ++v)
        {
            numAllElements += maps[ v ]->getGlobalNumElements();
        }
        TEUCHOS_TEST_FOR_EXCEPTION(fullmap->getGlobalNumElements() != numAllElements,
                                   std::logic_error,
                                   "logic error. full map and sub maps have not same number of elements ("
                                     << fullmap->getGlobalNumElements() << " versus " << numAllElements
                                     << "). We cannot build MapExtractor with Xpetra-style numbering. Please make sure that you want "
                                        "Xpetra-style numbering instead of Thyra-style numbering.");

        fullmap_ = fullmap;
        maps_    = maps;
    }
    else
    {
        // std::cout << "Create Map Extractor in Thyra Mode!!! " << std::endl;
        // use Thyra-style numbering for sub-block maps
        // That is, all sub-block maps start with zero as GID and are contiguous

        // plausibility check
        for(size_t v = 0; v < maps.size(); ++v)
        {
            TEUCHOS_TEST_FOR_EXCEPTION(maps[ v ]->getMinAllGlobalIndex() != 0,
                                       std::logic_error,
                                       "logic error. When using Thyra-style numbering all sub-block maps must start with zero as GID. Map block "
                                         << v << " starts with GID " << maps[ v ]->getMinAllGlobalIndex());
        }

        // store submaps in Thyra-style ordering
        thyraMaps_ = maps;

        // get offsets
        std::vector<GlobalOrdinal> gidOffsets(maps.size(), 0);
        for(size_t v = 1; v < maps.size(); ++v)
        {
            gidOffsets[ v ] = maps[ v - 1 ]->getMaxAllGlobalIndex() + gidOffsets[ v - 1 ] + 1;
        }

        // build submaps
        maps_.resize(maps.size());
        std::vector<GlobalOrdinal> fullMapGids;
        const GO                   INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
        for(size_t v = 0; v < maps.size(); ++v)
        {
            size_t                     myNumElements = maps[ v ]->getNodeNumElements();
            std::vector<GlobalOrdinal> subMapGids(myNumElements, 0);
            for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(myNumElements); ++l)
            {
                GlobalOrdinal myGid = maps[ v ]->getGlobalElement(l);
                subMapGids[ l ]     = myGid + gidOffsets[ v ];
                fullMapGids.push_back(myGid + gidOffsets[ v ]);
            }
            // std::sort(subMapGids.begin(), subMapGids.end());
            // subMapGids.erase(std::unique(subMapGids.begin(), subMapGids.end()), subMapGids.end());

            Teuchos::ArrayView<GlobalOrdinal> subMapGidsView(&subMapGids[ 0 ], subMapGids.size());

            Teuchos::RCP<Map> mySubMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(
              maps[ v ]->lib(), INVALID, subMapGidsView, maps[ v ]->getIndexBase(), maps[ v ]->getComm());
            maps_[ v ] = mySubMap;
        }

        // const GO INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
        // std::sort(coarseMapGids.begin(), coarseMapGids.end());
        // coarseMapGids.erase(std::unique(coarseMapGids.begin(), coarseMapGids.end()), coarseMapGids.end());
        // Teuchos::ArrayView<GO> coarseMapGidsView(&coarseMapGids[0], coarseMapGids.size());
        // std::sort(fullMapGids.begin(), fullMapGids.end());
        // fullMapGids.erase(std::unique(fullMapGids.begin(), fullMapGids.end()), fullMapGids.end());

        Teuchos::ArrayView<GlobalOrdinal> fullMapGidsView(&fullMapGids[ 0 ], fullMapGids.size());

        fullmap_ = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(
          fullmap->lib(), INVALID, fullMapGidsView, fullmap->getIndexBase(), fullmap->getComm());

        // plausibility check
        size_t numAllElements = 0;
        for(size_t v = 0; v < maps_.size(); ++v)
        {
            numAllElements += maps_[ v ]->getGlobalNumElements();
        }
        TEUCHOS_TEST_FOR_EXCEPTION(
          fullmap_->getGlobalNumElements() != numAllElements,
          std::logic_error,
          "logic error. full map and sub maps have not same number of elements. This cannot be. Please report the bug to the Xpetra developers!");
    }

    // build importers for sub maps
    importers_.resize(maps_.size());
    for(unsigned i = 0; i < maps_.size(); ++i)
    {
        if(maps[ i ] != null)
        {
            importers_[ i ] = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(fullmap_, maps_[ i ]);
        }
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
      CheckConsistency() == false, std::logic_error, "logic error. full map and sub maps are inconsistently distributed over the processors.");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
BlockedMap(const std::vector<RCP<const Map>>& maps, const std::vector<RCP<const Map>>& thyramaps)
{
    bThyraMode_ = true;

    // plausibility check
    TEUCHOS_TEST_FOR_EXCEPTION(thyramaps.size() != maps.size(), std::logic_error, "logic error. The number of submaps must be identical!");
    for(size_t v = 0; v < thyramaps.size(); ++v)
    {
        TEUCHOS_TEST_FOR_EXCEPTION(thyramaps[ v ]->getMinAllGlobalIndex() != 0,
                                   std::logic_error,
                                   "logic error. When using Thyra-style numbering all sub-block maps must start with zero as GID.");

        XPETRA_TEST_FOR_EXCEPTION(thyramaps[ v ]->getNodeNumElements() != maps[ v ]->getNodeNumElements(),
                                  std::logic_error,
                                  "logic error. The size of the submaps must be identical (same distribution, just different GIDs)");
    }

    // store user-provided maps and thyramaps
    thyraMaps_ = thyramaps;
    maps_      = maps;
    fullmap_   = this->concatenateMaps(maps);

    // plausibility check
    size_t numAllElements = 0;
    for(size_t v = 0; v < maps_.size(); ++v)
    {
        numAllElements += maps_[ v ]->getGlobalNumElements();
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
      fullmap_->getGlobalNumElements() != numAllElements,
      std::logic_error,
      "logic error. full map and sub maps have not same number of elements. This cannot be. Please report the bug to the Xpetra developers!");

    // build importers for sub maps
    importers_.resize(maps_.size());
    for(unsigned i = 0; i < maps_.size(); ++i)
    {
        if(maps[ i ] != null)
        {
            importers_[ i ] = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(fullmap_, maps_[ i ]);
        }
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
      CheckConsistency() == false, std::logic_error, "logic error. full map and sub maps are inconsistently distributed over the processors.");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
BlockedMap(const BlockedMap& input)
{
    bThyraMode_ = input.getThyraMode();
    fullmap_    = Teuchos::null;
    maps_.resize(input.getNumMaps(), Teuchos::null);
    thyraMaps_.resize(input.getNumMaps(), Teuchos::null);
    this->assign(input);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
~BlockedMap()
{
    // make sure all RCP's are freed
    for(size_t v = 0; v < maps_.size(); ++v)
    {
        maps_[ v ] = Teuchos::null;
        if(bThyraMode_ == true)
            thyraMaps_[ v ] = Teuchos::null;
        importers_[ v ] = Teuchos::null;
    }

    fullmap_ = Teuchos::null;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
global_size_t
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getGlobalNumElements() const
{
    return fullmap_->getGlobalNumElements();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
size_t
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::getNodeNumElements() const
{
    return fullmap_->getNodeNumElements();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getIndexBase() const
{
    return fullmap_->getIndexBase();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMinLocalIndex() const
{
    return fullmap_->getMinLocalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMaxLocalIndex() const
{
    return fullmap_->getMaxLocalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMinGlobalIndex() const
{
    return fullmap_->getMinGlobalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMaxGlobalIndex() const
{
    return fullmap_->getMaxGlobalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMinAllGlobalIndex() const
{
    return fullmap_->getMinAllGlobalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMaxAllGlobalIndex() const
{
    return fullmap_->getMaxAllGlobalIndex();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getLocalElement(GlobalOrdinal globalIndex) const
{
    return fullmap_->getLocalElement(globalIndex);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
GlobalOrdinal
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getGlobalElement(LocalOrdinal localIndex) const
{
    return fullmap_->getGlobalElement(localIndex);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
LookupStatus
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getRemoteIndexList(const Teuchos::ArrayView<const GlobalOrdinal>& /* GIDList    */,
                                            const Teuchos::ArrayView<int>& /* nodeIDList */,
                                            const Teuchos::ArrayView<LocalOrdinal>& /* LIDList    */) const
{
    throw Xpetra::Exceptions::RuntimeError("BlockedMap::getRemoteIndexList: routine not implemented.");
    TEUCHOS_UNREACHABLE_RETURN(IDNotPresent);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
LookupStatus
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getRemoteIndexList(const Teuchos::ArrayView<const GlobalOrdinal>& /* GIDList */,
                   const Teuchos::ArrayView<int>& /* nodeIDList */) const
{
    throw Xpetra::Exceptions::RuntimeError("BlockedMap::getRemoteIndexList: routine not implemented.");
    TEUCHOS_UNREACHABLE_RETURN(IDNotPresent);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::ArrayView<const GlobalOrdinal>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getNodeElementList() const
{
    return fullmap_->getNodeElementList();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isNodeLocalElement(LocalOrdinal localIndex) const
{
    return fullmap_->isNodeLocalElement(localIndex);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isNodeGlobalElement(GlobalOrdinal globalIndex) const
{
    return fullmap_->isNodeGlobalElement(globalIndex);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isContiguous() const
{
    throw Xpetra::Exceptions::RuntimeError("BlockedMap::isContiguous: routine not implemented.");
    TEUCHOS_UNREACHABLE_RETURN(false);
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isDistributed() const
{
    return fullmap_->isDistributed();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isCompatible(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) const
{
    RCP<const Map>        rcpMap  = Teuchos::rcpFromRef(map);
    RCP<const BlockedMap> rcpBMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rcpMap);
    if(rcpBMap.is_null() == true)
        return false;

    for(size_t v = 0; v < maps_.size(); ++v)
    {
        bool bSame = getMap(v, false)->isCompatible(*(rcpBMap->getMap(v, false)));
        if(bSame == false)
            return false;
        if(bThyraMode_)
        {
            bSame = getMap(v, true)->isCompatible(*(rcpBMap->getMap(v, true)));
        }
    }
    return true;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
isSameAs(const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>& map) const
{
    RCP<const Map>        rcpMap  = Teuchos::rcpFromRef(map);
    RCP<const BlockedMap> rcpBMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rcpMap);
    if(rcpBMap.is_null() == true)
    {
        // If this is a blocked map with > 1 blocks but "map" is a plain map they can't be the same
        if(this->getNumMaps() > 1)
        {
            return false;
        }

        // special case: this is a single blocked map and "map" is a plain map object
        bool bSame = getMap(0, bThyraMode_)->isSameAs(*rcpMap);
        return bSame;
    }

    for(size_t v = 0; v < maps_.size(); ++v)
    {
        bool bSame = getMap(v, false)->isSameAs(*(rcpBMap->getMap(v, false)));
        if(bSame == false)
        {
            return false;
        }
        if(bThyraMode_)
        {
            bSame = getMap(v, true)->isSameAs(*(rcpBMap->getMap(v, true)));
            if(bSame == false)
            {
                return false;
            }
        }
    }
    return true;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Teuchos::Comm<int>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getComm() const
{
    return fullmap_->getComm();
}


#ifdef TPETRA_ENABLE_DEPRECATED_CODE
template<class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getNode() const
{
    return fullmap_->getNode();
}
#endif      // TPETRA_ENABLE_DEPRECATED_CODE


template<class LocalOrdinal, class GlobalOrdinal, class Node>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>&
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
operator=(const BlockedMap& rhs)
{
    assign(rhs);      // dispatch to protected virtual method
    return *this;
}

template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::getThyraMode() const
{
    return bThyraMode_;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
removeEmptyProcesses() const
{
    throw Xpetra::Exceptions::RuntimeError("BlockedMap::removeEmptyProcesses: routine not implemented.");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
replaceCommWithSubset(const Teuchos::RCP<const Teuchos::Comm<int>>& /* newComm */) const
{
    throw Xpetra::Exceptions::RuntimeError("BlockedMap::replaceCommWithSubset: routine not implemented.");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
UnderlyingLib
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::lib() const
{
    return fullmap_->lib();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMap() const
{
    return getFullMap();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
size_t
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getNumMaps() const
{
    return maps_.size();
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
const RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMap(size_t i,
       bool   bThyraMode) const
{
    XPETRA_TEST_FOR_EXCEPTION(i >= getNumMaps(),
                              Xpetra::Exceptions::RuntimeError,
                              "BlockedMap::getMap: tried to access block " << i << ", but BlockedMap has only " << getNumMaps()
                                                                           << " blocks! Block indices must be between 0 and " << getNumMaps() - 1
                                                                           << ".");
    if(bThyraMode_ == true && bThyraMode == true)
    {
        return thyraMaps_[ i ];
    }

    XPETRA_TEST_FOR_EXCEPTION(bThyraMode_ == false && bThyraMode == true,
                              Xpetra::Exceptions::RuntimeError,
                              "BlockedMap::getMap: cannot return sub map in Thyra-style numbering if BlockedMap object is not created using "
                              "Thyra-style numbered submaps.");
    return maps_[ i ];
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
const RCP<Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getImporter(size_t i) const
{
    XPETRA_TEST_FOR_EXCEPTION(i >= getNumMaps(),
                              Xpetra::Exceptions::RuntimeError,
                              "BlockedMap::getImporter: tried to access block " << i << ", but BlockedMap has only " << getNumMaps()
                                                                                << " blocks! Block indices must be between 0 and " << getNumMaps() - 1
                                                                                << ".");
    return importers_[ i ];
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getFullMap() const
{
    return fullmap_;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
size_t
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
getMapIndexForGID(GlobalOrdinal gid) const
{
    for(size_t i = 0; i < getNumMaps(); i++)
        if(getMap(i)->isNodeGlobalElement(gid) == true)
            return i;

    TEUCHOS_TEST_FOR_EXCEPTION(
      false, Xpetra::Exceptions::RuntimeError, "getMapIndexForGID: GID " << gid << " is not contained by a map in mapextractor.");
    return 0;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
std::string
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
description() const
{
    return std::string("BlockedMap");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
void
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const
{
    out << "------------- Blocked Map -----------" << std::endl;
    out << description() << std::endl;
    out << "Thyra mode: " << getThyraMode() << std::endl;
    out << "No of submaps: " << getNumMaps() << std::endl;
    Teuchos::OSTab tab(out);
    for(size_t r = 0; r < getNumMaps(); r++)
    {
        std::cout << "MAP " << r << "/" << getNumMaps() - 1 << std::endl;
        getMap(r, false)->describe(out, verbLevel);
    }
    if(getThyraMode() == true)
    {
        for(size_t r = 0; r < getNumMaps(); r++)
        {
            std::cout << "Thyra MAP " << r << "/" << getNumMaps() - 1 << std::endl;
            getMap(r, true)->describe(out, verbLevel);
        }
    }
    out << "-------------------------------------" << std::endl;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
void
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
assign(const BlockedMap& input)
{
    // TODO check implementation, simplify copy constructor
    bThyraMode_ = input.getThyraMode();

    fullmap_ = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(input.getFullMap(), 1);

    maps_.resize(input.getNumMaps(), Teuchos::null);
    if(bThyraMode_ == true)
        thyraMaps_.resize(input.getNumMaps(), Teuchos::null);
    for(size_t i = 0; i < input.getNumMaps(); ++i)
    {
        maps_[ i ] = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(input.getMap(i, false), 1);
        if(bThyraMode_ == true)
            thyraMaps_[ i ] = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(input.getMap(i, true), 1);
    }

    // plausibility check
    size_t numAllElements = 0;
    for(size_t v = 0; v < maps_.size(); ++v)
    {
        numAllElements += maps_[ v ]->getGlobalNumElements();
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
      fullmap_->getGlobalNumElements() != numAllElements,
      std::logic_error,
      "logic error. full map and sub maps have not same number of elements. This cannot be. Please report the bug to the Xpetra developers!");

    // build importers for sub maps
    importers_.resize(maps_.size());
    for(unsigned i = 0; i < maps_.size(); ++i)
        if(maps_[ i ] != null)
            importers_[ i ] = Xpetra::ImportFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(fullmap_, maps_[ i ]);
    TEUCHOS_TEST_FOR_EXCEPTION(
      CheckConsistency() == false, std::logic_error, "logic error. full map and sub maps are inconsistently distributed over the processors.");
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
concatenateMaps(const std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>>>& subMaps)
{

    // merge submaps to global map
    std::vector<GlobalOrdinal> gids;
    for(size_t tt = 0; tt < subMaps.size(); ++tt)
    {
        Teuchos::RCP<const Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> subMap = subMaps[ tt ];

#if 1      // WCMCLEN : IS THIS NECESSARY TO HANG ONTO?
        Teuchos::ArrayView<const GlobalOrdinal> subMapGids = subMap->getNodeElementList();
        gids.insert(gids.end(), subMapGids.begin(), subMapGids.end());
#else
        size_t myNumElements = subMap->getNodeNumElements();
        for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(myNumElements); ++l)
        {
            GlobalOrdinal gid = subMap->getGlobalElement(l);
            gids.push_back(gid);
        }
#endif
    }

    const GlobalOrdinal INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
    // std::sort(gids.begin(), gids.end());
    // gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
    Teuchos::ArrayView<GlobalOrdinal>                            gidsView(&gids[ 0 ], gids.size());

    Teuchos::RCP<Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node>> fullMap = Xpetra::MapFactory<LocalOrdinal, GlobalOrdinal, Node>::
        Build(subMaps[ 0 ]->lib(), INVALID, gidsView, subMaps[ 0 ]->getIndexBase(), subMaps[ 0 ]->getComm());

    return fullMap;
}


template<class LocalOrdinal, class GlobalOrdinal, class Node>
bool
BlockedMap<LocalOrdinal, GlobalOrdinal, Node>::
CheckConsistency() const
{
    const RCP<const Map> fullMap = getFullMap();

    for(size_t i = 0; i < getNumMaps(); i++)
    {
        const RCP<const Map> map = getMap(i);

        ArrayView<const GlobalOrdinal> mapGids = map->getNodeElementList();
        for(typename ArrayView<const GlobalOrdinal>::const_iterator it = mapGids.begin(); it != mapGids.end(); it++)
        {
            if(fullMap->isNodeGlobalElement(*it) == false)
            {
                return false;      // Global ID (*it) not found locally on this proc in fullMap -> error
            }
        }
    }
    return true;
}



}      // namespace Xpetra
#endif /* PACKAGES_XPETRA_SUP_BLOCKEDMAP_XPETRA_BLOCKEDMAP_DECL_HPP_ */