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   //
00151   // The commented-out line below is not correct for negative
00152   // GID values.  This is because a GlobalOrdinal (typically
00153   // signed) divided by a global_size_t (unsigned) turns into
00154   // unsigned, then gets cast back to (signed) int.  This turns
00155   // small negative values of GID into bogus values.  We
00156   // observed this in practice.
00157   //
00158         //int curimg = TEUCHOS_MIN((int)(GID / TEUCHOS_MAX(nOverP, 1)) + 2, numImages - 1);
00159   int curimg;
00160   {
00161     const GlobalOrdinal one = Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00162     const GlobalOrdinal two = one + one;
00163     const GlobalOrdinal nOverP_GID = static_cast<GlobalOrdinal> (nOverP);
00164     const GlobalOrdinal lowerBound = GID / TEUCHOS_MAX(nOverP_GID, one) + two;
00165     // It's probably not OK to cast this to int in general.  It
00166     // works as long as |GID| <= the global number of entries
00167     // and nOverP is appropriately sized for int.  Trouble may
00168     // ensue if the index base has an exotic value.
00169     const int lowerBound_int = static_cast<int> (lowerBound);
00170     curimg = TEUCHOS_MIN(lowerBound_int, numImages - 1);
00171   }
00172         bool found = false;
00173         while (curimg >= 0 && curimg < numImages) {
00174           if (allMinGIDs_[curimg] <= GID) {
00175             if (GID < allMinGIDs_[curimg + 1]) {
00176               found = true;
00177               break;
00178             }
00179             else {
00180               curimg++;
00181             }
00182           }
00183           else {
00184             curimg--;
00185           }
00186         }
00187         if (found) {
00188           image = curimg;
00189           LID = Teuchos::as<LocalOrdinal>(GID - allMinGIDs_[image]);
00190         }
00191         else {
00192           res = IDNotPresent;
00193         }
00194         *imgptr++ = image;
00195         if (computeLIDs) {
00196           *lidptr++ = LID;
00197         }
00198       }
00199     }
00200     else {
00201       // General Case: Map is distributed and allocated arbitrarily
00202       // Here we need to set up an actual directory structure
00203       using Teuchos::as;
00204       int packetSize = 2;
00205       if (computeLIDs) {
00206         packetSize = 3;
00207       }
00208 
00209       Distributor distor(comm_);
00210 
00211       // Get directory locations for the requested list of entries
00212       Teuchos::Array<int> dirImages(numEntries);
00213       res = directoryMap_->getRemoteIndexList(globalIDs, dirImages());
00214       // Check for unfound globalIDs and set corresponding nodeIDs to -1
00215       size_t numMissing = 0;
00216       if (res == IDNotPresent) {
00217         for (size_t i=0; i < numEntries; ++i) {
00218           if (dirImages[i] == -1) {
00219             nodeIDs[i] = -1;
00220             if (computeLIDs) {
00221               localIDs[i] = LINVALID;
00222             }
00223             numMissing++;
00224           }
00225         }
00226       }
00227 
00228       Teuchos::ArrayRCP<GlobalOrdinal> sendGIDs; 
00229       Teuchos::ArrayRCP<int> sendImages;
00230       distor.createFromRecvs(globalIDs, dirImages(), sendGIDs, sendImages);
00231       size_t numSends = sendGIDs.size();
00232 
00233       //    global_size_t >= GlobalOrdinal
00234       //    global_size_t >= size_t >= int
00235       //    global_size_t >= size_t >= LocalOrdinal
00236       // Therefore, we can safely stored all of these in a global_size_t
00237       Teuchos::Array<global_size_t> exports(packetSize*numSends);
00238       {
00239         LocalOrdinal curLID;
00240         typename Teuchos::Array<global_size_t>::iterator ptr = exports.begin();
00241         typename Teuchos::ArrayRCP<GlobalOrdinal>::const_iterator gidptr;
00242         for (gidptr = sendGIDs.begin(); gidptr != sendGIDs.end(); ++gidptr) {
00243           *ptr++ = as<global_size_t>(*gidptr);
00244           curLID = directoryMap_->getLocalElement(*gidptr);
00245           TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
00246               Teuchos::typeName(*this) << "::getEntries(): Internal logic error. Please contact Tpetra team.");
00247           *ptr++ = as<global_size_t>(nodeIDs_[curLID]);
00248           if (computeLIDs) {
00249             *ptr++ = as<global_size_t>(LIDs_[curLID]);
00250           }
00251         }
00252       }
00253 
00254       Teuchos::Array<global_size_t> imports(packetSize*distor.getTotalReceiveLength());
00255       distor.doPostsAndWaits(exports().getConst(), packetSize, imports());
00256 
00257       typename Teuchos::Array<global_size_t>::iterator ptr = imports.begin();
00258       const size_t numRecv = numEntries - numMissing;
00259 
00260       Teuchos::Array<GlobalOrdinal> sortedIDs(globalIDs);
00261       Teuchos::ArrayRCP<GlobalOrdinal> offset = arcp<GlobalOrdinal>(numEntries);
00262       GlobalOrdinal ii=0;
00263       for (typename Teuchos::ArrayRCP<GlobalOrdinal>::iterator oo = offset.begin(); oo != offset.end(); ++oo,++ii)
00264         *oo = ii;
00265       sort2(sortedIDs.begin(),sortedIDs.begin()+numEntries,offset.begin());
00266 
00267       typedef typename Teuchos::Array<GlobalOrdinal>::iterator IT;
00268       // we know these conversions are in range, because we loaded this data
00269       for (size_t i = 0; i < numRecv; ++i) {
00270         GlobalOrdinal curGID = as<GlobalOrdinal>(*ptr++);
00271         std::pair< IT, IT> p1 = std::equal_range(sortedIDs.begin(),sortedIDs.end(),curGID);
00272         if (p1.first != p1.second) {
00273           //found it
00274           size_t j = p1.first - sortedIDs.begin();
00275           nodeIDs[offset[j]] = as<int>(*ptr++);
00276           if (computeLIDs) {
00277             localIDs[offset[j]] = as<LocalOrdinal>(*ptr++);
00278           }
00279           if (nodeIDs[offset[j]] == -1) res = IDNotPresent;
00280         }
00281       }
00282     }
00283     return res;
00284   }
00285 
00286 
00287   // directory setup for non-contiguous Map
00288   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00289   void Directory<LocalOrdinal,GlobalOrdinal,Node>::generateDirectory() {
00290     using Teuchos::as;
00291     const LocalOrdinal  LINVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00292           
00293     const GlobalOrdinal minAllGID = map_->getMinAllGlobalIndex();
00294     const GlobalOrdinal maxAllGID = map_->getMaxAllGlobalIndex();
00295 
00296     // DirectoryMap will have a range of elements from the minimum to the maximum
00297     // GID of the user Map, and an indexBase of minAllGID from the user Map
00298     const global_size_t numGlobalEntries = maxAllGID - minAllGID + 1;
00299 
00300     // Obviously, we can't afford to store the whole directory on each node
00301     // Create a uniform linear map to contain the directory to split up the storage among all nodes
00302     directoryMap_ = Teuchos::rcp(new Map<LocalOrdinal,GlobalOrdinal,Node>(numGlobalEntries, minAllGID, comm_, GloballyDistributed, map_->getNode() ));
00303 
00304     const size_t dir_numMyEntries = directoryMap_->getNodeNumElements();
00305 
00306     // Allocate imageID List and LID List.  Initialize to -1s.
00307     // Initialize values to -1 in case the user global element list does
00308     // fill all IDs from minAllGID to maxAllGID (e.g., allows global indices to be 
00309     // all even integers).
00310     nodeIDs_.resize(dir_numMyEntries, -1);
00311     LIDs_.resize(dir_numMyEntries, LINVALID);
00312 
00313     // Get list of nodeIDs owning the directory entries for the Map GIDs
00314     const int myImageID = comm_->getRank();
00315     const size_t numMyEntries = map_->getNodeNumElements();
00316     Teuchos::Array<int> sendImageIDs(numMyEntries);
00317     Teuchos::ArrayView<const GlobalOrdinal> myGlobalEntries = map_->getNodeElementList();
00318     // an ID not present in this lookup indicates that it lies outside of the range [minAllGID,maxAllGID] (from map_). 
00319     // this means something is wrong with map_, our fault.
00320     TEST_FOR_EXCEPTION( directoryMap_->getRemoteIndexList(myGlobalEntries, sendImageIDs) == IDNotPresent, std::logic_error,
00321         Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please contact Tpetra team.");
00322 
00323     // Create distributor & call createFromSends
00324     size_t numReceives = 0;
00325     Distributor distor(comm_);      
00326     numReceives = distor.createFromSends(sendImageIDs);
00327 
00328     // Execute distributor plan
00329     // Transfer GIDs, ImageIDs, and LIDs that we own to all nodeIDs
00330     // End result is all nodeIDs have list of all GIDs and corresponding ImageIDs and LIDs
00331     int packetSize = 3; // We will send GIDs, ImageIDs, and LIDs.
00332     // GlobalOrdinal >= LocalOrdinal and int, so this is safe
00333     Teuchos::Array<GlobalOrdinal> exportEntries(packetSize*numMyEntries);
00334     {
00335       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = exportEntries.begin();
00336       for (size_t i=0; i < numMyEntries; ++i) {
00337         *ptr++ = myGlobalEntries[i];
00338         *ptr++ = as<GlobalOrdinal>(myImageID);
00339         *ptr++ = as<GlobalOrdinal>(i);
00340       }
00341     }
00342 
00343     Teuchos::Array<GlobalOrdinal> importElements(packetSize*distor.getTotalReceiveLength());
00344     distor.doPostsAndWaits(exportEntries().getConst(), packetSize, importElements());
00345 
00346     {// begin scoping block
00347       typename Teuchos::Array<GlobalOrdinal>::iterator ptr = importElements.begin();
00348       for (size_t i = 0; i < numReceives; ++i) {
00349         LocalOrdinal currLID = directoryMap_->getLocalElement(*ptr++); // Convert incoming GID to Directory LID
00350         TEST_FOR_EXCEPTION(currLID == LINVALID, std::logic_error,
00351             Teuchos::typeName(*this) << "::generateDirectory(): logic error. Please notify the Tpetra team.");
00352         nodeIDs_[currLID] = *ptr++;
00353         LIDs_[currLID] = *ptr++;
00354       }
00355     }
00356   } // end generateDirectory()
00357     
00358 } // namespace Tpetra
00359 
00360 //
00361 // Explicit instantiation macro
00362 //
00363 // Must be expanded from within the Tpetra namespace!
00364 //
00365 
00366 #define TPETRA_DIRECTORY_INSTANT(LO,GO,NODE) \
00367   \
00368   template class Directory< LO , GO , NODE >; \
00369 
00370 #endif // TPETRA_DIRECTORY_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines