Xpetra_ReorderedBlockedMultiVector.hpp

Go to the documentation of this file.

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

#include <Kokkos_DefaultNode.hpp>

#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"

#include "Xpetra_MapUtils.hpp"

#include "Xpetra_BlockReorderManager.hpp"
#include "Xpetra_BlockedMap.hpp"
#include "Xpetra_BlockedMultiVector.hpp"


namespace Xpetra {

  typedef std::string viewLabel_t;

  template <class Scalar,
            class LocalOrdinal,
            class GlobalOrdinal,
            class Node>
  class ReorderedBlockedMultiVector :
    public BlockedMultiVector<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_REORDEREDBLOCKEDMULTIVECTOR_SHORT
#include "Xpetra_UseShortNames.hpp"

  public:




    ReorderedBlockedMultiVector
        (Teuchos::RCP<const BlockedMap>& rangeMap,
         Teuchos::RCP<const Xpetra::BlockReorderManager> brm,
         Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec)
  : Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rangeMap,bvec->getNumVectors(), false) {
      brm_ = brm;
      fullVec_ = bvec;
    }

  //protected:

    virtual ~ReorderedBlockedMultiVector() {
      brm_ = Teuchos::null;
      fullVec_ = Teuchos::null;
    }


  private:
    Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm) {
      RCP<const BlockedMap> bmap = fullVec_->getBlockedMap();

      // number of sub blocks
      size_t numBlocks = brm->GetNumBlocks();

      Teuchos::RCP<const Map> map = Teuchos::null;

      if(numBlocks == 0) {
        // it is a leaf node
        Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);

        // never extract Thyra style maps (since we have to merge them)
        map = bmap->getMap(Teuchos::as<size_t>(leaf->GetIndex()), false);
      } else {
        // initialize vector for sub maps
        std::vector<Teuchos::RCP<const Map> > subMaps (numBlocks, Teuchos::null);

        for(size_t i = 0; i < numBlocks; i++) {
          Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
          subMaps[i] = mergeSubBlockMaps(blkMgr);
          TEUCHOS_ASSERT(subMaps[i].is_null()==false);
        }

        map = MapUtils::concatenateMaps(subMaps);
      }
      TEUCHOS_ASSERT(map.is_null()==false);
      return map;
    }

  public:


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

    void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
      TEUCHOS_ASSERT(brm_ != Teuchos::null);
      out << description() << ": " << brm_->toString() << std::endl;
      fullVec_->describe(out,verbLevel);
    }


  private:
    Teuchos::RCP<const Xpetra::BlockReorderManager > brm_;
    Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullVec_;


};

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec, bool bThyraMode) {
  typedef Xpetra::MapUtils<LocalOrdinal,GlobalOrdinal,Node> MapUtils;

  // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
  RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > bmap = bvec->getBlockedMap();

  // number of sub blocks
  size_t numBlocks = brm->GetNumBlocks();

  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> map = Teuchos::null;

  if(numBlocks == 0) {
    // it is a leaf node
    Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);

    map = bmap->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
  } else {
    // initialize vector for sub maps
    std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);

    for(size_t i = 0; i < numBlocks; i++) {
      Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
      subMaps[i] = mergeSubBlockMaps(blkMgr,bvec,bThyraMode);
      TEUCHOS_ASSERT(subMaps[i].is_null()==false);
    }
#if 1
    // concatenate submaps
    // for Thyra mode this map isn't important
    Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = MapUtils::concatenateMaps(subMaps);

    // create new BlockedMap (either in Thyra Mode or Xpetra mode)
    map = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(fullMap, subMaps, bThyraMode));
#else
    // TAW: 11/27/16 we just concatenate the submaps to one monolithic Map object.
    // Alternatively, we could create a new BlockedMap using the concatenated map and the submaps
    // However, the block smoothers only need the concatenated map for creating MultiVectors...
    // But for the Thyra mode version concatenating would not be ok for the whole map
    map = MapUtils::concatenateMaps(subMaps);
