Tpetra_Directory.hpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2004) 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_OrdinalTraits.hpp>
00033 #include <Teuchos_CommHelpers.hpp>
00034 #include <Teuchos_as.hpp>
00035 #include "Tpetra_DirectoryDecl.hpp"
00036 #include "Tpetra_Distributor.hpp"
00037 #include "Tpetra_Map.hpp"
00038 
00039 namespace Tpetra {
00040 
00041   template<typename Ordinal>
00042   Directory<Ordinal>::Directory(const Map<Ordinal> & map)  
00043     : Teuchos::Object("Tpetra::Directory") 
00044     , map_(map) 
00045   {
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         const Ordinal one = Teuchos::OrdinalTraits<Ordinal>::one();
00055         // make room for the min on each proc, plus one entry at the end for the max cap
00056         allMinGIDs_.resize(comm_->getSize() + one);
00057         // get my min
00058         Ordinal minMyGID = map.getMinGlobalIndex();
00059         // gather all of the mins into the first getSize() entries of allMinDIGs_
00060         Teuchos::gatherAll(*comm_,one,&minMyGID,(Ordinal)(map.getComm()->getSize()),&allMinGIDs_.front());
00061         // put the max cap at the end
00062         allMinGIDs_.back() = map.getMaxAllGlobalIndex() + one; // FINISH: is this right?
00063       }
00064       // Otherwise we have to generate the directory using MPI calls.
00065       else {
00066         generateDirectory();
00067       }
00068     }
00069   }
00070 
00071   template<typename Ordinal>
00072   Directory<Ordinal>::Directory(const Directory<Ordinal> & directory)
00073     : Teuchos::Object(directory.label()) 
00074     , map_(directory.map_) 
00075     , comm_(directory.comm_)
00076     , allMinGIDs_(directory.allMinGIDs_)
00077     , imageIDs_(directory.imageIDs_)
00078     , LIDs_(directory.LIDs_)
00079   {}
00080 
00081   template<typename Ordinal>
00082   Directory<Ordinal>::~Directory() {}
00083 
00084   template<typename Ordinal>
00085   bool Directory<Ordinal>::getDirectoryEntries(
00086       const Teuchos::ArrayView<const Ordinal> &globalEntries, 
00087       const Teuchos::ArrayView<Ordinal> &images) const 
00088   {
00089     return getEntries(globalEntries, images, Teuchos::ArrayView<Ordinal>(Teuchos::null), false);
00090   }
00091 
00092   template<typename Ordinal>
00093   bool Directory<Ordinal>::getDirectoryEntries(
00094       const Teuchos::ArrayView<const Ordinal> &globalEntries, 
00095       const Teuchos::ArrayView<Ordinal> &images, 
00096       const Teuchos::ArrayView<Ordinal> &localEntries) const 
00097   {
00098     return getEntries(globalEntries, images, localEntries, true);
00099   }
00100 
00101   template<typename Ordinal>
00102   bool Directory<Ordinal>::getEntries(
00103       const Teuchos::ArrayView<const Ordinal> &globalEntries, 
00104       const Teuchos::ArrayView<Ordinal> &images, 
00105       const Teuchos::ArrayView<Ordinal> &localEntries, 
00106             bool computeLIDs) const 
00107   {
00108     const Ordinal ZERO = Teuchos::OrdinalTraits<Ordinal>::zero();
00109     const Ordinal ONE  = Teuchos::OrdinalTraits<Ordinal>::one();
00110     const Ordinal NEGONE = ZERO - ONE;
00111 
00112     bool invalidGIDs = false;
00113 
00114     // fill images and localEntries with -1s
00115     TEST_FOR_EXCEPTION(images.size() != globalEntries.size(), std::runtime_error,
00116         "Tpetra::Directory<" << Teuchos::TypeNameTraits<Ordinal>::name() << 
00117         ">::getEntries(): Output buffers are not allocated properly.");
00118     std::fill(images.begin(),images.end(),NEGONE);
00119     if (computeLIDs) {
00120       TEST_FOR_EXCEPTION(localEntries.size() != globalEntries.size(), std::runtime_error,
00121           "Tpetra::Directory<" << Teuchos::TypeNameTraits<Ordinal>::name() << 
00122           ">::getEntries(): Output buffers are not allocated properly.");
00123       std::fill(localEntries.begin(),localEntries.end(),NEGONE);
00124     }
00125 
00126     const Ordinal numImages  = comm_->getSize();
00127     const Ordinal myImageID  = comm_->getRank();
00128     const Ordinal numEntries = globalEntries.size();
00129     const Ordinal nOverP     = map_.getNumGlobalEntries() / numImages;
00130 
00131     if (map_.isDistributed() == false) {
00132       // Easiest case: Map is serial or locally-replicated
00133       typename Teuchos::ArrayView<Ordinal>::iterator imgptr = images.begin(),    
00134                                                      lidptr = localEntries.begin();
00135       for (typename Teuchos::ArrayView<const Ordinal>::iterator gid = globalEntries.begin(); gid != globalEntries.end(); ++gid) {
00136         if (map_.isMyGlobalIndex(*gid)) {
00137           *imgptr++ = myImageID;
00138           if (computeLIDs) {
00139             *lidptr++ = map_.getLocalIndex(*gid);
00140           }
00141         }
00142         else {
00143           imgptr++;
00144           if (computeLIDs) {
00145             lidptr++;
00146           }
00147           invalidGIDs = true;
00148         }
00149       }
00150     }
00151     else if (map_.isContiguous()) {
00152       // Next Easiest Case: Map is distributed but allocated contiguously
00153       typename Teuchos::ArrayView<Ordinal>::iterator imgptr = images.begin(),    
00154                                                      lidptr = localEntries.begin();
00155       for (typename Teuchos::ArrayView<const Ordinal>::iterator gid = globalEntries.begin(); gid != globalEntries.end(); ++gid) {
00156         Ordinal LID = NEGONE; // Assume not found
00157         Ordinal image = NEGONE;
00158         Ordinal GID = *gid;
00159         // Guess uniform distribution and start a little above it
00160         Ordinal image1 = TEUCHOS_MIN((GID / TEUCHOS_MAX(nOverP, ONE)) + Teuchos::as<Ordinal>(2), numImages - ONE);
00161         bool found = false;
00162         while (image1 >= ZERO && image1 < numImages) {
00163           if (allMinGIDs_[image1] <= GID) {
00164             if (GID < allMinGIDs_[image1 + ONE]) {
00165               found = true;
00166               break;
00167             }
00168             else {
00169               image1++;
00170             }
00171           }
00172           else {
00173             image1--;
00174           }
00175         }
00176         if (found) {
00177           image = image1;
00178           LID = GID - allMinGIDs_[image];
00179         }
00180         else {
00181           invalidGIDs = true;
00182         }
00183         *imgptr++ = image;
00184         if (computeLIDs) {
00185           *lidptr++ = LID;
00186         }
00187       }
00188     }
00189     else {
00190       // General Case: Map is distributed and allocated arbitrarily
00191       // Here we need to set up an actual directory structure
00192       Ordinal packetSize = Teuchos::as<Ordinal>(2);
00193       if (computeLIDs) {
00194         ++packetSize;
00195       }
00196 
00197       Distributor<Ordinal> distor(comm_);
00198 
00199       // Get directory locations for the requested list of entries
00200       Teuchos::Array<Ordinal> dirImages(numEntries);
00201       invalidGIDs = directoryMap_->getRemoteIndexList(globalEntries, dirImages());
00202       // Check for unfound globalEntries and set corresponding images to -1
00203       Ordinal numMissing = ZERO;
00204       if (invalidGIDs) {
00205         for (Ordinal i = ZERO; i < numEntries; ++i) {
00206           if (dirImages[i] == NEGONE) {
00207             images[i] = NEGONE;
00208             if (computeLIDs) {
00209               localEntries[i] = NEGONE;
00210             }
00211             numMissing++;
00212           }
00213         }
00214       }
00215 
00216       Teuchos::ArrayRCP<Ordinal> sendGIDs, sendImages;
00217       distor.createFromRecvs(globalEntries, dirImages(), sendGIDs, sendImages);
00218       Ordinal numSends = Teuchos::as<Ordinal>(sendGIDs.size());
00219 
00220       Ordinal currLID;
00221       Teuchos::Array<Ordinal> exports(packetSize*numSends);
00222       {
00223         typename Teuchos::Array<Ordinal>::iterator ptr = exports.begin();
00224         for(Ordinal i = ZERO; i < numSends; i++) {
00225           Ordinal currGID = sendGIDs[i];
00226           *ptr++ = currGID;
00227           currLID = directoryMap_->getLocalIndex(currGID);
00228           assert(currLID != NEGONE); // Internal error
00229           *ptr++ = imageIDs_[currLID];
00230           if(computeLIDs) {
00231             *ptr++ = LIDs_[currLID];
00232           }
00233         }
00234       }
00235 
00236       Teuchos::Array<Ordinal> imports(packetSize*distor.getTotalReceiveLength());
00237       distor.doPostsAndWaits(exports().getConst(), packetSize, imports());
00238 
00239       typename Teuchos::Array<Ordinal>::iterator ptr = imports.begin();
00240       const Ordinal numRecv = numEntries - numMissing;
00241       for(Ordinal i = ZERO; i < numRecv; i++) {
00242         currLID = *ptr++;
00243         for(Ordinal j = ZERO; j < numEntries; j++) {
00244           if(currLID == globalEntries[j]) {
00245             images[j] = *ptr++;
00246             if(computeLIDs) {
00247               localEntries[j] = *ptr++;
00248             }
00249             break;
00250           }
00251         }
00252       }
00253     }
00254     return invalidGIDs;
00255   }
00256 
00257 
00258   // directory setup for non-contiguous Map
00259   template<typename Ordinal>
00260   void Directory<Ordinal>::generateDirectory() 
00261   {
00262     const Ordinal ONE = Teuchos::OrdinalTraits<Ordinal>::one();
00263     const Ordinal ZERO = Teuchos::OrdinalTraits<Ordinal>::zero();
00264     const Ordinal NEGONE = ZERO - ONE;
00265           
00266     const Ordinal minAllGID = map_.getMinAllGlobalIndex();
00267     const Ordinal maxAllGID = map_.getMaxAllGlobalIndex();
00268 
00269     comm_ = map_.getComm();
00270 
00271     // DirectoryMap will have a range of elements from the minimum to the maximum
00272     // GID of the user Map, and an indexBase of minAllGID from the user Map
00273     Ordinal numGlobalEntries = maxAllGID - minAllGID + ONE;
00274 
00275     // Create a uniform linear map to contain the directory
00276     directoryMap_ = Teuchos::rcp(new Map<Ordinal>(numGlobalEntries, minAllGID, *map_.getPlatform()));
00277 
00278     Ordinal dir_numMyEntries = directoryMap_->getNumMyEntries();
00279 
00280     // Allocate imageID List and LID List.  Initialize to -1s.
00281     // Initialize values to -1 in case the user global element list does
00282     // fill all IDs from minAllGID to maxAllGID (e.g., allows global indices to be 
00283     // all even integers).
00284     imageIDs_.resize(dir_numMyEntries, NEGONE);
00285     LIDs_.resize(dir_numMyEntries, NEGONE);
00286 
00287     // Get list of images owning the directory entries for the Map GIDs
00288     Ordinal myImageID = comm_->getRank();
00289     Ordinal numMyEntries = map_.getNumMyEntries();
00290     std::vector<Ordinal> sendImageIDs(numMyEntries);
00291     Teuchos::ArrayView<const Ordinal> myGlobalEntries = map_.getMyGlobalEntries();
00292     // a "true" return here indicates that one of myGlobalEntries (from map_) is not on the map directoryMap_, indicating that 
00293     // it lies outside of the range [minAllGID,maxAllGID] (from map_). this means something is wrong with map_.
00294     TEST_FOR_EXCEPTION( directoryMap_->getRemoteIndexList(myGlobalEntries, sendImageIDs) == true, std::logic_error,
00295         "Tpetra::Directory::generateDirectory(): logic error. Please contact Tpetra team.");
00296 
00297     // Create distributor & call createFromSends
00298     Ordinal numReceives = ZERO;
00299     Distributor<Ordinal> distor(comm_);      
00300     distor.createFromSends(sendImageIDs, numReceives);
00301 
00302     // Execute distributor plan
00303     // Transfer GIDs, ImageIDs, and LIDs that we own to all images
00304     // End result is all images have list of all GIDs and corresponding ImageIDs and LIDs
00305     Ordinal packetSize = Teuchos::as<Ordinal>(3); // We will send GIDs, ImageIDs, and LIDs.
00306     Teuchos::Array<Ordinal> exportEntries(packetSize*numMyEntries);
00307     {
00308       typename Teuchos::Array<Ordinal>::iterator ptr = exportEntries.begin();
00309       for(Ordinal i = ZERO; i < numMyEntries; ++i) {
00310         *ptr++ = myGlobalEntries[i];
00311         *ptr++ = myImageID;
00312         *ptr++ = i;
00313       }
00314     }
00315 
00316     Teuchos::Array<Ordinal> importElements(packetSize*distor.getTotalReceiveLength());
00317     distor.doPostsAndWaits(exportEntries().getConst(), packetSize, importElements());
00318 
00319     {
00320       typename Teuchos::Array<Ordinal>::iterator ptr = importElements.begin();
00321       for(Ordinal i = ZERO; i < numReceives; i++) {
00322         Ordinal currLID = directoryMap_->getLocalIndex(*ptr++); // Convert incoming GID to Directory LID
00323         assert(currLID != NEGONE); // Internal error
00324         TEST_FOR_EXCEPTION(currLID == NEGONE, std::logic_error,
00325             "Tpetra::Directory<" << Teuchos::OrdinalTraits<Ordinal>::name() << ">::generateDirectory(): logic error. Please notify the Tpetra team.");
00326         imageIDs_[currLID] = *ptr++;
00327         LIDs_[currLID] = *ptr++;
00328       }
00329     }
00330   }
00331     
00332 } // namespace Tpetra
00333 
00334 #endif // TPETRA_DIRECTORY_HPP

Generated on Wed May 12 21:59:41 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7