Tpetra_Map.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_MAP_HPP
00030 #define TPETRA_MAP_HPP
00031 
00032 #include <Teuchos_OrdinalTraits.hpp>
00033 #include <Teuchos_CommHelpers.hpp>
00034 #include "Tpetra_MapDecl.hpp"
00035 #include "Tpetra_MapData.hpp"
00036 #include "Tpetra_Directory.hpp"
00037 
00038 //
00039 // compute isDistributed. it will be global.
00040 // min/max GID are always computed (from indexBase and numGlobal), and are therefore globally coherent as long as deps are.
00041 // indexBase and numGlobal must always be verified.
00042 // isContiguous is true for the "easy" constructors, assume false for the expert constructor
00043 // 
00044 // so we explicitly check    : isCont, numGlobal, indexBase
00045 // then we explicitly compute: MinAllGID, MaxAllGID
00046 // Data explicitly checks    : isDistributed
00047 
00048 namespace Tpetra {
00049 
00050   template<typename Ordinal>
00051   Map<Ordinal>::Map(Ordinal numGlobalEntries, Ordinal indexBase, const Platform<Ordinal> &platform, bool local) 
00052     : Teuchos::Object(local == false ? "Tpetra::Map" : "Teptra::Map(local)")
00053     , MapData_()
00054   {
00055     // distribute the entries across the nodes so that they are 
00056     // - non-overlapping
00057     // - contiguous
00058     // - as evenly distributed as possible
00059     const Ordinal one = Teuchos::OrdinalTraits<Ordinal>::one();
00060     const Ordinal zero = Teuchos::OrdinalTraits<Ordinal>::zero();
00061     const Ordinal negOne = zero - one;
00062 
00063     std::string errPrefix;
00064     errPrefix = "Tpetra::Map<" + Teuchos::OrdinalTraits<Ordinal>::name() 
00065               + ">::constructor(numGlobal,indexBase,platform,local): ";
00066 
00067     // get a internodal communicator from the Platform
00068     Teuchos::RCP< Teuchos::Comm<Ordinal> > comm = platform.createComm();
00069     if (local == false) 
00070     {
00071       Ordinal numImages = comm->getSize();
00072       Ordinal myImageID = comm->getRank();
00073 
00074       // check that numGlobalEntries,indexBase is equivalent across images
00075       std::vector<Ordinal> root_entries(2);
00076       root_entries[0] = numGlobalEntries;
00077       root_entries[1] = indexBase;
00078       Teuchos::broadcast(*comm,(Ordinal)0,(Ordinal)2,&root_entries[0]);   // broadcast 2 ordinals from node 0
00079       std::vector<Ordinal> localChecks(2);
00080       localChecks[0] = negOne;  // fail or pass
00081       localChecks[1] = zero;    // fail reason
00082       if (numGlobalEntries != root_entries[0]) {
00083         localChecks[0] = myImageID;
00084         localChecks[1] = one;
00085       }
00086       else if (indexBase != root_entries[1]) {
00087         localChecks[0] = myImageID;
00088         localChecks[1] = one + one;
00089       }
00090       // if checkPassed == negOne, then it passed on this proc
00091       // otherwise, checkPassed == myImageID signifies it did not pass
00092       // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass, as well as the reason
00093       // this will be negOne and zero if all procs passed
00094       Ordinal globalChecks[2];
00095       Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,(Ordinal)2,&localChecks[0],&globalChecks[0]);
00096       if (globalChecks[0] != negOne) {
00097         if (globalChecks[1] == one) {
00098           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00099               errPrefix << "numGlobal must be the same on all nodes (examine node " << globalChecks[0] << ").");
00100         }
00101         else if (globalChecks[1] == one+one) {
00102           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00103               errPrefix << "indexBase must be the same on all nodes (examine node " << globalChecks[0] << ").");
00104         }
00105         else {
00106           // logic error on my part
00107           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00108               errPrefix << "logic error. Please contact the Tpetra team.");
00109         }
00110       }
00111       // numgGlobalEntries is coherent, but is it valid?
00112       TEST_FOR_EXCEPTION(numGlobalEntries < zero, std::invalid_argument,
00113           errPrefix << "numGlobal == " << numGlobalEntries << ". Must be >= 0.");
00114 
00115       /* compute numLocalEntries
00116          We can write numGlobalEntries as follows:
00117          numGlobalEntries == numImages * B + remainder
00118          Each image is allocated entries as follows:
00119          [ B+1    iff myImageID <  remainder
00120          numLocalEntries = [
00121          [ B      iff myImageID >= remainder
00122          In the case that remainder == 0, then all images fall into the 
00123          latter case: numLocalEntries == B == numGlobalEntries / numImages
00124          It can then be shown that 
00125          numImages
00126          \Sum      numLocalEntries_i  == numGlobalEntries
00127          i=0
00128          This strategy is simple, requires no communication, and is optimal vis-a-vis
00129          uniform distribution of entries.
00130          This strategy is valid for any value of numGlobalEntries and numImages, 
00131          including the following border cases:
00132          - numImages == 1         -> remainder == 0 && numGlobalEntries == numLocalEntries
00133          - numEntries < numImages -> remainder == numGlobalEntries && numLocalEntries \in [0,1]
00134        */
00135       Ordinal numLocalEntries = numGlobalEntries / numImages;    
00136       Ordinal remainder = numGlobalEntries % numImages;
00137 #ifdef TEUCHOS_DEBUG
00138       // the above code assumes truncation. is that safe?
00139       SHARED_TEST_FOR_EXCEPTION(numLocalEntries * numImages + remainder != numGlobalEntries,
00140           std::logic_error, "Tpetra::Map::constructor(numGlobal,indexBase,platform): Ordinal does not implement division with truncation."
00141           << " Please contact Tpetra team.",*comm);
00142 #endif
00143       Ordinal start_index;
00144       if (myImageID < remainder) {
00145         ++numLocalEntries;
00146         /* the myImageID images before were each allocated 
00147            numGlobalEntries/numImages+1
00148            ergo, my offset is:
00149            myImageID * (numGlobalEntries/numImages+1)
00150            == myImageID * numLocalEntries
00151          */
00152         start_index = myImageID * numLocalEntries;
00153       }
00154       else {
00155         /* a quantity (equal to remainder) of the images before me
00156            were each allocated 
00157            numGlobalEntries/numImages+1
00158            entries. a quantity (equal to myImageID-remainder) of the remaining 
00159            images before me were each allocated 
00160            numGlobalEntries/numImages
00161            entries. ergo, my offset is:
00162            remainder*(numGlobalEntries/numImages+1) + (myImageID-remainder)*numGlobalEntries/numImages
00163            == remainder*numLocalEntries + remainder + myImageID*numLocalEntries - remainder*numLocalEntries
00164            == myImageID*numLocalEntries + remainder
00165          */
00166         start_index = myImageID*numLocalEntries + remainder;
00167       }
00168 
00169       // create empty maps between local and global entries: let the MapData constructor fill them
00170       Teuchos::ArrayRCP<Ordinal> lgMap = Teuchos::null;
00171       std::map<Ordinal,Ordinal> glMap;
00172 
00173       // compute the min/max global IDs
00174       Ordinal minAllGID = indexBase;
00175       Ordinal maxAllGID = indexBase + numGlobalEntries - one;
00176       Ordinal minMyGID  = start_index + indexBase;
00177       Ordinal maxMyGID  = minMyGID + numLocalEntries - one;
00178 
00179       Teuchos::RCP< Platform<Ordinal> > platform_clone = platform.clone();
00180 
00181       // create a MapData structure 
00182       MapData_ = Teuchos::rcp( new MapData<Ordinal>(indexBase,numGlobalEntries,numLocalEntries,
00183             minAllGID,maxAllGID,minMyGID,maxMyGID,
00184             lgMap,glMap,true,  // we are contiguous
00185             platform_clone, comm) );
00186 
00187       // initialize the directory
00188       directorySetup();
00189     }
00190     else 
00191     {
00192       Teuchos::ArrayRCP<Ordinal> lgMap = Teuchos::null;
00193       std::map<Ordinal,Ordinal>  glMap;
00194 
00195       // compute the min/max global IDs
00196       Ordinal minAllGID = indexBase;
00197       Ordinal maxAllGID = indexBase + numGlobalEntries - one;
00198       Ordinal minMyGID  = minAllGID;
00199       Ordinal maxMyGID  = maxAllGID;
00200 
00201       Teuchos::RCP< Platform<Ordinal> > platform_clone = platform.clone();
00202 
00203       // create a MapData structure 
00204       MapData_ = Teuchos::rcp( new MapData<Ordinal>(indexBase,numGlobalEntries,numGlobalEntries,
00205             minAllGID,maxAllGID,minMyGID,maxMyGID,
00206             lgMap,glMap,true,  // we are contiguous
00207             platform_clone, comm,true) );
00208     }
00209   }
00210 
00211   template<typename Ordinal>
00212   Map<Ordinal>::Map (Ordinal numGlobalEntries, Ordinal numLocalEntries, Ordinal indexBase, 
00213                          const Platform<Ordinal> &platform)
00214     : Teuchos::Object("Tpetra::Map")
00215     , MapData_()
00216   {
00217     // Distribute the entries across the nodes so that they are 
00218     // - non-overlapping
00219     // - contiguous
00220     // This differs from Map(Ord,Ord,Plat) in that the user has specified the number of entries 
00221     // per node, so that they are not (necessarily) evenly distributed
00222 
00223     const Ordinal one = Teuchos::OrdinalTraits<Ordinal>::one();
00224     const Ordinal zero = Teuchos::OrdinalTraits<Ordinal>::zero();
00225     const Ordinal negOne = zero - one;
00226 
00227     std::string errPrefix;
00228     errPrefix = "Tpetra::Map<" + Teuchos::OrdinalTraits<Ordinal>::name() 
00229               + ">::constructor(numGlobal,numLocal,indexBase,platform): ";
00230 
00231     // get a internodal communicator from the Platform
00232     Teuchos::RCP< Teuchos::Comm<Ordinal> > comm = platform.createComm();
00233     Ordinal myImageID = comm->getRank();
00234     // for communicating failures 
00235     std::vector<Ordinal> localChecks(2), globalChecks(2);
00236 
00237     /* compute the global size 
00238        we are computing the number of global entries because exactly ONE of the following is true:
00239        - the user didn't specify it, and we need it
00240        - the user did specify it, but we need to 
00241          + validate it against the sum of the local sizes, and
00242          + ensure that it is the same on all nodes
00243      */
00244     Ordinal global_sum;
00245     Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,numLocalEntries,&global_sum);
00246     /* there are three errors we should be detecting:
00247        - numGlobalEntries != -1 and it is incorrect/invalid
00248        - numLocalEntries invalid (<0)
00249     */
00250     localChecks[0] = negOne;
00251     if (numLocalEntries < zero) {
00252       localChecks[0] = myImageID;
00253       localChecks[1] = one;
00254     }
00255     else if (numGlobalEntries < negOne) {
00256       localChecks[0] = myImageID;
00257       localChecks[1] = one+one;
00258     }
00259     else if (numGlobalEntries > negOne && numGlobalEntries != global_sum) {
00260       localChecks[0] = myImageID;
00261       localChecks[1] = one+one+one;
00262     }
00263     // now check that indexBase is equivalent across images
00264     Ordinal root_indexbase = indexBase;
00265     Teuchos::broadcast(*comm,0,&root_indexbase);   // broadcast one ordinal from node 0
00266     if (indexBase != root_indexbase) {
00267       localChecks[0] = myImageID;
00268       localChecks[1] = one+one+one+one;
00269     }
00270     // if checkPassed == negOne, then it passed on this proc
00271     // otherwise, checkPassed == myImageID signifies it did not pass
00272     // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass
00273     // this will be negOne if all procs passed
00274     Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,(Ordinal)2,&localChecks[0],&globalChecks[0]);
00275     if (globalChecks[0] != negOne) {
00276       if (globalChecks[1] == one) {
00277         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00278             errPrefix << "numLocal is invalid (< 0) on at least one node (possibly node " 
00279             << globalChecks[0] << ").");
00280       }
00281       else if (globalChecks[1] == one+one) {
00282         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00283             errPrefix << "numGlobal is invalid (< -1) on at least one node (possibly node " 
00284             << globalChecks[0] << ").");
00285       }
00286       else if (globalChecks[1] == one+one+one) {
00287         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00288             errPrefix << "numGlobal doesn't match sum of numLocal (== " 
00289             << global_sum << ") on at least one node (possibly node " 
00290             << globalChecks[0] << ").");
00291       }
00292       else if (globalChecks[1] == one+one+one+one) {
00293         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00294             errPrefix << "indexBase is not the same on all nodes (examine node " 
00295             << globalChecks[0] << ").");
00296       }
00297       else {
00298         // logic error on my part
00299         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00300             errPrefix << "logic error. Please contact the Tpetra team.");
00301       }
00302     }
00303     // set numGlobalEntries
00304     if (numGlobalEntries == negOne) {
00305       numGlobalEntries = global_sum;
00306     }
00307 
00308     // compute my local offset
00309     Ordinal start_index;
00310     Teuchos::scan(*comm,Teuchos::REDUCE_SUM,numLocalEntries,&start_index);
00311     start_index -= numLocalEntries;
00312 
00313     // create empty maps between local and global entries: let the MapData constructor fill them
00314     Teuchos::ArrayRCP<Ordinal> lgMap = Teuchos::null;
00315     std::map<Ordinal,Ordinal> glMap;
00316 
00317     // compute the min/max global IDs
00318     Ordinal minAllGID = indexBase;
00319     Ordinal maxAllGID = indexBase + numGlobalEntries - one;
00320     Ordinal minMyGID = start_index + indexBase;
00321     Ordinal maxMyGID = minMyGID + numLocalEntries - one;
00322 
00323     Teuchos::RCP< Platform<Ordinal> > platform_clone = platform.clone();
00324 
00325     // create a MapData structure 
00326     MapData_ = Teuchos::rcp( new MapData<Ordinal>(indexBase,numGlobalEntries,numLocalEntries,
00327                                                       minAllGID,maxAllGID,minMyGID,maxMyGID,
00328                                                       lgMap,glMap,true,  // we are contiguous
00329                                                       platform_clone, comm) );
00330 
00331     // initialize the directory
00332     directorySetup();
00333   }
00334 
00335   template<typename Ordinal>
00336   Map<Ordinal>::Map (Ordinal numGlobalEntries, const Teuchos::ArrayView<const Ordinal> &entryList, 
00337                      Ordinal indexBase, const Platform<Ordinal> &platform) 
00338     : Teuchos::Object("Tpetra::Map")
00339     , MapData_()
00340   {
00341     // Distribute the entries across the nodes in an arbitrary user-specified manner
00342     // They are not necessarily contiguous or evenly distributed
00343     const Ordinal one = Teuchos::OrdinalTraits<Ordinal>::one();
00344     const Ordinal zero = Teuchos::OrdinalTraits<Ordinal>::zero();
00345     const Ordinal negOne = zero - one;
00346 
00347     Ordinal numLocalEntries = entryList.size();
00348 
00349     std::string errPrefix;
00350     errPrefix = "Tpetra::Map<" + Teuchos::OrdinalTraits<Ordinal>::name() 
00351               + ">::constructor(numGlobal,entryList,indexBase,platform): ";
00352 
00353     // platform & comm setup
00354     Teuchos::RCP< Teuchos::Comm<Ordinal> > comm = platform.createComm();
00355     Ordinal myImageID = comm->getRank();
00356     // for communicating failures 
00357     std::vector<Ordinal> localChecks(2), globalChecks(2);
00358 
00359     /* compute the global size 
00360        we are computing the number of global entries because exactly ONE of the following is true:
00361        - the user didn't specify it, and we need it
00362        - the user did specify it, but we need to 
00363          + validate it against the sum of the local sizes, and
00364          + ensure that it is the same on all nodes
00365      */
00366     Ordinal global_sum;
00367     Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,numLocalEntries,&global_sum);
00368     localChecks[0] = negOne;
00369     if (numGlobalEntries < negOne) {
00370       localChecks[0] = myImageID;
00371       localChecks[1] = one;
00372     }
00373     else if (numGlobalEntries > negOne && numGlobalEntries != global_sum) {
00374       localChecks[0] = myImageID;
00375       localChecks[1] = one + one;
00376     }
00377     // now check that indexBase is equivalent across images
00378     Ordinal root_indexbase = indexBase;
00379     Teuchos::broadcast(*comm,0,&root_indexbase);   // broadcast one ordinal from node 0
00380     if (indexBase != root_indexbase) {
00381       localChecks[0] = myImageID;
00382       localChecks[1] = one + one + one;
00383     }
00384     // if checkPassed == negOne, then it passed on this proc
00385     // otherwise, checkPassed == myImageID signifies it did not pass
00386     // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass
00387     // this will be negOne if all procs passed
00388     Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,(Ordinal)2,&localChecks[0],&globalChecks[0]);
00389     if (globalChecks[0] != negOne) {
00390       if (globalChecks[1] == one) {
00391         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00392             errPrefix << "numGlobal is invalid (< -1) on at least one node (possibly node "
00393             << globalChecks[0] << ").");
00394       }
00395       else if (globalChecks[1] == one+one) {
00396         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00397             errPrefix << "numGlobal doesn't match sum of numLocal (" 
00398             << global_sum << ") on at least one node (possibly node "
00399             << globalChecks[0] << ").");
00400       }
00401       else if (globalChecks[1] == one+one+one) {
00402         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00403             errPrefix << "indexBase is not the same on all nodes (possibly node "
00404             << globalChecks[0] << ").");
00405       }
00406       else {
00407         // logic error on my part
00408         TEST_FOR_EXCEPTION(true,std::invalid_argument,
00409             errPrefix << "logic error. Please contact the Tpetra team.");
00410       }
00411     }
00412     // set numGlobalEntries
00413     if (numGlobalEntries == negOne) {
00414       numGlobalEntries = global_sum;
00415     }
00416 
00417     // setup lgmap and glmap, and min/maxMyGIDs
00418     std::map<Ordinal,Ordinal> glMap;
00419     // assume for now that there are numLocalEntries (there may be less, if some
00420     // GIDs are duplicated in entryList)
00421     Teuchos::ArrayRCP<Ordinal> lgMap;
00422     Ordinal minMyGID = indexBase;
00423     Ordinal maxMyGID = indexBase;
00424     // create the GID to LID map; do not assume GID in entryList are distinct.
00425     // in the case that a GID is duplicated, keep the previous LID
00426     // this is necessary so that LIDs are in [0,numLocal)
00427     // FINISH: make sure this is legal
00428     {
00429       Ordinal numLIDs = zero;
00430       if (numLocalEntries > zero) {
00431         lgMap = Teuchos::arcp<Ordinal>(numLocalEntries);
00432         for(Ordinal i = zero; i < numLocalEntries; i++) {
00433           if (glMap.find(entryList[i]) == glMap.end()) {
00434             lgMap[numLIDs] = entryList[i];   // lgMap: LID to GID
00435             glMap[entryList[i]] = numLIDs;   // glMap: GID to LID
00436             numLIDs++;
00437           }
00438         }
00439         minMyGID = *std::min_element(entryList.begin(), entryList.end());
00440         maxMyGID = *std::max_element(entryList.begin(), entryList.end());
00441         // shrink lgMap appropriately
00442         if (numLocalEntries != numLIDs) {
00443           numLocalEntries = numLIDs;
00444           lgMap = lgMap.persistingView(0,numLocalEntries);
00445         }
00446       }
00447     } 
00448 
00449     // set min/maxAllGIDs
00450     Ordinal minAllGID, maxAllGID;
00451     {
00452       Ordinal minmaxAllGIDlcl[2], minmaxAllGIDgbl[2];
00453       minmaxAllGIDlcl[0] = -minMyGID;  // negative allows us to do a single
00454       minmaxAllGIDlcl[1] =  maxMyGID;  // reduction below
00455       Teuchos::reduceAll<Ordinal>(*comm,Teuchos::REDUCE_MAX,2,minmaxAllGIDlcl,minmaxAllGIDgbl);
00456       minAllGID = -minmaxAllGIDgbl[0];
00457       maxAllGID =  minmaxAllGIDgbl[1];
00458     }
00459     TEST_FOR_EXCEPTION(minAllGID < indexBase, std::invalid_argument,
00460         errPrefix << "minimum GID == " << minAllGID << " is less than indexBase == " << indexBase);
00461 
00462     Teuchos::RCP< Platform<Ordinal> > platform_clone = platform.clone();
00463 
00464     // call ESData constructor
00465     MapData_ = Teuchos::rcp(new MapData<Ordinal>(indexBase, numGlobalEntries, 
00466           numLocalEntries, minAllGID, maxAllGID, 
00467           minMyGID, maxMyGID, lgMap, glMap, 
00468           false, platform_clone, comm));
00469 
00470     // initialize directory
00471     directorySetup();
00472   }
00473 
00474   template<typename Ordinal>
00475   Map<Ordinal>::Map (const Map<Ordinal> &map) 
00476     : Teuchos::Object(map.label())
00477     , MapData_(map.MapData_)
00478   {}
00479 
00480   template<typename Ordinal>
00481   Map<Ordinal>::~Map () 
00482   { }
00483 
00484   template<typename Ordinal>
00485   Ordinal Map<Ordinal>::getNumGlobalEntries() const {
00486     return MapData_->numGlobalEntries_;
00487   }
00488 
00489   template<typename Ordinal>
00490   Ordinal Map<Ordinal>::getNumMyEntries() const {
00491     return MapData_->numMyEntries_;
00492   }
00493 
00494   template<typename Ordinal>
00495   Ordinal Map<Ordinal>::getIndexBase() const {
00496     return MapData_->indexBase_;
00497   }
00498 
00499   template<typename Ordinal>
00500   Ordinal Map<Ordinal>::getMinLocalIndex() const {
00501     return Teuchos::OrdinalTraits<Ordinal>::zero();
00502   }
00503 
00504   template<typename Ordinal>
00505   Ordinal Map<Ordinal>::getMaxLocalIndex() const {
00506     return MapData_->numMyEntries_-1;
00507   }
00508 
00509   template<typename Ordinal>
00510   Ordinal Map<Ordinal>::getMinGlobalIndex() const {
00511     return MapData_->minMyGID_;
00512   }
00513 
00514   template<typename Ordinal>
00515   Ordinal Map<Ordinal>::getMaxGlobalIndex() const {
00516     return MapData_->maxMyGID_;
00517   }
00518 
00519   template<typename Ordinal>
00520   Ordinal Map<Ordinal>::getMinAllGlobalIndex() const {
00521     return MapData_->minAllGID_;
00522   }
00523 
00524   template<typename Ordinal>
00525   Ordinal Map<Ordinal>::getMaxAllGlobalIndex() const {
00526     return MapData_->maxAllGID_;
00527   }
00528 
00529   template<typename Ordinal>
00530   Ordinal Map<Ordinal>::getLocalIndex(Ordinal globalIndex) const {
00531     if (MapData_->contiguous_) {
00532       TEST_FOR_EXCEPTION(
00533         globalIndex < getMinGlobalIndex() || globalIndex > getMaxGlobalIndex(), 
00534         std::invalid_argument,
00535         "Tpetra::Map<" << Teuchos::OrdinalTraits<Ordinal>::name() 
00536                        << ">::getLocalIndex(gid): gid does not belong to the map."
00537       );
00538       return globalIndex - getMinGlobalIndex();
00539     }
00540     else {
00541       typename std::map<Ordinal,Ordinal>::const_iterator i;
00542       i = MapData_->glMap_.find(globalIndex);
00543       TEST_FOR_EXCEPTION(
00544         i == MapData_->glMap_.end(), std::invalid_argument,
00545         "Tpetra::Map<" << Teuchos::OrdinalTraits<Ordinal>::name() 
00546                        << ">::getLocalIndex(gid): gid does not belong to the map."
00547       );
00548       return i->second;
00549     }
00550   }
00551 
00552   template<typename Ordinal>
00553   Ordinal Map<Ordinal>::getGlobalIndex(Ordinal localIndex) const {
00554     TEST_FOR_EXCEPTION(
00555         localIndex < 0 || localIndex > MapData_->numMyEntries_-1,
00556         std::invalid_argument,
00557         "Tpetra::Map<" << Teuchos::OrdinalTraits<Ordinal>::name() 
00558            << ">::getGlobalIndex(lid): lid not valid.");
00559     if (MapData_->contiguous_) {
00560       return localIndex + getMinGlobalIndex();
00561     }
00562     else {
00563       return MapData_->lgMap_[localIndex];
00564     }
00565   }
00566 
00567   template<typename Ordinal>
00568   bool Map<Ordinal>::isMyLocalIndex(Ordinal localIndex) const {
00569     // we have local indices in [0,numLocalEntries)
00570     if (localIndex >= 0 && localIndex < MapData_->numMyEntries_) {
00571       return true;
00572     }
00573     return false;
00574   }
00575 
00576   template<typename Ordinal>
00577   bool Map<Ordinal>::isMyGlobalIndex(Ordinal globalIndex) const {
00578     if (MapData_->contiguous_) {
00579       return (MapData_->minMyGID_ <= globalIndex) && (globalIndex <= MapData_->maxMyGID_);
00580     }
00581     else {
00582       typename std::map<Ordinal,Ordinal>::iterator i;
00583       i = MapData_->glMap_.find(globalIndex);
00584       return (i != MapData_->glMap_.end());
00585     }
00586   }
00587 
00588   template<typename Ordinal>
00589   bool Map<Ordinal>::isContiguous() const {
00590     return MapData_->contiguous_;
00591   }
00592 
00593   template<typename Ordinal>
00594   bool Map<Ordinal>::isCompatible (const Map< Ordinal> &map) const {
00595     // this is wrong: the non-compat procs will return false, then the compat 
00596     // procs will stall on this communication. this is demonstrated by tpetra/test/Map/Map_test.exe
00597     /*
00598     if (getNumGlobalEntries() != map.getNumGlobalEntries() ||
00599         getNumMyEntries() != map.getNumMyEntries()) 
00600     { return false; }
00601     */
00602 
00603     // check to make sure distribution is the same
00604     char iscompat_lcl;
00605     if (getNumGlobalEntries() != map.getNumGlobalEntries() ||
00606         getNumMyEntries() != map.getNumMyEntries()) 
00607     {
00608       // NOT compat on this node
00609       iscompat_lcl = Teuchos::ScalarTraits<char>::zero();
00610     }
00611     else {
00612       // compat on this node
00613       iscompat_lcl = Teuchos::ScalarTraits<char>::one();
00614     }
00615     char iscompat_gbl;
00616     Teuchos::reduceAll(*MapData_->comm_,Teuchos::REDUCE_MIN,iscompat_lcl,&iscompat_gbl);
00617     return(iscompat_gbl == Teuchos::ScalarTraits<char>::one());
00618   }
00619 
00620   template<typename Ordinal>
00621   bool Map<Ordinal>::isSameAs (const Map<Ordinal> &map) const {
00622     if (MapData_ == map.MapData_) {
00623       // we should assume that this is globally coherent
00624       // if they share the same underlying MapData, then they must be equivalent maps
00625       return true;
00626     }
00627 
00628     // check all other globally coherent properties
00629     // if they do not share each of these properties, then they cannot be 
00630     // equivalent maps
00631     if ( (getMinGlobalIndex()   != map.getMinGlobalIndex()) ||
00632          (getMaxGlobalIndex()   != map.getMaxGlobalIndex()) ||
00633          (getNumGlobalEntries() != map.getNumGlobalEntries()) ||
00634          (isDistributed()       != map.isDistributed()) ) {
00635       return false;
00636     }
00637 
00638     // If we get this far, we need to check local properties and the 
00639     // communicate same-ness across all nodes
00640     // we prefer local work over communication, ergo, we will perform all
00641     // comparisons and conduct a single communication
00642     char isSame_lcl = Teuchos::ScalarTraits<char>::one();
00643 
00644     // check number of entries owned by this node
00645     if (getNumMyEntries() != map.getNumMyEntries()) {
00646       isSame_lcl = Teuchos::ScalarTraits<char>::zero();
00647     }
00648 
00649     // check the identity of the entries owned by this node
00650     // only do this if we haven't already determined not-same-ness
00651     if (isSame_lcl == Teuchos::ScalarTraits<char>::one()) {
00652       // if they are contiguous, we can check the ranges easily
00653       // if they are not contiguous, we must check the individual LID -> GID mappings
00654       // the latter approach is valid in either case, but the former is faster
00655       if (MapData_->contiguous_ == true && map.MapData_->contiguous_ == true) {
00656         if (MapData_->minMyGID_ != map.MapData_->minMyGID_ ||
00657             MapData_->maxMyGID_ != map.MapData_->maxMyGID_) {
00658           isSame_lcl = Teuchos::ScalarTraits<char>::zero();
00659         }
00660       }
00661       else {
00662         if (!std::equal(    MapData_->lgMap_.begin(),MapData_->lgMap_.end(),
00663                         map.MapData_->lgMap_.begin())) {
00664           isSame_lcl = Teuchos::ScalarTraits<char>::zero();
00665         }
00666       }
00667     }
00668 
00669     // now, determine if we detected not-same-ness on any node
00670     char isSame_gbl;
00671     Teuchos::reduceAll(*MapData_->comm_,Teuchos::REDUCE_MIN,isSame_lcl,&isSame_gbl);
00672     return(isSame_gbl == Teuchos::ScalarTraits<char>::one());
00673   }
00674 
00675   template<typename Ordinal>
00676   bool Map<Ordinal>::operator== (const Map< Ordinal > &map) const {
00677     return isSameAs(map);
00678   }
00679 
00680   template<typename Ordinal>
00681   bool Map<Ordinal>::operator!= (const Map< Ordinal > &map) const {
00682     return !isSameAs(map);
00683   }
00684 
00685   template<typename Ordinal>
00686   Map<Ordinal>& Map<Ordinal>::operator = (const Map<Ordinal> & Source) {
00687     MapData_ = Source.MapData_;
00688     return *this;
00689   }
00690 
00691   template<typename Ordinal>
00692   Teuchos::ArrayView<const Ordinal>
00693   Map<Ordinal>::getMyGlobalEntries() const {
00694     const Ordinal ZERO = Teuchos::OrdinalTraits<Ordinal>::zero();
00695     // check to see if lgMap is empty
00696     // if so (and we have local entries), then fill it.
00697     if (MapData_->lgMap_ == Teuchos::null && MapData_->numMyEntries_ > ZERO) {
00698       // this would have been set up for a non-contiguous map
00699 #ifdef TEUCHOS_DEBUG
00700       TEST_FOR_EXCEPTION(MapData_->contiguous_ != true, std::logic_error,
00701           "Tpetra::Map::getMyGlobalEntries: logic error. Please notify the Tpetra team.");
00702 #endif
00703       Teuchos::ArrayRCP<Ordinal> lgMap = Teuchos::arcp<Ordinal>(MapData_->numMyEntries_);
00704       MapData_->lgMap_ = lgMap;
00705       for (Ordinal lid=MapData_->minMyGID_; lid <= MapData_->maxMyGID_; ++lid) {
00706         *(lgMap++) = lid;
00707       }
00708     }
00709     return MapData_->lgMap_();
00710   }
00711 
00712   template<typename Ordinal>
00713   bool Map<Ordinal>::isDistributed() const {
00714     return MapData_->distributed_;
00715   }
00716 
00717   template<typename Ordinal>
00718   void Map<Ordinal>::print(std::ostream& os) const {
00719     const Ordinal zero = Teuchos::OrdinalTraits<Ordinal>::zero();
00720     const Ordinal nME = getNumMyEntries();
00721     
00722     using std::endl;
00723 
00724     getMyGlobalEntries(); // neglect the output, we call this here only to ensure that sure list is generated
00725     int myImageID = MapData_->comm_->getRank();
00726     int numImages = MapData_->comm_->getSize();
00727 
00728     for(int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00729       if(myImageID == imageCtr) {
00730         if(myImageID == 0) { // this is the root node (only output this info once)
00731           os << endl 
00732              << "Number of Global Entries = " << getNumGlobalEntries()  << endl
00733              << "Maximum of all GIDs      = " << getMaxAllGlobalIndex() << endl
00734              << "Minimum of all GIDs      = " << getMinAllGlobalIndex() << endl
00735              << "Index Base               = " << getIndexBase()         << endl;
00736         }
00737         os << endl;
00738         os << "Number of Local Elements   = " << nME           << endl
00739            << "Maximum of my GIDs         = " << getMaxGlobalIndex() << endl
00740            << "Minimum of my GIDs         = " << getMinGlobalIndex() << endl;
00741         os << endl;
00742         os << std::setw(16) << "ImageID"
00743            << std::setw(16) << "Local Index"
00744            << std::setw(16) << "Global Index"
00745            << endl;
00746         Teuchos::ArrayView<const Ordinal> myEntries = getMyGlobalEntries();
00747         for(Ordinal i = zero; i < nME; i++) {
00748           os << std::setw(16) << myImageID 
00749              << std::setw(16) << i
00750              << std::setw(16) << myEntries[i]
00751              << endl;
00752         }
00753         os << std::flush;
00754       }
00755       // Do a few global ops to give I/O a chance to complete
00756       MapData_->comm_->barrier();
00757       MapData_->comm_->barrier();
00758       MapData_->comm_->barrier();
00759     }
00760   }
00761 
00762   template<typename Ordinal>
00763   void Map<Ordinal>::directorySetup() {
00764     if (getNumGlobalEntries() != Teuchos::OrdinalTraits<Ordinal>::zero()) {
00765       if (MapData_->directory_ == Teuchos::null) {
00766         MapData_->directory_ = Teuchos::rcp( new Directory<Ordinal>(*this) );
00767       }
00768     }
00769   }
00770 
00771   template<typename Ordinal>
00772   bool Map<Ordinal>::getRemoteIndexList(
00773       const Teuchos::ArrayView<const Ordinal> & GIDList, 
00774       const Teuchos::ArrayView<Ordinal> & imageIDList, 
00775       const Teuchos::ArrayView<Ordinal> & LIDList) const 
00776   {
00777     if (GIDList.size() == 0) return false;
00778     TEST_FOR_EXCEPTION(getNumGlobalEntries() == 0, std::runtime_error,
00779         "Tpetra::Map<" + Teuchos::OrdinalTraits<Ordinal>::name() 
00780         + ">::getRemoteIndexList(): getRemoteIndexList() cannot be called, zero entries on node.");
00781     return MapData_->directory_->getDirectoryEntries(GIDList, imageIDList, LIDList);
00782   }
00783 
00784   template<typename Ordinal>
00785   bool Map<Ordinal>::getRemoteIndexList(
00786       const Teuchos::ArrayView<const Ordinal> & GIDList, 
00787       const Teuchos::ArrayView<Ordinal> & imageIDList) const 
00788   {
00789     if (GIDList.size() == 0) return false;
00790     TEST_FOR_EXCEPTION(getNumGlobalEntries() == 0, std::runtime_error,
00791         "Tpetra::Map<" + Teuchos::OrdinalTraits<Ordinal>::name() 
00792         + ">::getRemoteIndexList(): getRemoteIndexList() cannot be called, zero entries on node.");
00793     return MapData_->directory_->getDirectoryEntries(GIDList, imageIDList);
00794   }
00795 
00796 
00797   template<typename Ordinal>
00798   Teuchos::RCP<const Platform<Ordinal> >
00799   Map<Ordinal>::getPlatform() const {
00800     return MapData_->platform_;
00801   }
00802 
00803   template<typename Ordinal>
00804   Teuchos::RCP<const Teuchos::Comm<Ordinal> >
00805   Map<Ordinal>::getComm() const {
00806     return MapData_->comm_;
00807   }
00808 
00809 
00810 } // Tpetra namespace
00811 
00812 #endif // TPETRA_MAP_HPP
00813 

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