#endif
  }
  TEUCHOS_ASSERT(map.is_null()==false);
  return map;
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocks(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {

  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
  typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> BlockedMap;
  typedef Xpetra::MapUtils<LocalOrdinal,GlobalOrdinal,Node> MapUtils;
  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVector;
  typedef Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedMultiVector;
  typedef Xpetra::ReorderedBlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> ReorderedBlockedMultiVector;

  // number of sub blocks
  size_t rowSz = rowMgr->GetNumBlocks();

  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;

  if(rowSz == 0) {

    // it is a leaf node
    Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);

    // extract leaf node
    Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),false);

    TEUCHOS_ASSERT(vec != Teuchos::null);

    // check, whether leaf node is of type Xpetra::CrsMatrixWrap
    Teuchos::RCP<BlockedMultiVector> subBVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
    if(subBVec == Teuchos::null) {

      // DEBUG
      /*{
        std::cout << "MultiVector:" << std::endl;
        Teuchos::ArrayRCP<const Scalar> vData = vec->getData(0);
        for(size_t j=0; j< vec->getMap()->getNodeNumElements(); j++) {
            std::cout << j << ": " << vec->getMap()->getGlobalElement(j) << ": " << vData[j] << std::endl;
        }
      }*/
      // END DEBUG

      // If the leaf node is of type Xpetra::MultiVector. Wrap it into a ReorderedBlockMultiVector
      // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
      RCP<const BlockedMap> fullBlockedMap = bvec->getBlockedMap();
      Teuchos::RCP<const Map> submap = fullBlockedMap->getMap(rowleaf->GetIndex(),false);
      std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
      Teuchos::RCP<const BlockedMap> bbmap = Teuchos::rcp(new BlockedMap(submap, rowSubMaps, false));

      rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(bbmap, rowMgr, bvec));
      rbvec->setMultiVector(0,Teuchos::rcp_const_cast<MultiVector>(vec),false);

    } else {
      // If leaf node is already wrapped into a blocked matrix do not wrap it again.
      rbvec = subBVec;
      TEUCHOS_ASSERT(rbvec != Teuchos::null);
    }
  } else {
    // create the map extractors
    // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
    Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::null;
    std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
    for(size_t i = 0; i < rowSz; i++) {
      Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
      rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false /*xpetra*/);
      TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
    }
    Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
    rgBlockedMap = Teuchos::rcp(new BlockedMap(rgMergedSubMaps, rowSubMaps, false));
    rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));

    for(size_t i = 0; i < rowSz; i++) {
      Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
      Teuchos::RCP<const MultiVector> subvec = mergeSubBlocks(rowSubMgr, bvec);
      rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),false);
    }

  }
  return rbvec;
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocksThyra(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {
  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
  typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> BlockedMap;
  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVector;
  typedef Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedMultiVector;
  typedef Xpetra::ReorderedBlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> ReorderedBlockedMultiVector;

  TEUCHOS_ASSERT(bvec->getBlockedMap()->getThyraMode() == true);

  // number of sub blocks
  size_t rowSz = rowMgr->GetNumBlocks();

  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;

  if(rowSz == 0) {
    // it is a leaf node
    Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);

    // this MultiVector uses Thyra style GIDs as global row indices
    Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),true);

    TEUCHOS_ASSERT(vec.is_null() == false);

    // check, whether leaf node is of type Xpetra::CrsMatrixWrap
    Teuchos::RCP<BlockedMultiVector> bbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
    if(bbvec == Teuchos::null) {
      // build blocked map
      RCP<const BlockedMap> fullBlockedRangeMap = bvec->getBlockedMap();
      // extract Xpetra and Thyra based GIDs
      Teuchos::RCP<const Map> xpsubmap  = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),false);
      Teuchos::RCP<const Map> thysubmap = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),true);
      std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
      std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
      // use expert constructor
      Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));

      // build reordered blocked multi vector
      rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr, bvec));
      rbvec->setMultiVector(0,vec,true);
    } else {
      // If leaf node is already wrapped into a blocked matrix do not wrap it again.
      rbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
    }
  } else {
    // create the blocked map
    // we cannot create block multivector in thyra mode since merged maps might not start with 0 GID

    std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
    std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
    for(size_t i = 0; i < rowSz; i++) {
      Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
      // extract Xpetra and Thyra based merged GIDs
      rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false);
      rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,true);
      TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
      TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
    }
    // use expert constructor
    Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));

    rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));

    for(size_t i = 0; i < rowSz; i++) {
      Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
      Teuchos::RCP<const MultiVector> subvec = mergeSubBlocksThyra(rowSubMgr, bvec);
      rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),true);
    }

  }
  return rbvec;
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > buildReorderedBlockedMultiVector(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {

  Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbvec = Teuchos::null;
  if(bvec->getBlockedMap()->getThyraMode() == false) {
    rbvec = mergeSubBlocks(brm,bvec);
  } else {
    rbvec = mergeSubBlocksThyra(brm,bvec);
  }
  TEUCHOS_ASSERT(rbvec.is_null() == false);
  return rbvec;
}

template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > buildReorderedBlockedMultiVector(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {

  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVector;
  typedef Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedMultiVector;
  Teuchos::RCP<const MultiVector> rbvec = Teuchos::null;
  if(bvec->getBlockedMap()->getThyraMode() == false) {
    rbvec = mergeSubBlocks(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
  } else {
    rbvec = mergeSubBlocksThyra(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
  }
  TEUCHOS_ASSERT(rbvec.is_null() == false);
  return Teuchos::rcp_const_cast<MultiVector>(rbvec);
}


} //namespace Xpetra

#define XPETRA_REORDEREDBLOCKEDMULTIVECTOR_SHORT
#endif /* XPETRA_REORDEREDBLOCKEDMULTIVECTOR_HPP */