Tpetra Matrix/Vector Services Version of the Day
Tpetra_Directory_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_DIRECTORY_HPP
00043 #define TPETRA_DIRECTORY_HPP
00044 
00045 #include <Teuchos_as.hpp>
00046 #include "Tpetra_ConfigDefs.hpp"
00047 #include "Tpetra_Distributor.hpp"
00048 #include "Tpetra_Map.hpp"
00049 
00050 #ifdef DOXYGEN_USE_ONLY
00051   #include "Tpetra_Directory_decl.hpp"
00052 #endif
00053 
00054 namespace Tpetra {
00055 
00056   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00057   Directory<LocalOrdinal,GlobalOrdinal,Node>::Directory(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map_in)
00058   : map_(map_in) 
00059   {
00060     // initialize Comm instance
00061     comm_ = map_->getComm();
00062 
00063     // A directory is not necessary for a non-global ES.
00064     if (map_->isDistributed()) {
00065       // If map_ is contiguously allocated, we can construct the 
00066       // directory from the minMyGID value from each image.
00067       if (map_->isContiguous()) {
00068         // make room for the min on each proc, plus one entry at the end for the max cap
00069         allMinGIDs_.resize(comm_->getSize() + 1);
00070         // get my min
00071         GlobalOrdinal minMyGID = map_->getMinGlobalIndex();
00072         // gather all of the mins into the first getSize() entries of allMinDIGs_
00073         Teuchos::gatherAll<int,GlobalOrdinal>(*comm_,1,&minMyGID,comm_->getSize(),&allMinGIDs_.front());
00074         // put the max cap at the end
00075         allMinGIDs_.back() = map_->getMaxAllGlobalIndex() + Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00076       }
00077       // Otherwise we have to generate the directory using MPI calls.
00078       else {
00079         generateDirectory();
00080       }
00081     }
00082   }
00083 
00084   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00085   Directory<LocalOrdinal,GlobalOrdinal,Node>::~Directory() {}
00086 
00087   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00088   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getDirectoryEntries(
00089               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00090               const Teuchos::ArrayView<int> &nodeIDs) const {
00091     const bool computeLIDs = false;
00092     return getEntries(globalIDs, nodeIDs, Teuchos::null, computeLIDs);
00093   }
00094 
00095   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00096   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getDirectoryEntries(
00097               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00098               const Teuchos::ArrayView<int> &nodeIDs, 
00099               const Teuchos::ArrayView<LocalOrdinal> &localIDs) const {
00100     const bool computeLIDs = true;
00101     return getEntries(globalIDs, nodeIDs, localIDs, computeLIDs);
00102   }
00103 
00104   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00105   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getEntries(
00106               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00107               const Teuchos::ArrayView<int> &nodeIDs, 
00108               const Teuchos::ArrayView<LocalOrdinal> &localIDs, 
00109               bool computeLIDs) const {
00110     const LocalOrdinal LINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00111 
00112     LookupStatus res = AllIDsPresent;
00113 
00114     // fill nodeIDs and localIDs with -1s
00115     TEUCHOS_TEST_FOR_EXCEPTION(nodeIDs.size() != globalIDs.size(), std::runtime_error,
00116         Teuchos::typeName(*this) << "::getEntries(): Output buffers are not allocated properly.");
00117     std::fill(nodeIDs.begin(),nodeIDs.end(), -1);
00118     if (computeLIDs) {
00119       TEUCHOS_TEST_FOR_EXCEPTION(localIDs.size() != globalIDs.size(), std::runtime_error,
00120           Teuchos::typeName(*this) << "::getEntries(): Output buffers are not allocated properly.");
00121       std::fill(localIDs.begin(),localIDs.end(), LINVALID);
00122     }
00123 
00124     const int numImages  = comm_->getSize();
00125     const int myImageID  = comm_->getRank();
00126     const size_t numEntries = globalIDs.size();
00127     const global_size_t nOverP = map_->getGlobalNumElements() / numImages;
00128 
00129     if (map_->isDistributed() == false) {
00130       // Easiest case: Map is serial or locally-replicated
00131       typename Teuchos::ArrayView<int>::iterator imgptr = nodeIDs.begin();
00132       typename Teuchos::ArrayView<LocalOrdinal>::iterator lidptr = localIDs.begin();
00133       typename Teuchos::ArrayView<const GlobalOrdinal>::iterator gid;
00134       for (gid = globalIDs.begin(); gid != globalIDs.end(); ++gid) 
00135       {
00136         if (map_->isNodeGlobalElement(*gid)) {
00137           *imgptr++ = myImageID;
00138           if (computeLIDs) {
00139             *lidptr++ = map_->getLocalElement(*gid);
00140           }
00141         }
00142         else {
00143           // advance the pointers, leaving these values set to invalid
00144           imgptr++;
00145           if (computeLIDs) {
00146             lidptr++;
00147           }
00148           res = IDNotPresent;
00149         }
00150       }
00151     }
00152     else if (map_->isContiguous()) {
00153       // Next Easiest Case: Map is distributed but allocated contiguously
00154       typename Teuchos::ArrayView<int>::iterator imgptr = nodeIDs.begin();
00155       typename Teuchos::ArrayView<LocalOrdinal>::iterator lidptr = localIDs.begin();
00156       typename Teuchos::ArrayView<const GlobalOrdinal>::iterator gid;
00157       for (gid = globalIDs.begin(); gid != globalIDs.end(); ++gid) 
00158       {
00159         LocalOrdinal LID = LINVALID; // Assume not found
00160         int image = -1;
00161         GlobalOrdinal GID = *gid;
00162         // Guess uniform distribution and start a little above it
00163         // TODO: replace by a binary search
00164   //
00165   // The commented-out line below is not correct for negative
00166   // GID values.  This is because a GlobalOrdinal (typically
00167   // signed) divided by a global_size_t (unsigned) turns into
00168   // unsigned, then gets cast back to (signed) int.  This turns
00169   // small negative values of GID into bogus values.  We
00170   // observed this in practice.
00171   //
00172         //int curimg = TEUCHOS_MIN((int)(GID / TEUCHOS_MAX(nOverP, 1)) + 2, numImages - 1);
00173   int curimg;
00174   {
00175     const GlobalOrdinal one = Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00176     const GlobalOrdinal two = one + one;
00177     const GlobalOrdinal nOverP_GID = static_cast<GlobalOrdinal> (nOverP);
00178     const GlobalOrdinal lowerBound = GID / TEUCHOS_MAX(nOverP_GID, one) + two;
00179     // It's probably not OK to cast this to int in general.  It
00180     // works as long as |GID| <= the global number of entries
00181     // and nOverP is appropriately sized for int.  Trouble may
00182     // ensue if the index base has an exotic value.
00183     const int lowerBound_int = static_cast<int> (lowerBound);
00184     curimg = TEUCHOS_MIN(lowerBound_int, numImages - 1);
00185   }
00186         bool found = false;
00187         while (curimg >= 0 && curimg < numImages) {
00188           if (allMinGIDs_[curimg] <= GID) {
00189             if (GID < allMinGIDs_[curimg + 1]) {
00190               found = true;
00191               break;
00192             }
00193             else {
00194               curimg++;
00195             }
00196           }
00197           else {
00198             curimg--;
00199           }
00200         }
00201         if (found) {
00202           image = curimg;
00203           LID = Teuchos::as<LocalOrdinal>(GID - allMinGIDs_[image]);
00204         }
00205         else {
00206           res = IDNotPresent;
00207         }
00208         *imgptr++ = image;
00209         if (computeLIDs) {
00210           *lidptr++ = LID;
00211         }
00212       }
00213     }
00214     else {
00215       // General Case: Map is distributed and allocated arbitrarily
00216       // Here we need to set up an actual directory structure
00217       using Teuchos::as;
00218       int packetSize = 2;
00219       if (computeLIDs) {
00220         packetSize = 3;
00221       }
00222 
00223       Distributor distor(comm_);
00224 
00225       // Get directory locations for the requested list of entries
00226       Teuchos::Array<int> dirImages(numEntries);
00227       res = directoryMap_->getRemoteIndexList(globalIDs, dirImages());
00228       // Check for unfound globalIDs and set corresponding nodeIDs to -1
00229       size_t numMissing = 0;
00230       if (res == IDNotPresent) {
00231         for (size_t i=0; i < numEntries; ++i) {
00232           if (dirImages[i] == -1) {
00233             nodeIDs[i] = -1;
00234             if (computeLIDs) {
00235               localIDs[i] = LINVALID;
00236             }
00237             numMissing++;
00238           }
00239         }
00240       }
00241 
00242       Teuchos::ArrayRCP<GlobalOrdinal> sendGIDs; 
00243       Teuchos::ArrayRCP<int> sendImages;
00244       distor.createFromRecvs(globalIDs, dirImages(), sendGIDs, sendImages);
00245       size_t numSends = sendGIDs.size();
00246 
00247       //    global_size_t >= GlobalOrdinal
00248       //    global_size_t >= size_t >= int
00249       //    global_size_t >= size_t >= LocalOrdinal
00250       // Therefore, we can safely stored all of these in a global_size_t
00251       Teuchos::Array<global_size_t> exports(packetSize*numSends);
00252       {
00253         LocalOrdinal curLID;
00254         typename Teuchos::Array<global_size_t>::iterator ptr = exports.begin();
00255         typename Teuchos::ArrayRCP<GlobalOrdinal>::const_iterator gidptr;
00256         for (gidptr = sendGIDs.begin(); gidptr != sendGIDs.end(); ++gidptr) {
00257           *ptr++ = as<global_size_t>(*gidptr);
00258           curLID = directoryMap_->getLocalElement(*gidptr);
00259           TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
00260               Teuchos::typeName(*this) << "::getEntries(): Internal logic error. Please contact Tpetra team.");
00261           *ptr++ = as<global_size_t>(nodeIDs_[curLID]);
00262           if (computeLIDs) {
00263             *ptr++ = as<global_size_t>(LIDs_[curLID]);
00264           }
00265         }
00266       }
00267 
00268       Teuchos::Array<global_size_t> imports(packetSize*distor.getTotalReceiveLength());
00269       distor.doPostsAndWaits(exports().getConst(), packetSize, imports());
00270 
00271       typename Teuchos::Array<global_size_t>::iterator ptr = imports.begin();
00272       const size_t numRecv = numEntries - numMissing;
00273 
00274       Teuchos::Array<GlobalOrdinal> sortedIDs(globalIDs);
00275       Teuchos::ArrayRCP<GlobalOrdinal> offset = arcp<GlobalOrdinal>(numEntries);
00276       GlobalOrdinal ii=0;
00277       for (typename Teuchos::ArrayRCP<GlobalOrdinal>::iterator oo = offset.begin(); oo != offset.end(); ++oo,++ii)
00278         *oo = ii;
00279       sort2(sortedIDs.begin(),sortedIDs.begin()+numEntries,offset.begin());
00280 
00281       typedef typename Teuchos::Array<GlobalOrdinal>::iterator IT;
00282       // we know these conversions are in range, because we loaded this data
00283       for (size_t i = 0; i < numRecv; ++i) {
00284         GlobalOrdinal curGID = as<GlobalOrdinal>(*ptr++);
00285         std::pair< IT, IT> p1 = std::equal_range(sortedIDs.begin(),sortedIDs.end(),curGID);
00286         if (p1.first != p1.second) {
00287           //found it
00288           size_t j = p1.first - sortedIDs.begin();
00289           nodeIDs[offset[j]] = as<int>(*ptr++);
00290           if (computeLIDs) {
00291             localIDs[offset[j]] = as<LocalOrdinal>(*ptr++);
00292           }
00293           if (nodeIDs[offset[j]] == -1) res = IDNotPresent;
00294         }
00295       }
00296     }
00297     return res;
00298   }
00299 
00300 
00301   // directory setup for non-contiguous Map
00302   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00303   void Directory<LocalOrdinal,GlobalOrdinal,Node>::generateDirectory() {
00304     using Teuchos::as;
00305     const LocalOrdinal  LINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00306           
00307     const GlobalOrdinal minAllGID = map_->getMinAllGlobalIndex();
00308     const GlobalOrdinal maxAllGID = map_->getMaxAllGlobalIndex();
00309 
00310     // DirectoryMap will have a range of elements from the minimum to the maximum
00311     // GID of the user Map, and an indexBase of minAllGID from the user Map
00312     const global_size_t numGlobalEntries = maxAllGID - minAllGID + 1;
00313 
00314     // Obviously, we can't afford to store the whole directory on each node
00315     // Create a uniform linear map to contain the directory to split up the storage among all nodes
00316     directoryMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(numGlobalEntries, minAllGID, comm_, GloballyDistributed, map_->getNode() ));
00317 
00318     const size_t dir_numMyEntries = directoryMap_->getNodeNumElements();
00319 
00320     // Allocate imageID List and LID List.  Initialize to -1s.
00321     // Initialize values to -1 in case the user global element list does
00322     // fill all IDs from minAllGID to maxAllGID (e.g., allows global indices to be 
00323     // all even integers).
00324     nodeIDs_.resize(dir_numMyEntries, -1);
00325     LIDs_.resize(dir_numMyEntries, LINVALID);
00326 
00327     // Get list of nodeIDs owning the directory entries for the Map GIDs
00328     const int myImageID = comm_->getRank();
00329     const size_t numMyEntries = map_->getNodeNumElements();
00330     Teuchos::Array<int> sendImageIDs(numMyEntries);
00331     Teuchos::ArrayView<const GlobalOrdinal> myGlobalEntries = map_->getNodeElementList();
00332     // an ID not present in this lookup indicates that it lies outside of the range [minAllGID,maxAllGID] (from map_). 
00333     // this means something is wrong with map_, our fault.
00334     TEUCHOS_TEST_FOR_EXCEPTION( directoryMap_->getRemoteIndexList(myGlobalEntries, sendImageIDs) == IDNotPresent, std::logic_error,
00335         Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please contact Tpetra team.");
00336 
00337     // Create distributor & call createFromSends
00338     size_t numReceives = 0;
00339     Distributor distor(comm_);      
00340     numReceives = distor.createFromSends(sendImageIDs);
00341 
00342     // Execute distributor plan
00343     // Transfer GIDs, ImageIDs, and LIDs that we own to all nodeIDs
00344     // End result is all nodeIDs have list of all GIDs and corresponding ImageIDs and LIDs
00345     int packetSize = 3; // We will send GIDs, ImageIDs, and LIDs.
00346     // GlobalOrdinal >= LocalOrdinal and int, so this is safe
00347     Teuchos::Array<GlobalOrdinal> exportEntries(packetSize*numMyEntries);
00348     {
00349       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = exportEntries.begin();
00350       for (size_t i=0; i < numMyEntries; ++i) {
00351         *ptr++ = myGlobalEntries[i];
00352         *ptr++ = as<GlobalOrdinal>(myImageID);
00353         *ptr++ = as<GlobalOrdinal>(i);
00354       }
00355     }
00356 
00357     Teuchos::Array<GlobalOrdinal> importElements(packetSize*distor.getTotalReceiveLength());
00358     distor.doPostsAndWaits(exportEntries().getConst(), packetSize, importElements());
00359 
00360     {// begin scoping block
00361       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = importElements.begin();
00362       for (size_t i = 0; i < numReceives; ++i) {
00363   // Each "packet" (contiguous chunk of importElements) contains
00364   // a triple: (GID, imageID, LID).
00365         LocalOrdinal currLID = directoryMap_->getLocalElement(*ptr++); // Convert incoming GID to Directory LID
00366         TEUCHOS_TEST_FOR_EXCEPTION(currLID == LINVALID, std::logic_error,
00367             Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please notify the Tpetra team.");
00368         nodeIDs_[currLID] = *ptr++;
00369         LIDs_[currLID] = *ptr++;
00370       }
00371     }
00372   } // end generateDirectory()
00373     
00374 } // namespace Tpetra
00375 
00376 //
00377 // Explicit instantiation macro
00378 //
00379 // Must be expanded from within the Tpetra namespace!
00380 //
00381 
00382 #define TPETRA_DIRECTORY_INSTANT(LO,GO,NODE) \
00383   \
00384   template class Directory< LO , GO , NODE >; \
00385 
00386 #endif // TPETRA_DIRECTORY_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines