Tpetra_Directory_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 // FINISH: HERE: Directory should take an RCP<const Map> in constructor
00030 
00031 #ifndef TPETRA_DIRECTORY_HPP
00032 #define TPETRA_DIRECTORY_HPP
00033 
00034 #include <Teuchos_TypeNameTraits.hpp>
00035 #include <Teuchos_OrdinalTraits.hpp>
00036 #include <Teuchos_CommHelpers.hpp>
00037 #include <Teuchos_as.hpp>
00038 #include "Tpetra_ConfigDefs.hpp"
00039 #include "Tpetra_Distributor.hpp"
00040 #include "Tpetra_Map.hpp"
00041 
00042 namespace Tpetra {
00043 
00044   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00045   Directory<LocalOrdinal,GlobalOrdinal,Node>::Directory(const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map_in)
00046   : map_(map_in) {
00047     // initialize Comm instance
00048     comm_ = map_->getComm();
00049 
00050     // A directory is not necessary for a non-global ES.
00051     if (map_->isDistributed()) {
00052       // If map_ is contiguously allocated, we can construct the 
00053       // directory from the minMyGID value from each image.
00054       if (map_->isContiguous()) {
00055         // make room for the min on each proc, plus one entry at the end for the max cap
00056         allMinGIDs_.resize(comm_->getSize() + 1);
00057         // get my min
00058         GlobalOrdinal minMyGID = map_->getMinGlobalIndex();
00059         // gather all of the mins into the first getSize() entries of allMinDIGs_
00060         Teuchos::gatherAll<int,GlobalOrdinal>(*comm_,1,&minMyGID,comm_->getSize(),&allMinGIDs_.front());
00061         // put the max cap at the end
00062         allMinGIDs_.back() = map_->getMaxAllGlobalIndex() + Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00063       }
00064       // Otherwise we have to generate the directory using MPI calls.
00065       else {
00066         generateDirectory();
00067       }
00068     }
00069   }
00070 
00071   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00072   Directory<LocalOrdinal,GlobalOrdinal,Node>::~Directory() {}
00073 
00074   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00075   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getDirectoryEntries(
00076               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00077               const Teuchos::ArrayView<int> &nodeIDs) const {
00078     const bool computeLIDs = false;
00079     return getEntries(globalIDs, nodeIDs, Teuchos::null, computeLIDs);
00080   }
00081 
00082   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00083   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getDirectoryEntries(
00084               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00085               const Teuchos::ArrayView<int> &nodeIDs, 
00086               const Teuchos::ArrayView<LocalOrdinal> &localIDs) const {
00087     const bool computeLIDs = true;
00088     return getEntries(globalIDs, nodeIDs, localIDs, computeLIDs);
00089   }
00090 
00091   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00092   LookupStatus Directory<LocalOrdinal,GlobalOrdinal,Node>::getEntries(
00093               const Teuchos::ArrayView<const GlobalOrdinal> &globalIDs, 
00094               const Teuchos::ArrayView<int> &nodeIDs, 
00095               const Teuchos::ArrayView<LocalOrdinal> &localIDs, 
00096               bool computeLIDs) const {
00097     const LocalOrdinal LINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00098 
00099     LookupStatus res = AllIDsPresent;
00100 
00101     // fill nodeIDs and localIDs with -1s
00102     TEST_FOR_EXCEPTION(nodeIDs.size() != globalIDs.size(), std::runtime_error,
00103         Teuchos::typeName(*this) << "::getEntries(): Output buffers are not allocated properly.");
00104     std::fill(nodeIDs.begin(),nodeIDs.end(), -1);
00105     if (computeLIDs) {
00106       TEST_FOR_EXCEPTION(localIDs.size() != globalIDs.size(), std::runtime_error,
00107           Teuchos::typeName(*this) << "::getEntries(): Output buffers are not allocated properly.");
00108       std::fill(localIDs.begin(),localIDs.end(), LINVALID);
00109     }
00110 
00111     const int numImages  = comm_->getSize();
00112     const int myImageID  = comm_->getRank();
00113     const size_t numEntries = globalIDs.size();
00114     const global_size_t nOverP = map_->getGlobalNumElements() / numImages;
00115 
00116     if (map_->isDistributed() == false) {
00117       // Easiest case: Map is serial or locally-replicated
00118       typename Teuchos::ArrayView<int>::iterator imgptr = nodeIDs.begin();
00119       typename Teuchos::ArrayView<LocalOrdinal>::iterator lidptr = localIDs.begin();
00120       typename Teuchos::ArrayView<const GlobalOrdinal>::iterator gid;
00121       for (gid = globalIDs.begin(); gid != globalIDs.end(); ++gid) 
00122       {
00123         if (map_->isNodeGlobalElement(*gid)) {
00124           *imgptr++ = myImageID;
00125           if (computeLIDs) {
00126             *lidptr++ = map_->getLocalElement(*gid);
00127           }
00128         }
00129         else {
00130           // advance the pointers, leaving these values set to invalid
00131           imgptr++;
00132           if (computeLIDs) {
00133             lidptr++;
00134           }
00135           res = IDNotPresent;
00136         }
00137       }
00138     }
00139     else if (map_->isContiguous()) {
00140       // Next Easiest Case: Map is distributed but allocated contiguously
00141       typename Teuchos::ArrayView<int>::iterator imgptr = nodeIDs.begin();
00142       typename Teuchos::ArrayView<LocalOrdinal>::iterator lidptr = localIDs.begin();
00143       typename Teuchos::ArrayView<const GlobalOrdinal>::iterator gid;
00144       for (gid = globalIDs.begin(); gid != globalIDs.end(); ++gid) 
00145       {
00146         LocalOrdinal LID = LINVALID; // Assume not found
00147         int image = -1;
00148         GlobalOrdinal GID = *gid;
00149         // Guess uniform distribution and start a little above it
00150         // TODO: replace by a binary search
00151         int curimg = TEUCHOS_MIN((int)(GID / TEUCHOS_MAX(nOverP, 1)) + 2, numImages - 1);
00152         bool found = false;
00153         while (curimg >= 0 && curimg < numImages) {
00154           if (allMinGIDs_[curimg] <= GID) {
00155             if (GID < allMinGIDs_[curimg + 1]) {
00156               found = true;
00157               break;
00158             }
00159             else {
00160               curimg++;
00161             }
00162           }
00163           else {
00164             curimg--;
00165           }
00166         }
00167         if (found) {
00168           image = curimg;
00169           LID = Teuchos::as<LocalOrdinal>(GID - allMinGIDs_[image]);
00170         }
00171         else {
00172           res = IDNotPresent;
00173         }
00174         *imgptr++ = image;
00175         if (computeLIDs) {
00176           *lidptr++ = LID;
00177         }
00178       }
00179     }
00180     else {
00181       // General Case: Map is distributed and allocated arbitrarily
00182       // Here we need to set up an actual directory structure
00183       using Teuchos::as;
00184       int packetSize = 2;
00185       if (computeLIDs) {
00186         packetSize = 3;
00187       }
00188 
00189       Distributor distor(comm_);
00190 
00191       // Get directory locations for the requested list of entries
00192       Teuchos::Array<int> dirImages(numEntries);
00193       res = directoryMap_->getRemoteIndexList(globalIDs, dirImages());
00194       // Check for unfound globalIDs and set corresponding nodeIDs to -1
00195       size_t numMissing = 0;
00196       if (res == IDNotPresent) {
00197         for (size_t i=0; i < numEntries; ++i) {
00198           if (dirImages[i] == -1) {
00199             nodeIDs[i] = -1;
00200             if (computeLIDs) {
00201               localIDs[i] = LINVALID;
00202             }
00203             numMissing++;
00204           }
00205         }
00206       }
00207 
00208       Teuchos::ArrayRCP<GlobalOrdinal> sendGIDs; 
00209       Teuchos::ArrayRCP<int> sendImages;
00210       distor.createFromRecvs(globalIDs, dirImages(), sendGIDs, sendImages);
00211       size_t numSends = sendGIDs.size();
00212 
00213       //    global_size_t >= GlobalOrdinal
00214       //    global_size_t >= size_t >= int
00215       //    global_size_t >= size_t >= LocalOrdinal
00216       // Therefore, we can safely stored all of these in a global_size_t
00217       Teuchos::Array<global_size_t> exports(packetSize*numSends);
00218       {
00219         LocalOrdinal curLID;
00220         typename Teuchos::Array<global_size_t>::iterator ptr = exports.begin();
00221         typename Teuchos::ArrayRCP<GlobalOrdinal>::const_iterator gidptr;
00222         for (gidptr = sendGIDs.begin(); gidptr != sendGIDs.end(); ++gidptr) {
00223           *ptr++ = as<global_size_t>(*gidptr);
00224           curLID = directoryMap_->getLocalElement(*gidptr);
00225           TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
00226               Teuchos::typeName(*this) << "::getEntries(): Internal logic error. Please contact Tpetra team.");
00227           *ptr++ = as<global_size_t>(nodeIDs_[curLID]);
00228           if (computeLIDs) {
00229             *ptr++ = as<global_size_t>(LIDs_[curLID]);
00230           }
00231         }
00232       }
00233 
00234       Teuchos::Array<global_size_t> imports(packetSize*distor.getTotalReceiveLength());
00235       distor.doPostsAndWaits(exports().getConst(), packetSize, imports());
00236 
00237       typename Teuchos::Array<global_size_t>::iterator ptr = imports.begin();
00238       const size_t numRecv = numEntries - numMissing;
00239       // we know these conversions are in range, because we loaded this data
00240       for (size_t i = 0; i < numRecv; ++i) {
00241         GlobalOrdinal curGID = as<GlobalOrdinal>(*ptr++);
00242         for (size_t j = 0; j < numEntries; ++j) {
00243           if (curGID == globalIDs[j]) {
00244             nodeIDs[j] = as<int>(*ptr++);
00245             if (computeLIDs) {
00246               localIDs[j] = as<LocalOrdinal>(*ptr++);
00247             }
00248             break;
00249           }
00250         }
00251       }
00252     }
00253     return res;
00254   }
00255 
00256 
00257   // directory setup for non-contiguous Map
00258   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00259   void Directory<LocalOrdinal,GlobalOrdinal,Node>::generateDirectory() {
00260     using Teuchos::as;
00261     const GlobalOrdinal GONE     = Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00262     const LocalOrdinal  LINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00263           
00264     const GlobalOrdinal minAllGID = map_->getMinAllGlobalIndex();
00265     const GlobalOrdinal maxAllGID = map_->getMaxAllGlobalIndex();
00266 
00267     // DirectoryMap will have a range of elements from the minimum to the maximum
00268     // GID of the user Map, and an indexBase of minAllGID from the user Map
00269     global_size_t numGlobalEntries = maxAllGID - minAllGID + GONE;
00270 
00271     // Obviously, we can't afford to store the whole directory on each node
00272     // Create a uniform linear map to contain the directory to split up the storage among all nodes
00273     directoryMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(numGlobalEntries, minAllGID, comm_, GloballyDistributed, map_->getNode() ));
00274 
00275     const size_t dir_numMyEntries = directoryMap_->getNodeNumElements();
00276 
00277     // Allocate imageID List and LID List.  Initialize to -1s.
00278     // Initialize values to -1 in case the user global element list does
00279     // fill all IDs from minAllGID to maxAllGID (e.g., allows global indices to be 
00280     // all even integers).
00281     nodeIDs_.resize(dir_numMyEntries, -1);
00282     LIDs_.resize(dir_numMyEntries, LINVALID);
00283 
00284     // Get list of nodeIDs owning the directory entries for the Map GIDs
00285     const int myImageID = comm_->getRank();
00286     const size_t numMyEntries = map_->getNodeNumElements();
00287     Teuchos::Array<int> sendImageIDs(numMyEntries);
00288     Teuchos::ArrayView<const GlobalOrdinal> myGlobalEntries = map_->getNodeElementList();
00289     // an ID not present in this lookup indicates that it lies outside of the range [minAllGID,maxAllGID] (from map_). 
00290     // this means something is wrong with map_, our fault.
00291     TEST_FOR_EXCEPTION( directoryMap_->getRemoteIndexList(myGlobalEntries, sendImageIDs) == IDNotPresent, std::logic_error,
00292         Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please contact Tpetra team.");
00293 
00294     // Create distributor & call createFromSends
00295     size_t numReceives = 0;
00296     Distributor distor(comm_);      
00297     numReceives = distor.createFromSends(sendImageIDs);
00298 
00299     // Execute distributor plan
00300     // Transfer GIDs, ImageIDs, and LIDs that we own to all nodeIDs
00301     // End result is all nodeIDs have list of all GIDs and corresponding ImageIDs and LIDs
00302     int packetSize = 3; // We will send GIDs, ImageIDs, and LIDs.
00303     // GlobalOrdinal >= LocalOrdinal and int, so this is safe
00304     Teuchos::Array<GlobalOrdinal> exportEntries(packetSize*numMyEntries);
00305     {
00306       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = exportEntries.begin();
00307       for (size_t i=0; i < numMyEntries; ++i) {
00308         *ptr++ = myGlobalEntries[i];
00309         *ptr++ = as<GlobalOrdinal>(myImageID);
00310         *ptr++ = as<GlobalOrdinal>(i);
00311       }
00312     }
00313 
00314     Teuchos::Array<GlobalOrdinal> importElements(packetSize*distor.getTotalReceiveLength());
00315     distor.doPostsAndWaits(exportEntries().getConst(), packetSize, importElements());
00316 
00317     {// begin scoping block
00318       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = importElements.begin();
00319       for (size_t i = 0; i < numReceives; ++i) {
00320         LocalOrdinal currLID = directoryMap_->getLocalElement(*ptr++); // Convert incoming GID to Directory LID
00321         TEST_FOR_EXCEPTION(currLID == LINVALID, std::logic_error,
00322             Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please notify the Tpetra team.");
00323         nodeIDs_[currLID] = *ptr++;
00324         LIDs_[currLID] = *ptr++;
00325       }
00326     }
00327   } // end generateDirectory()
00328     
00329 } // namespace Tpetra
00330 
00331 //
00332 // Explicit instantiation macro
00333 //
00334 // Must be expanded from within the Tpetra namespace!
00335 //
00336 
00337 #define TPETRA_DIRECTORY_INSTANT(LO,GO,NODE) \
00338   \
00339   template class Directory< LO , GO , NODE >; \
00340 
00341 #endif // TPETRA_DIRECTORY_HPP

Generated on Tue Jul 13 09:39:06 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7