Tpetra Matrix/Vector Services Version of the Day
Tpetra_BlockMap_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef TPETRA_BLOCKMAP_DEF_HPP
00030 #define TPETRA_BLOCKMAP_DEF_HPP
00031 
00032 #include "Tpetra_ConfigDefs.hpp"
00033 
00034 #ifdef DOXYGEN_USE_ONLY
00035 #include "Tpetra_BlockMap_decl.hpp"
00036 #endif
00037 
00038 namespace Tpetra {
00039 
00040 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00041 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::BlockMap(
00042     global_size_t numGlobalBlocks,
00043     LocalOrdinal blockSize,
00044     GlobalOrdinal indexBase,
00045     const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
00046     const Teuchos::RCP<Node> &node)
00047  : pointMap_(),
00048    globalNumBlocks_(numGlobalBlocks),
00049    myGlobalBlockIDs_(),
00050    pbuf_firstPointInBlock_(),
00051    view_firstPointInBlock_(),
00052    blockIDsAreContiguous_(true),
00053    constantBlockSize_(blockSize)
00054 {
00055   TEST_FOR_EXCEPTION( blockSize <= 0, std::runtime_error,
00056        "Tpetra::BlockMap::BlockMap ERROR: blockSize must be greater than 0.");
00057 
00058   global_size_t numGlobalPoints = numGlobalBlocks*blockSize;
00059 
00060   //we compute numLocalPoints to make sure that Tpetra::Map doesn't split
00061   //numGlobalPoints in a way that would separate the points within a block
00062   //onto different mpi processors.
00063 
00064   size_t numLocalPoints = numGlobalBlocks/comm->getSize();
00065   int remainder = numGlobalBlocks%comm->getSize();
00066   int localProc = comm->getRank();
00067   if (localProc < remainder) ++numLocalPoints;
00068   numLocalPoints *= blockSize;
00069 
00070   //now create the point-map specifying both numGlobalPoints and numLocalPoints:
00071   pointMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(numGlobalPoints, numLocalPoints, indexBase, comm, node));
00072 
00073   size_t numLocalBlocks = pointMap_->getNodeNumElements()/blockSize;
00074   size_t checkLocalBlocks = numLocalPoints/blockSize;
00075   //can there be an inconsistency here???
00076   TEST_FOR_EXCEPTION(numLocalBlocks != checkLocalBlocks, std::runtime_error,
00077        "Tpetra::BlockMap::BlockMap ERROR: internal failure, numLocalBlocks not consistent with point-map.");
00078   
00079   myGlobalBlockIDs_.resize(numLocalBlocks);
00080   pbuf_firstPointInBlock_ = node->template allocBuffer<LocalOrdinal>(numLocalBlocks+1);
00081   Teuchos::ArrayRCP<LocalOrdinal> v_firstPoints = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, numLocalBlocks+1, pbuf_firstPointInBlock_);
00082 
00083   LocalOrdinal firstPoint = pointMap_->getMinLocalIndex();
00084   GlobalOrdinal blockID = pointMap_->getMinGlobalIndex()/blockSize;
00085   for(size_t i=0; i<numLocalBlocks; ++i) {
00086     myGlobalBlockIDs_[i] = blockID++;
00087     v_firstPoints[i] = firstPoint;
00088     firstPoint += blockSize;
00089   }
00090   v_firstPoints[numLocalBlocks] = firstPoint;
00091   v_firstPoints = Teuchos::null;
00092   view_firstPointInBlock_ = node->template viewBuffer<LocalOrdinal>(numLocalBlocks+1, pbuf_firstPointInBlock_);
00093 }
00094 
00095 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00096 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::BlockMap(
00097     global_size_t numGlobalBlocks,
00098     size_t numLocalBlocks,
00099     LocalOrdinal blockSize,
00100     GlobalOrdinal indexBase,
00101     const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
00102     const Teuchos::RCP<Node> &node)
00103  : pointMap_(),
00104    globalNumBlocks_(numGlobalBlocks),
00105    myGlobalBlockIDs_(),
00106    pbuf_firstPointInBlock_(),
00107    view_firstPointInBlock_(),
00108    blockIDsAreContiguous_(true),
00109    constantBlockSize_(blockSize)
00110 {
00111   TEST_FOR_EXCEPTION( blockSize <= 0, std::runtime_error,
00112        "Tpetra::BlockMap::BlockMap ERROR: blockSize must be greater than 0.");
00113 
00114   global_size_t numGlobalPoints = numGlobalBlocks*blockSize;
00115   if (numGlobalBlocks == Teuchos::OrdinalTraits<global_size_t>::invalid()) {
00116     numGlobalPoints = Teuchos::OrdinalTraits<global_size_t>::invalid();
00117   }
00118 
00119   size_t numLocalPoints = numLocalBlocks*blockSize;
00120 
00121   //now create the point-map specifying both numGlobalPoints and numLocalPoints:
00122   pointMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(numGlobalPoints, numLocalPoints, indexBase, comm, node));
00123 
00124   myGlobalBlockIDs_.resize(numLocalBlocks);
00125   pbuf_firstPointInBlock_ = node->template allocBuffer<LocalOrdinal>(numLocalBlocks+1);
00126   Teuchos::ArrayRCP<LocalOrdinal> v_firstPoints = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, numLocalBlocks+1, pbuf_firstPointInBlock_);
00127 
00128   LocalOrdinal firstPoint = pointMap_->getMinLocalIndex();
00129   GlobalOrdinal blockID = pointMap_->getMinGlobalIndex()/blockSize;
00130   for(size_t i=0; i<numLocalBlocks; ++i) {
00131     myGlobalBlockIDs_[i] = blockID++;
00132     v_firstPoints[i] = firstPoint;
00133     firstPoint += blockSize;
00134   }
00135   v_firstPoints[numLocalBlocks] = firstPoint;
00136   v_firstPoints = Teuchos::null;
00137   view_firstPointInBlock_ = node->template viewBuffer<LocalOrdinal>(numLocalBlocks+1, pbuf_firstPointInBlock_);
00138 }
00139 
00140 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00141 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::BlockMap(
00142     global_size_t numGlobalBlocks,
00143     const Teuchos::ArrayView<const GlobalOrdinal>& myGlobalBlockIDs,
00144     const Teuchos::ArrayView<const GlobalOrdinal>& myFirstGlobalPointInBlocks,
00145     const Teuchos::ArrayView<const LocalOrdinal>& blockSizes,
00146     GlobalOrdinal indexBase,
00147     const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
00148     const Teuchos::RCP<Node> &node)
00149  : pointMap_(),
00150    globalNumBlocks_(numGlobalBlocks),
00151    myGlobalBlockIDs_(myGlobalBlockIDs),
00152    pbuf_firstPointInBlock_(),
00153    view_firstPointInBlock_(),
00154    blockIDsAreContiguous_(false),
00155    constantBlockSize_(0)
00156 {
00157   TEST_FOR_EXCEPTION(myGlobalBlockIDs_.size()!=blockSizes.size(), std::runtime_error,
00158              "Tpetra::BlockMap::BlockMap ERROR: input myGlobalBlockIDs and blockSizes arrays must have the same length.");
00159 
00160   size_t sum_blockSizes = 0;
00161   Teuchos::Array<GlobalOrdinal> myGlobalPoints;
00162   typename Teuchos::ArrayView<const LocalOrdinal>::const_iterator
00163     iter = blockSizes.begin(), iend = blockSizes.end();
00164   size_t i = 0;
00165   for(; iter!=iend; ++iter) {
00166     LocalOrdinal bsize = *iter;
00167     sum_blockSizes += bsize;
00168     GlobalOrdinal firstPoint = myFirstGlobalPointInBlocks[i++];
00169     for(LocalOrdinal j=0; j<bsize; ++j) {
00170       myGlobalPoints.push_back(firstPoint+j);
00171     }
00172   }
00173 
00174   pointMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(Teuchos::OrdinalTraits<global_size_t>::invalid(), myGlobalPoints(), indexBase, comm, node));
00175 
00176   iter = blockSizes.begin();
00177   LocalOrdinal firstBlockSize = *iter;
00178   LocalOrdinal firstPoint = pointMap_->getMinLocalIndex();
00179   LocalOrdinal numLocalBlocks = myGlobalBlockIDs.size();
00180   pbuf_firstPointInBlock_ = node->template allocBuffer<LocalOrdinal>(numLocalBlocks+1);
00181   Teuchos::ArrayRCP<LocalOrdinal> v_firstPoints = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, numLocalBlocks+1, pbuf_firstPointInBlock_);
00182 
00183   bool blockSizesAreConstant = true;
00184   i=0;
00185   for(; iter!=iend; ++iter, ++i) {
00186     v_firstPoints[i] = firstPoint;
00187     firstPoint += *iter;
00188     if (*iter != firstBlockSize) {
00189       blockSizesAreConstant = false;
00190     }
00191   }
00192   v_firstPoints[i] = firstPoint;
00193   v_firstPoints = Teuchos::null;
00194   if (blockSizesAreConstant) constantBlockSize_ = firstBlockSize;
00195 
00196   size_t num_points = pointMap_->getNodeNumElements();
00197   TEST_FOR_EXCEPTION(sum_blockSizes != num_points, std::runtime_error,
00198             "Tpetra::BlockMap::BlockMap ERROR: internal failure, sum of block-sizes must equal pointMap->getNodeNumElements().");
00199 
00200   typename Teuchos::Array<GlobalOrdinal>::const_iterator
00201     b_iter = myGlobalBlockIDs_.begin(), b_end = myGlobalBlockIDs_.end();
00202   GlobalOrdinal id = *b_iter;
00203   ++b_iter;
00204   blockIDsAreContiguous_ = true;
00205   for(; b_iter != b_end; ++b_iter) {
00206     if (*b_iter != id+1) {
00207       blockIDsAreContiguous_ = false;
00208       break;
00209     }
00210     ++id;
00211   }
00212   if (blockIDsAreContiguous_ == false) {
00213     setup_noncontig_mapping();
00214   }
00215 
00216   view_firstPointInBlock_ = node->template viewBuffer<LocalOrdinal>(numLocalBlocks+1, pbuf_firstPointInBlock_);
00217 }
00218 
00219 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00220 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::BlockMap(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& pointMap, const Teuchos::ArrayView<const GlobalOrdinal>& myGlobalBlockIDs, const Teuchos::ArrayView<const LocalOrdinal>& blockSizes, const Teuchos::RCP<Node>& node)
00221  : pointMap_(pointMap),
00222    globalNumBlocks_(Teuchos::OrdinalTraits<LocalOrdinal>::invalid()),
00223    myGlobalBlockIDs_(myGlobalBlockIDs),
00224    pbuf_firstPointInBlock_(),
00225    view_firstPointInBlock_(),
00226    blockIDsAreContiguous_(false),
00227    constantBlockSize_(0)
00228 {
00229   TEST_FOR_EXCEPTION(myGlobalBlockIDs_.size()!=blockSizes.size(), std::runtime_error,
00230              "Tpetra::BlockMap::BlockMap ERROR: input myGlobalBlockIDs and blockSizes arrays must have the same length.");
00231 
00232   global_size_t numLocalBlocks = myGlobalBlockIDs.size();
00233   Teuchos::reduceAll<int,global_size_t>(*pointMap->getComm(), Teuchos::REDUCE_SUM, 1, &numLocalBlocks, &globalNumBlocks_);
00234 
00235   LocalOrdinal firstPoint = pointMap->getMinLocalIndex();
00236   size_t sum_blockSizes = 0;
00237   typename Teuchos::ArrayView<const LocalOrdinal>::const_iterator
00238     iter = blockSizes.begin(), iend = blockSizes.end();
00239   LocalOrdinal firstBlockSize = *iter;
00240 
00241   pbuf_firstPointInBlock_ = node->template allocBuffer<LocalOrdinal>(myGlobalBlockIDs.size()+1);
00242   Teuchos::ArrayRCP<LocalOrdinal> v_firstPoints = node->template viewBufferNonConst<LocalOrdinal>(Kokkos::WriteOnly, numLocalBlocks+1, pbuf_firstPointInBlock_);
00243 
00244   bool blockSizesAreConstant = true;
00245   size_t i=0;
00246   for(; iter!=iend; ++iter, ++i) {
00247     sum_blockSizes += *iter;
00248     v_firstPoints[i] = firstPoint;
00249     firstPoint += *iter;
00250     if (*iter != firstBlockSize) {
00251       blockSizesAreConstant = false;
00252     }
00253   }
00254   v_firstPoints[i] = firstPoint;
00255   v_firstPoints = Teuchos::null;
00256   if (blockSizesAreConstant) constantBlockSize_ = firstBlockSize;
00257 
00258   size_t num_points = pointMap->getNodeNumElements();
00259   TEST_FOR_EXCEPTION(sum_blockSizes != num_points, std::runtime_error,
00260             "Tpetra::BlockMap::BlockMap ERROR: sum of block-sizes must equal pointMap->getNodeNumElements().");
00261 
00262   typename Teuchos::Array<GlobalOrdinal>::const_iterator
00263     b_iter = myGlobalBlockIDs_.begin(), b_end = myGlobalBlockIDs_.end();
00264   GlobalOrdinal id = *b_iter;
00265   ++b_iter;
00266   blockIDsAreContiguous_ = true;
00267   for(; b_iter != b_end; ++b_iter) {
00268     if (*b_iter != id+1) {
00269       blockIDsAreContiguous_ = false;
00270       break;
00271     }
00272     ++id;
00273   }
00274   if (blockIDsAreContiguous_ == false) {
00275     setup_noncontig_mapping();
00276   }
00277 
00278   view_firstPointInBlock_ = node->template viewBuffer<LocalOrdinal>(numLocalBlocks+1, pbuf_firstPointInBlock_);
00279 }
00280 
00281 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00282 void
00283 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::setup_noncontig_mapping()
00284 {
00285   //ouch, need to use a hash (unordered_map) here...
00286   typedef typename Teuchos::Array<GlobalOrdinal>::size_type Tsize_t;
00287   for(Tsize_t i=0; i<myGlobalBlockIDs_.size(); ++i) {
00288     LocalOrdinal li = i;
00289     map_global_to_local_.insert(std::make_pair(myGlobalBlockIDs_[i],li));
00290   }
00291 }
00292 
00293 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00294 global_size_t
00295 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getGlobalNumBlocks() const
00296 {
00297   return globalNumBlocks_;
00298 }
00299 
00300 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00301 size_t
00302 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getNodeNumBlocks() const
00303 {
00304   return myGlobalBlockIDs_.size();
00305 }
00306 
00307 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00308 Teuchos::ArrayView<const GlobalOrdinal>
00309 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getNodeBlockIDs() const
00310 {
00311   return myGlobalBlockIDs_();
00312 }
00313 
00314 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00315 bool
00316 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::isBlockSizeConstant() const
00317 { return constantBlockSize_ != 0; }
00318 
00319 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00320 Teuchos::ArrayRCP<const LocalOrdinal>
00321 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getNodeFirstPointInBlocks() const
00322 {
00323   return view_firstPointInBlock_;
00324 }
00325 
00326 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00327 Teuchos::ArrayRCP<const LocalOrdinal>
00328 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getNodeFirstPointInBlocks_Device() const
00329 {
00330   return pbuf_firstPointInBlock_;
00331 }
00332 
00333 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00334 GlobalOrdinal
00335 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getGlobalBlockID(LocalOrdinal localBlockID) const
00336 {
00337   LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00338   if (localBlockID < 0 || localBlockID >= myGlobalBlockIDs_.size()) {
00339     return invalid;
00340   }
00341 
00342   return myGlobalBlockIDs_[localBlockID];
00343 }
00344 
00345 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00346 LocalOrdinal
00347 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getLocalBlockID(GlobalOrdinal globalBlockID) const
00348 {
00349   LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00350   if (myGlobalBlockIDs_.size() == 0) {
00351     return invalid;
00352   }
00353 
00354   if (blockIDsAreContiguous_ == false) {
00355     //need to use a hash (unordered_map) here instead of a map...
00356     typename std::map<GlobalOrdinal,LocalOrdinal>::const_iterator iter =
00357       map_global_to_local_.find(globalBlockID);
00358     if (iter == map_global_to_local_.end()) return invalid;
00359     return iter->second;
00360   }
00361 
00362   LocalOrdinal localBlockID = globalBlockID - myGlobalBlockIDs_[0];
00363   LocalOrdinal numLocalBlocks = myGlobalBlockIDs_.size();
00364 
00365   if (localBlockID < 0 || localBlockID >= numLocalBlocks) {
00366     return invalid;
00367   }
00368 
00369   return localBlockID;
00370 }
00371 
00372 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00373 LocalOrdinal
00374 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getLocalBlockSize(LocalOrdinal localBlockID) const
00375 {
00376   if (constantBlockSize_ != 0) {
00377     return constantBlockSize_;
00378   }
00379 
00380   //should this be a debug-mode-only range check?
00381   if (localBlockID < 0 || localBlockID >= view_firstPointInBlock_.size()) {
00382     throw std::runtime_error("Tpetra::BlockMap::getLocalBlockSize ERROR: localBlockID out of range.");
00383   }
00384 
00385   return view_firstPointInBlock_[localBlockID+1]-view_firstPointInBlock_[localBlockID];
00386 }
00387 
00388 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00389 LocalOrdinal
00390 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getFirstLocalPointInLocalBlock(LocalOrdinal localBlockID) const
00391 {
00392   //should this be a debug-mode-only range check?
00393   if (localBlockID < 0 || localBlockID >= view_firstPointInBlock_.size()) {
00394     throw std::runtime_error("Tpetra::BlockMap::getFirstLocalPointInLocalBlock ERROR: localBlockID out of range.");
00395   }
00396 
00397   return view_firstPointInBlock_[localBlockID];
00398 }
00399 
00400 template<class LocalOrdinal,class GlobalOrdinal,class Node>
00401 GlobalOrdinal
00402 BlockMap<LocalOrdinal,GlobalOrdinal,Node>::getFirstGlobalPointInLocalBlock(LocalOrdinal localBlockID) const
00403 {
00404   //should this be a debug-mode-only range check?
00405   if (localBlockID < 0 || localBlockID >= view_firstPointInBlock_.size()) {
00406     throw std::runtime_error("Tpetra::BlockMap::getFirstGlobalPointInLocalBlock ERROR: localBlockID out of range.");
00407   }
00408 
00409   return pointMap_->getGlobalElement(view_firstPointInBlock_[localBlockID]);
00410 }
00411 
00412 //
00413 // Explicit instantiation macro
00414 //
00415 // Must be expanded from within the Tpetra namespace!
00416 //
00417 
00418 #define TPETRA_BLOCKMAP_INSTANT(LO,GO,NODE) \
00419   \
00420   template class BlockMap< LO , GO , NODE >;
00421 
00422 
00423 }//namespace Tpetra
00424 
00425 #endif
00426 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines