Xpetra_StridedMap.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
//                    Jeremie Gaidamour (jngaida@sandia.gov)
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

// WARNING: This code is experimental. Backwards compatibility should not be expected.

#ifndef XPETRA_STRIDEDMAP_HPP
#define XPETRA_STRIDEDMAP_HPP

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

#include <Kokkos_DefaultNode.hpp>
#include <Teuchos_Describable.hpp>
#include "Xpetra_Map.hpp"
#include "Xpetra_ConfigDefs.hpp"

#include "Xpetra_Exceptions.hpp"

// MPI helper
#define sumAll(rcpComm, in, out)                                        \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out));
#define minAll(rcpComm, in, out)                                        \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out));
#define maxAll(rcpComm, in, out)                                        \
  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out));

namespace Xpetra {

  template <class LocalOrdinal, class GlobalOrdinal = LocalOrdinal, class Node = Kokkos::DefaultNode::DefaultNodeType>
  class StridedMap
    : public virtual Map<LocalOrdinal, GlobalOrdinal, Node>
  {

  public:



    StridedMap(global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId, GlobalOrdinal offset = 0)
    : stridingInfo_(stridingInfo), stridedBlockId_(stridedBlockId), offset_(offset)
    { 
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() == 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: stridingInfo.size() = 0?");
      //TEUCHOS_TEST_FOR_EXCEPTION(numGlobalElements % getFixedBlockSize() != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: getFixedBlockSize is not an integer multiple of numGlobalElements.");

    }

    StridedMap(global_size_t numGlobalElements, size_t numLocalElements, GlobalOrdinal indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId, GlobalOrdinal offset = 0)
    : stridingInfo_(stridingInfo), stridedBlockId_(stridedBlockId), offset_(offset)
    { 
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() == 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: stridingInfo.size() = 0?");
      //TEUCHOS_TEST_FOR_EXCEPTION(numGlobalElements % getFixedBlockSize() != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: getFixedBlockSize is not an integer multiple of numGlobalElements.");
      //TEUCHOS_TEST_FOR_EXCEPTION(numLocalElements % getFixedBlockSize() != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: getFixedBlockSize is not an integer multiple of numLocalElements.");
    
    }

    StridedMap(std::vector<size_t>& stridingInfo, LocalOrdinal stridedBlockId, GlobalOrdinal offset = 0)
    : stridingInfo_(stridingInfo), stridedBlockId_(stridedBlockId), offset_(offset)
    {
    }

    StridedMap(global_size_t numGlobalElements, const Teuchos::ArrayView< const GlobalOrdinal > &elementList, GlobalOrdinal indexBase, std::vector<size_t>& stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1)
    : stridingInfo_(stridingInfo), stridedBlockId_(stridedBlockId)
    {
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo.size() == 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: stridingInfo.size() = 0?");
      TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId < -1, Exceptions::RuntimeError, "StridedMap::StridedMap: stridedBlockId must not be smaller than -1.");

      // the following tests are not valid if stridedBlockId != 1
      if(stridedBlockId == -1) {
        TEUCHOS_TEST_FOR_EXCEPTION(numGlobalElements % getFixedBlockSize() != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: getFixedBlockSize is not an integer multiple of numGlobalElements.");
        TEUCHOS_TEST_FOR_EXCEPTION(elementList.size() % getFixedBlockSize() != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: getFixedBlockSize is not an integer multiple of elementList.size().");
      } else {
        // numGlobalElements can be -1! FIXME
        //TEUCHOS_TEST_FOR_EXCEPTION(numGlobalElements % stridingInfo[stridedBlockId] != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: stridingBlockInfo[stridedBlockId] is not an integer multiple of numGlobalElements.");
        TEUCHOS_TEST_FOR_EXCEPTION(elementList.size() % stridingInfo[stridedBlockId] != 0, Exceptions::RuntimeError, "StridedMap::StridedMap: stridingInfo not valid: stridingBlockInfo[stridedBlockId] is not an integer multiple of elementList.size().");
      }
       
      // calculate offset_

      // find minimum GID over all procs
      GlobalOrdinal minGidOnCurProc = 99999;  // TODO use scalar traits for max possible gid.
      for(Teuchos_Ordinal k=0; k<elementList.size(); ++k) { // TODO fix occurence of Teuchos_Ordinal
        if(elementList[k] < minGidOnCurProc) minGidOnCurProc = elementList[k];
      }
      Teuchos::reduceAll(*comm, Teuchos::REDUCE_MIN, minGidOnCurProc, Teuchos::outArg(offset_));

      // calculate striding index
      size_t nStridedOffset = 0;
      for(int j=0; j<stridedBlockId; j++) {
        nStridedOffset += stridingInfo[j];
      }
      const GlobalOrdinal goStridedOffset = Teuchos::as<GlobalOrdinal>(nStridedOffset);

      // adapt offset_
      offset_ -= goStridedOffset;
    }
     
    virtual ~StridedMap() { }



    
    std::vector<size_t> getStridingData() const { return stridingInfo_; }

    void setStridingData(std::vector<size_t> stridingInfo) { stridingInfo_ = stridingInfo; }

    size_t getFixedBlockSize() const {
      // sum up size of all strided blocks (= number of dofs per node)
      size_t blkSize = 0;
      std::vector<size_t>::const_iterator it;
      for(it = stridingInfo_.begin(); it != stridingInfo_.end(); ++it) {
        blkSize += *it;
      }
      return blkSize;
    }

    LocalOrdinal getStridedBlockId() const { return stridedBlockId_; }

    bool isStrided() { return stridingInfo_.size() > 1 ? true : false; }

    bool isBlocked() { return getFixedBlockSize() > 1 ? true : false; }
    
    GlobalOrdinal getOffset() const { return offset_; } 
    
    void setOffset( GlobalOrdinal offset ) { offset_ = offset; }
    
    // returns number of strided block id which gid belongs to.
    size_t GID2StridingBlockId( GlobalOrdinal gid ) const {
      GlobalOrdinal tgid = gid - offset_;
      tgid = tgid % getFixedBlockSize();

      size_t nStridedOffset = 0;
      size_t stridedBlockId = 0;
      for(size_t j=0; j<stridingInfo_.size(); j++) {
        nStridedOffset += stridingInfo_[j];
        if(Teuchos::as<size_t>(tgid) < nStridedOffset) {
          stridedBlockId = j;
          break;
        }
      }
      return stridedBlockId;
    }

    /* // function currently not needed but maybe useful
    std::vector<GlobalOrdinal> NodeId2GlobalDofIds(GlobalOrdinal nodeId) const {
      TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_.size() == 0, Exceptions::RuntimeError, "StridedMap::NodeId2GlobalDofIds: stridingInfo not valid: stridingInfo.size() = 0?");
      std::vector<GlobalOrdinal> dofs;
      if(stridedBlockId_ > -1) {
          TEUCHOS_TEST_FOR_EXCEPTION(stridingInfo_[stridedBlockId_] == 0, Exceptions::RuntimeError, "StridedMap::NodeId2GlobalDofIds: stridingInfo not valid: stridingInfo[stridedBlockId] = 0?");

          // determine nStridedOffset
          size_t nStridedOffset = 0;
          for(int j=0; j<stridedBlockId_; j++) {
            nStridedOffset += stridingInfo_[j];
          }

          for(size_t i = 0; i<stridingInfo_[stridedBlockId_]; i++) {
            GlobalOrdinal gid =
                nodeId * Teuchos::as<GlobalOrdinal>(getFixedBlockSize()) +
                offset_ +
                Teuchos::as<GlobalOrdinal>(nStridedOffset) +
                Teuchos::as<GlobalOrdinal>(i);
            dofs.push_back(gid);
          }
      } else {
        for(size_t i = 0; i<getFixedBlockSize(); i++) {
          GlobalOrdinal gid =
              nodeId * Teuchos::as<GlobalOrdinal>(getFixedBlockSize()) +
              offset_ +
              Teuchos::as<GlobalOrdinal>(i);
          dofs.push_back(gid);
        }
      }
      return dofs;
    }*/

  private:
    virtual bool CheckConsistency() = 0;

  protected:
    std::vector<size_t> stridingInfo_;   
    LocalOrdinal stridedBlockId_;        
                                         // stridedBlock == -1: the full map (with all strided block dofs)
                                         // stridedBlock >  -1: only dofs of strided block with index "stridedBlockId" are stored in this map
    GlobalOrdinal offset_;     

  }; // StridedMap class

} // Xpetra namespace

#define XPETRA_STRIDEDMAP_SHORT
#endif // XPETRA_STRIDEDMAP_HPP