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