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