Tpetra_Map_def.hpp

Go to the documentation of this file.
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 // TODO: make sure that Ordinal values in constructors aren't invalid()
00030 
00031 #ifndef TPETRA_MAP_DEF_HPP
00032 #define TPETRA_MAP_DEF_HPP
00033 
00034 #include <Teuchos_TypeNameTraits.hpp>
00035 #include <Teuchos_OrdinalTraits.hpp>
00036 #include <Teuchos_CommHelpers.hpp>
00037 #include <Teuchos_as.hpp>
00038 #include "Tpetra_Directory.hpp" // must include for implicit instantiation to work
00039 #include "Tpetra_Util.hpp"
00040 #include <stdexcept>
00041 
00042 #ifdef DOXYGEN_USE_ONLY
00043   #includee "Tpetra_Map_decl.hpp"
00044 #endif
00045 
00051 //
00052 // compute isDistributed. it will be global.
00053 // min/max GID are always computed (from indexBase and numGlobal), and are therefore globally coherent as long as deps are.
00054 // indexBase and numGlobal must always be verified.
00055 // isContiguous is true for the "easy" constructors, assume false for the expert constructor
00056 // 
00057 // so we explicitly check    : isCont, numGlobal, indexBase
00058 // then we explicitly compute: MinAllGID, MaxAllGID
00059 // Data explicitly checks    : isDistributed
00060 
00061 namespace Tpetra {
00062 
00063   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00064   Map<LocalOrdinal,GlobalOrdinal,Node>::Map(
00065                         global_size_t numGlobalElements_in, 
00066                         GlobalOrdinal indexBase_in, 
00067                         const Teuchos::RCP<const Teuchos::Comm<int> > &comm_in,
00068                         LocalGlobal lOrG, 
00069                         const Teuchos::RCP<Node> &node_in)
00070   : comm_(comm_in), 
00071     node_(node_in) {
00072     // distribute the elements across the nodes so that they are 
00073     // - non-overlapping
00074     // - contiguous
00075     // - as evenly distributed as possible
00076     using Teuchos::as;
00077     using Teuchos::outArg;
00078     const global_size_t GST0 = Teuchos::OrdinalTraits<global_size_t>::zero();
00079     const global_size_t GST1 = Teuchos::OrdinalTraits<global_size_t>::one();
00080     const global_size_t GSTI = Teuchos::OrdinalTraits<global_size_t>::invalid();
00081     const GlobalOrdinal G1 = Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00082 
00083     std::string errPrefix;
00084     errPrefix = Teuchos::typeName(*this) + "::constructor(numGlobal,indexBase,comm,lOrG): ";
00085 
00086     if (lOrG == GloballyDistributed) {
00087       const int numImages = comm_->getSize();
00088       const int myImageID = comm_->getRank();
00089 
00090       // check that numGlobalElements,indexBase is equivalent across images
00091       global_size_t rootNGE = numGlobalElements_in;
00092       GlobalOrdinal rootIB  = indexBase_in;
00093       Teuchos::broadcast<int,global_size_t>(*comm_,0,&rootNGE);
00094       Teuchos::broadcast<int,GlobalOrdinal>(*comm_,0,&rootIB);
00095       int localChecks[2], globalChecks[2];
00096       localChecks[0] = -1;   // fail or pass
00097       localChecks[1] = 0;    // fail reason
00098       if (numGlobalElements_in != rootNGE) {
00099         localChecks[0] = myImageID;
00100         localChecks[1] = 1;
00101       }
00102       else if (indexBase_in != rootIB) {
00103         localChecks[0] = myImageID;
00104         localChecks[1] = 2;
00105       }
00106       // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass, as well as the reason
00107       // these will be -1 and 0 if all procs passed
00108       Teuchos::reduceAll<int,int>(*comm_,Teuchos::REDUCE_MAX,2,localChecks,globalChecks);
00109       if (globalChecks[0] != -1) {
00110         if (globalChecks[1] == 1) {
00111           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00112               errPrefix << "numGlobal must be the same on all nodes (examine node " << globalChecks[0] << ").");
00113         }
00114         else if (globalChecks[1] == 2) {
00115           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00116               errPrefix << "indexBase must be the same on all nodes (examine node " << globalChecks[0] << ").");
00117         }
00118         else {
00119           // logic error on our part
00120           TEST_FOR_EXCEPTION(true,std::logic_error,
00121               errPrefix << "logic error. Please contact the Tpetra team.");
00122         }
00123       }
00124       // numGlobalElements is coherent, but is it valid? this comparison looks funny, but it avoids compiler warnings on unsigned types.
00125       TEST_FOR_EXCEPTION((numGlobalElements_in < GST1 && numGlobalElements_in != GST0) || numGlobalElements_in == GSTI, std::invalid_argument,
00126           errPrefix << "numGlobalElements (== " << rootNGE << ") must be >= 0.");
00127 
00128       indexBase_ = rootIB;
00129       numGlobalElements_ = rootNGE;
00130 
00131       /* compute numLocalElements
00132          We can write numGlobalElements as follows:
00133          numGlobalElements == numImages * B + remainder
00134          Each image is allocated elements as follows:
00135          [ B+1    iff myImageID <  remainder
00136          numLocalElements = [
00137          [ B      iff myImageID >= remainder
00138          In the case that remainder == 0, then all images fall into the 
00139          latter case: numLocalElements == B == numGlobalElements / numImages
00140          It can then be shown that 
00141          numImages
00142          \Sum      numLocalElements_i  == numGlobalElements
00143          i=0
00144          This strategy is simple, requires no communication, and is optimal vis-a-vis
00145          uniform distribution of elements.
00146          This strategy is valid for any value of numGlobalElements and numImages, 
00147          including the following border cases:
00148          - numImages == 1         -> remainder == 0 && numGlobalElements == numLocalElements
00149          - numelements < numImages -> remainder == numGlobalElements && numLocalElements \in [0,1]
00150        */
00151       numLocalElements_ = as<size_t>(numGlobalElements_ / as<global_size_t>(numImages));
00152       int remainder = as<int>(numGlobalElements_ % as<global_size_t>(numImages));
00153 #ifdef HAVE_TEUCHOS_DEBUG
00154       // the above code assumes truncation. is that safe?
00155       SHARED_TEST_FOR_EXCEPTION(numLocalElements_ * numImages + remainder != numGlobalElements_,
00156           std::logic_error, "Tpetra::Map::constructor(numGlobal,indexBase,platform): GlobalOrdinal does not implement division with truncation."
00157           << " Please contact Tpetra team.",*comm_);
00158 #endif
00159       GlobalOrdinal start_index;
00160       if (myImageID < remainder) {
00161         ++numLocalElements_;
00162         /* the myImageID images before were each allocated 
00163            numGlobalElements/numImages+1
00164            ergo, my offset is:
00165            myImageID * (numGlobalElements/numImages+1)
00166            == myImageID * numLocalElements
00167          */
00168         start_index = as<GlobalOrdinal>(myImageID) * as<GlobalOrdinal>(numLocalElements_);
00169       }
00170       else {
00171         /* a quantity (equal to remainder) of the images before me
00172            were each allocated 
00173            numGlobalElements/numImages+1
00174            elements. a quantity (equal to myImageID-remainder) of the remaining 
00175            images before me were each allocated 
00176            numGlobalElements/numImages
00177            elements. ergo, my offset is:
00178            remainder*(numGlobalElements/numImages+1) + (myImageID-remainder)*numGlobalElements/numImages
00179            == remainder*numLocalElements + remainder + myImageID*numLocalElements - remainder*numLocalElements
00180            == myImageID*numLocalElements + remainder
00181          */
00182         start_index = as<GlobalOrdinal>(myImageID)*as<GlobalOrdinal>(numLocalElements_) + as<GlobalOrdinal>(remainder);
00183       }
00184 
00185       // compute the min/max global IDs
00186       minMyGID_  = start_index + indexBase_;
00187       maxMyGID_  = minMyGID_ + numLocalElements_ - G1;
00188       minAllGID_ = indexBase_;
00189       maxAllGID_ = indexBase_ + numGlobalElements_ - G1;
00190       contiguous_ = true;
00191       distributed_ = (numImages > 1 ? true : false);
00192       setupDirectory();
00193     }
00194     else {  // lOrG == LocallyReplicated
00195       // compute the min/max global IDs
00196       indexBase_ = indexBase_in;
00197       numGlobalElements_ = numGlobalElements_in;
00198       numLocalElements_  = as<size_t>(numGlobalElements_in);
00199       minAllGID_ = indexBase_;
00200       maxAllGID_ = indexBase_ + numGlobalElements_ - G1;
00201       minMyGID_  = minAllGID_;
00202       maxMyGID_  = maxAllGID_;
00203       contiguous_ = true;
00204       distributed_ = false;
00205     }
00206   }
00207 
00208   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00209   Map<LocalOrdinal,GlobalOrdinal,Node>::Map(global_size_t numGlobalElements_in, size_t numLocalElements_in, GlobalOrdinal indexBase_in, 
00210                                             const Teuchos::RCP<const Teuchos::Comm<int> > &comm_in,
00211                                             const Teuchos::RCP<Node> &node_in) 
00212   : comm_(comm_in), 
00213     node_(node_in) {
00214     // Distribute the elements across the nodes so that they are 
00215     // - non-overlapping
00216     // - contiguous
00217     // This differs from Map(Ord,Ord,Plat) in that the user has specified the number of elements 
00218     // per node, so that they are not (necessarily) evenly distributed
00219 
00220     using Teuchos::outArg;
00221 
00222     const size_t  L0 = Teuchos::OrdinalTraits<size_t>::zero();
00223     const size_t  L1 = Teuchos::OrdinalTraits<size_t>::one();
00224     const global_size_t GST0 = Teuchos::OrdinalTraits<global_size_t>::zero();
00225     const global_size_t GST1 = Teuchos::OrdinalTraits<global_size_t>::one();
00226     const global_size_t GSTI = Teuchos::OrdinalTraits<global_size_t>::invalid();
00227     const GlobalOrdinal G1 = Teuchos::OrdinalTraits<GlobalOrdinal>::one();
00228 
00229     std::string errPrefix;
00230     errPrefix = Teuchos::typeName(*this) + "::constructor(numGlobal,numLocal,indexBase,platform): ";
00231 
00232     // get a internodal communicator from the Platform
00233     const int myImageID = comm_->getRank();
00234 
00235     { // begin scoping block
00236       // for communicating failures 
00237       int localChecks[2], globalChecks[2];
00238       /* compute the global size 
00239          we are computing the number of global elements because exactly ONE of the following is true:
00240          - the user didn't specify it, and we need it
00241          - the user did specify it, but we need to 
00242            + validate it against the sum of the local sizes, and
00243            + ensure that it is the same on all nodes
00244        */
00245       global_size_t global_sum;
00246       Teuchos::reduceAll<int,global_size_t>(*comm_,Teuchos::REDUCE_SUM,
00247         Teuchos::as<global_size_t>(numLocalElements_in),outArg(global_sum));
00248       /* there are three errors we should be detecting:
00249          - numGlobalElements != invalid() and it is incorrect/invalid
00250          - numLocalElements invalid (<0)
00251       */
00252       localChecks[0] = -1;
00253       localChecks[1] = 0;
00254       if (numLocalElements_in < L1 && numLocalElements_in != L0) {
00255         // invalid
00256         localChecks[0] = myImageID;
00257         localChecks[1] = 1;
00258       }
00259       else if (numGlobalElements_in < GST1 && numGlobalElements_in != GST0 && numGlobalElements_in != GSTI) {
00260         // invalid
00261         localChecks[0] = myImageID;
00262         localChecks[1] = 2;
00263       }
00264       else if (numGlobalElements_in != GSTI && numGlobalElements_in != global_sum) {
00265         // incorrect
00266         localChecks[0] = myImageID;
00267         localChecks[1] = 3;
00268       }
00269       // now check that indexBase is equivalent across images
00270       GlobalOrdinal rootIB = indexBase_in;
00271       Teuchos::broadcast<int,GlobalOrdinal>(*comm_,0,&rootIB);   // broadcast one ordinal from node 0
00272       if (indexBase_in != rootIB) {
00273         localChecks[0] = myImageID;
00274         localChecks[1] = 4;
00275       }
00276       // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass
00277       // this will be -1 if all procs passed
00278       Teuchos::reduceAll<int,int>(*comm_,Teuchos::REDUCE_MAX,2,localChecks,globalChecks);
00279       if (globalChecks[0] != -1) {
00280         if (globalChecks[1] == 1) {
00281           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00282               errPrefix << "numLocal is not valid on at least one node (possibly node " 
00283               << globalChecks[0] << ").");
00284         }
00285         else if (globalChecks[1] == 2) {
00286           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00287               errPrefix << "numGlobal is not valid on at least one node (possibly node " 
00288               << globalChecks[0] << ").");
00289         }
00290         else if (globalChecks[1] == 3) {
00291           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00292               errPrefix << "numGlobal doesn't match sum of numLocal (== " 
00293               << global_sum << ") on at least one node (possibly node " 
00294               << globalChecks[0] << ").");
00295         }
00296         else if (globalChecks[1] == 4) {
00297           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00298               errPrefix << "indexBase is not the same on all nodes (examine node " 
00299               << globalChecks[0] << ").");
00300         }
00301         else {
00302           // logic error on my part
00303           TEST_FOR_EXCEPTION(true,std::logic_error,
00304               errPrefix << "logic error. Please contact the Tpetra team.");
00305         }
00306       }
00307       // set numGlobalElements
00308       if (numGlobalElements_in == GSTI) {
00309         numGlobalElements_ = global_sum;
00310       }
00311       else {
00312         numGlobalElements_ = numGlobalElements_in;
00313       }
00314       numLocalElements_ = numLocalElements_in;
00315       indexBase_ = indexBase_in;
00316     } // end of scoping block
00317 
00318     // compute my local offset
00319     GlobalOrdinal start_index;
00320     Teuchos::scan<int,GlobalOrdinal>(*comm_,Teuchos::REDUCE_SUM,numLocalElements_,Teuchos::outArg(start_index));
00321     start_index -= numLocalElements_;
00322 
00323     minAllGID_ = indexBase_;
00324     maxAllGID_ = indexBase_ + numGlobalElements_ - G1;
00325     minMyGID_ = start_index + indexBase_;
00326     maxMyGID_ = minMyGID_ + numLocalElements_ - G1;
00327     contiguous_ = true;
00328     distributed_ = checkIsDist();
00329     setupDirectory();
00330   }
00331 
00332   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00333   Map<LocalOrdinal,GlobalOrdinal,Node>::Map (global_size_t numGlobalElements_in, const Teuchos::ArrayView<const GlobalOrdinal> &entryList, GlobalOrdinal indexBase_in, 
00334                                              const Teuchos::RCP<const Teuchos::Comm<int> > &comm_in,
00335                                              const Teuchos::RCP<Node> &node_in)
00336   : comm_(comm_in),
00337     node_(node_in) {
00338     using Teuchos::as;
00339     using Teuchos::outArg;
00340     // Distribute the elements across the nodes in an arbitrary user-specified manner
00341     // They are not necessarily contiguous or evenly distributed
00342     const size_t  L0 = Teuchos::OrdinalTraits<size_t>::zero();
00343     const global_size_t GST0 = Teuchos::OrdinalTraits<global_size_t>::zero();
00344     const global_size_t GST1 = Teuchos::OrdinalTraits<global_size_t>::one();
00345     const global_size_t GSTI = Teuchos::OrdinalTraits<global_size_t>::invalid();
00346 
00347     LocalOrdinal numLocalElements_in = Teuchos::as<LocalOrdinal>(entryList.size());
00348 
00349     std::string errPrefix;
00350     errPrefix = Teuchos::typeName(*this) + "::constructor(numGlobal,entryList,indexBase,platform): ";
00351 
00352     const int myImageID = comm_->getRank();
00353     { // begin scoping block
00354       // for communicating failures 
00355       int localChecks[2], globalChecks[2];
00356 
00357       /* compute the global size 
00358          we are computing the number of global elements because exactly ONE of the following is true:
00359          - the user didn't specify it, and we need it
00360          - the user did specify it, but we need to 
00361            + validate it against the sum of the local sizes, and
00362            + ensure that it is the same on all nodes
00363        */
00364       global_size_t global_sum;
00365       Teuchos::reduceAll<int,global_size_t>(*comm_,Teuchos::REDUCE_SUM,
00366         as<global_size_t>(numLocalElements_in),outArg(global_sum));
00367       localChecks[0] = -1;
00368       localChecks[1] = 0;
00369       if (numGlobalElements_in < GST1 && numGlobalElements_in != GST0 && numGlobalElements_in != GSTI) {
00370         // invalid
00371         localChecks[0] = myImageID;
00372         localChecks[1] = 1;
00373       }
00374       else if (numGlobalElements_in != GSTI && numGlobalElements_in != global_sum) {
00375         // incorrect
00376         localChecks[0] = myImageID;
00377         localChecks[1] = 2;
00378       }
00379       // now check that indexBase is equivalent across images
00380       GlobalOrdinal rootIB = indexBase_in;
00381       Teuchos::broadcast<int,GlobalOrdinal>(*comm_,0,&rootIB);   // broadcast one ordinal from node 0
00382       if (indexBase_in != rootIB) {
00383         localChecks[0] = myImageID;
00384         localChecks[1] = 3;
00385       }
00386       // REDUCE_MAX will give us the image ID of the highest rank proc that DID NOT pass
00387       // this will be -1 if all procs passed
00388       Teuchos::reduceAll<int,int>(*comm_,Teuchos::REDUCE_MAX,2,localChecks,globalChecks);
00389       if (globalChecks[0] != -1) {
00390         if (globalChecks[1] == 1) {
00391           TEST_FOR_EXCEPTION(true,std::invalid_argument,
00392               errPrefix << "numGlobal is not valid on at least one node (possibly node "
00393               << globalChecks[0] << ").");
00394         }
00395         else if (globalChecks[1] == 2) {
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] == 3) {
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::logic_error,
00409               errPrefix << "logic error. Please contact the Tpetra team.");
00410         }
00411       }
00412 
00413       // these are all validated/computed now
00414       if (numGlobalElements_in == GSTI) {
00415         numGlobalElements_ = global_sum;
00416       }
00417       else {
00418         numGlobalElements_ = numGlobalElements_in;
00419       }
00420       numLocalElements_ = numLocalElements_in;
00421       indexBase_ = indexBase_in;
00422     } // end scoping block
00423 
00424     // assume for now that there are numLocalElements (there may be less, if some
00425     // GIDs are duplicated in entryList)
00426     minMyGID_ = indexBase_;
00427     maxMyGID_ = indexBase_;
00428     // create the GID to LID map; do not assume GID in entryList are distinct.
00429     // in the case that a GID is duplicated, keep the previous LID
00430     // this is necessary so that LIDs are in [0,numLocal)
00431     size_t numUniqueGIDs = 0;
00432     if (numLocalElements_ > L0) {
00433       lgMap_ = Teuchos::arcp<GlobalOrdinal>(numLocalElements_);
00434       for (size_t i=0; i < numLocalElements_; i++) {
00435         if (glMap_.find(entryList[i]) == glMap_.end()) {
00436           lgMap_[numUniqueGIDs] = entryList[i];   // lgMap_:  LID to GID
00437           glMap_[entryList[i]] = numUniqueGIDs;   // glMap_: GID to LID
00438           numUniqueGIDs++;
00439         }
00440       }
00441       // shrink lgMap appropriately
00442       if (numLocalElements_ != numUniqueGIDs) {
00443         numLocalElements_ = numUniqueGIDs;
00444         lgMap_ = lgMap_.persistingView(0,numLocalElements_);
00445       }
00446       minMyGID_ = *std::min_element(lgMap_.begin(), lgMap_.end());
00447       maxMyGID_ = *std::max_element(lgMap_.begin(), lgMap_.end());
00448     }
00449 
00450     // set min/maxAllGIDs
00451     Teuchos::reduceAll<int,GlobalOrdinal>(*comm_,Teuchos::REDUCE_MIN,minMyGID_,Teuchos::outArg(minAllGID_));
00452     Teuchos::reduceAll<int,GlobalOrdinal>(*comm_,Teuchos::REDUCE_MAX,maxMyGID_,Teuchos::outArg(maxAllGID_));
00453     contiguous_  = false;
00454     distributed_ = checkIsDist();
00455     TEST_FOR_EXCEPTION(minAllGID_ < indexBase_, std::invalid_argument,
00456         errPrefix << "minimum GID (== " << minAllGID_ << ") is less than indexBase (== " << indexBase_ << ")");
00457     setupDirectory();
00458   }
00459 
00460   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00461   Map<LocalOrdinal,GlobalOrdinal,Node>::~Map () 
00462   {}
00463 
00464   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00465   LocalOrdinal Map<LocalOrdinal,GlobalOrdinal,Node>::getLocalElement(GlobalOrdinal globalIndex) const {
00466     if (contiguous_) {
00467       if (globalIndex < getMinGlobalIndex() || globalIndex > getMaxGlobalIndex()) {
00468         return Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00469       }
00470       return Teuchos::as<LocalOrdinal>(globalIndex - getMinGlobalIndex());
00471     }
00472     else {
00473       typename std::map<GlobalOrdinal,LocalOrdinal>::const_iterator i;
00474       i = glMap_.find(globalIndex);
00475       if (i == glMap_.end()) {
00476         return Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
00477       }
00478       return i->second;
00479     }
00480   }
00481 
00482   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00483   GlobalOrdinal Map<LocalOrdinal,GlobalOrdinal,Node>::getGlobalElement(LocalOrdinal localIndex) const {
00484     if (localIndex < getMinLocalIndex() || localIndex > getMaxLocalIndex()) {
00485       return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
00486     }
00487     if (contiguous_) {
00488       return getMinGlobalIndex() + localIndex;
00489     }
00490     else {
00491       return lgMap_[localIndex];
00492     }
00493   }
00494 
00495   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00496   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isNodeLocalElement(LocalOrdinal localIndex) const {
00497     if (localIndex < getMinLocalIndex() || localIndex > getMaxLocalIndex()) {
00498       return false;
00499     }
00500     return true;
00501   }
00502 
00503   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00504   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isNodeGlobalElement(GlobalOrdinal globalIndex) const {
00505     if (contiguous_) {
00506       return (getMinGlobalIndex() <= globalIndex) && (globalIndex <= getMaxGlobalIndex());
00507     }
00508     else {
00509       typename std::map<GlobalOrdinal,LocalOrdinal>::const_iterator i;
00510       i = glMap_.find(globalIndex);
00511       return (i != glMap_.end());
00512     }
00513   }
00514 
00515   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00516   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isContiguous() const {
00517     return contiguous_;
00518   }
00519 
00520   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00521   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isCompatible (const Map<LocalOrdinal,GlobalOrdinal,Node> &map) const {
00522     using Teuchos::outArg;
00523     // check to make sure distribution is the same
00524     char iscompat_lcl;
00525     if (getGlobalNumElements() != map.getGlobalNumElements() ||
00526           getNodeNumElements() != map.getNodeNumElements()) {
00527       // NOT compat on this node
00528       iscompat_lcl = 0;
00529     }
00530     else {
00531       // compat on this node
00532       iscompat_lcl = 1;
00533     }
00534     char iscompat_gbl;
00535     Teuchos::reduceAll<int,char>(*comm_,Teuchos::REDUCE_MIN,iscompat_lcl,
00536       outArg(iscompat_gbl));
00537     return (iscompat_gbl == 1);
00538   }
00539 
00540   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00541   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isSameAs(const Map<LocalOrdinal,GlobalOrdinal,Node> &map) const {
00542     using Teuchos::outArg;
00543     if (this == &map) {
00544       // we should assume that this is globally coherent
00545       // if they share the same underlying MapData, then they must be equivalent maps
00546       return true;
00547     }
00548 
00549     // check all other globally coherent properties
00550     // if they do not share each of these properties, then they cannot be 
00551     // equivalent maps
00552     if ( (getMinAllGlobalIndex()   != map.getMinAllGlobalIndex())   ||
00553          (getMaxAllGlobalIndex()   != map.getMaxAllGlobalIndex())   ||
00554          (getGlobalNumElements() != map.getGlobalNumElements()) ||
00555          (isDistributed()       != map.isDistributed())       || 
00556          (getIndexBase()        != map.getIndexBase())         )  {
00557       return false;
00558     }
00559 
00560     // If we get this far, we need to check local properties and the 
00561     // communicate same-ness across all nodes
00562     // we prefer local work over communication, ergo, we will perform all
00563     // comparisons and conduct a single communication
00564     char isSame_lcl = 1;
00565 
00566     // check number of entries owned by this node
00567     if (getNodeNumElements() != map.getNodeNumElements()) {
00568       isSame_lcl = 0;
00569     }
00570 
00571     // check the identity of the entries owned by this node
00572     // only do this if we haven't already determined not-same-ness
00573     if (isSame_lcl == 1) {
00574       // if they are contiguous, we can check the ranges easily
00575       // if they are not contiguous, we must check the individual LID -> GID mappings
00576       // the latter approach is valid in either case, but the former is faster
00577       if (contiguous_ == true && map.contiguous_ == true) {
00578         if ( (getMinGlobalIndex() != map.getMinGlobalIndex()) ||
00579              (getMaxGlobalIndex() != map.getMaxGlobalIndex()) ){
00580           isSame_lcl = 0;
00581         }
00582       }
00583       else {
00584         /* Note: std::equal requires that the latter range is as large as the former.
00585          * We know the ranges have equal length, because they have the same number of 
00586          * local entries. 
00587          * However, we only know that lgMap_ has been filled for the Map that is not
00588          * contiguous (one is potentially contiguous.) Calling getNodeElementList()
00589          * will create it. */
00590         Teuchos::ArrayView<const GlobalOrdinal> ge1, ge2;
00591         ge1 =     getNodeElementList();
00592         ge2 = map.getNodeElementList();
00593         if (!std::equal(ge1.begin(),ge1.end(),ge2.begin())) {
00594           isSame_lcl = 0;
00595         }
00596       }
00597     }
00598 
00599     // now, determine if we detected not-same-ness on any node
00600     char isSame_gbl;
00601     Teuchos::reduceAll<int,char>(*comm_,Teuchos::REDUCE_MIN,isSame_lcl,outArg(isSame_gbl));
00602     return (isSame_gbl == 1);
00603   }
00604 
00605   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00606   Teuchos::ArrayView<const GlobalOrdinal>
00607   Map<LocalOrdinal,GlobalOrdinal,Node>::getNodeElementList() const {
00608     // if so (and we have local entries), then fill it.
00609     if (lgMap_ == Teuchos::null && numLocalElements_ > 0) {
00610 #ifdef HAVE_TEUCHOS_DEBUG
00611       // this would have been set up for a non-contiguous map
00612       TEST_FOR_EXCEPTION(contiguous_ != true, std::logic_error,
00613           "Tpetra::Map::getNodeElementList: logic error. Please notify the Tpetra team.");
00614 #endif
00615       lgMap_ = Teuchos::arcp<GlobalOrdinal>(numLocalElements_);
00616       Teuchos::ArrayRCP<GlobalOrdinal> lgptr = lgMap_;
00617       for (GlobalOrdinal gid=minMyGID_; gid <= maxMyGID_; ++gid) {
00618         *(lgptr++) = gid;
00619       }
00620     }
00621     return lgMap_();
00622   }
00623 
00624   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00625   bool Map<LocalOrdinal,GlobalOrdinal,Node>::isDistributed() const {
00626     return distributed_;
00627   }
00628 
00629   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00630   std::string Map<LocalOrdinal,GlobalOrdinal,Node>::description() const {
00631     std::ostringstream oss;
00632     oss << Teuchos::Describable::description();
00633     oss << "{getGlobalNumElements() = " << getGlobalNumElements()
00634         << ", getNodeNumElements() = " << getNodeNumElements()
00635         << ", isContiguous() = " << isContiguous()
00636         << ", isDistributed() = " << isDistributed()
00637         << "}";
00638     return oss.str();
00639   }
00640 
00641   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00642   void Map<LocalOrdinal,GlobalOrdinal,Node>::describe( Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const {
00643     using std::endl;
00644     using std::setw;
00645     using Teuchos::VERB_DEFAULT;
00646     using Teuchos::VERB_NONE;
00647     using Teuchos::VERB_LOW;
00648     using Teuchos::VERB_MEDIUM;
00649     using Teuchos::VERB_HIGH;
00650     using Teuchos::VERB_EXTREME;
00651 
00652     const size_t nME = getNodeNumElements();
00653     Teuchos::ArrayView<const GlobalOrdinal> myEntries = getNodeElementList();
00654     int myImageID = comm_->getRank();
00655     int numImages = comm_->getSize();
00656 
00657     Teuchos::EVerbosityLevel vl = verbLevel;
00658     if (vl == VERB_DEFAULT) vl = VERB_LOW;
00659 
00660     size_t width = 1;
00661     for (size_t dec=10; dec<getGlobalNumElements(); dec *= 10) {
00662       ++width;
00663     }
00664     width = std::max<size_t>(width,12) + 2;
00665 
00666     Teuchos::OSTab tab(out);
00667 
00668     if (vl == VERB_NONE) {
00669       // do nothing
00670     }
00671     else if (vl == VERB_LOW) {
00672       out << this->description() << endl;
00673     }
00674     else {  // MEDIUM, HIGH or EXTREME
00675       for (int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
00676         if (myImageID == imageCtr) {
00677           if (myImageID == 0) { // this is the root node (only output this info once)
00678             out << endl 
00679                 << "Number of Global Entries = " << getGlobalNumElements()  << endl
00680                 << "Maximum of all GIDs      = " << getMaxAllGlobalIndex() << endl
00681                 << "Minimum of all GIDs      = " << getMinAllGlobalIndex() << endl
00682                 << "Index Base               = " << getIndexBase()         << endl;
00683           }
00684           out << endl;
00685           if (vl == VERB_HIGH || vl == VERB_EXTREME) {
00686             out << "Number of Local Elements   = " << nME           << endl
00687                 << "Maximum of my GIDs         = " << getMaxGlobalIndex() << endl
00688                 << "Minimum of my GIDs         = " << getMinGlobalIndex() << endl;
00689             out << endl;
00690           }
00691           if (vl == VERB_EXTREME) {
00692             out << std::setw(width) << "Node ID"
00693                 << std::setw(width) << "Local Index"
00694                 << std::setw(width) << "Global Index"
00695                 << endl;
00696             for (size_t i=0; i < nME; i++) {
00697               out << std::setw(width) << myImageID 
00698                   << std::setw(width) << i
00699                   << std::setw(width) << myEntries[i]
00700                   << endl;
00701             }
00702             out << std::flush;
00703           }
00704         }
00705         // Do a few global ops to give I/O a chance to complete
00706         comm_->barrier();
00707         comm_->barrier();
00708         comm_->barrier();
00709       }
00710     }
00711   }
00712 
00713   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00714   void Map<LocalOrdinal,GlobalOrdinal,Node>::setupDirectory() {
00715     if (getGlobalNumElements() != Teuchos::OrdinalTraits<global_size_t>::zero()) {
00716       if (directory_ == Teuchos::null) {
00717         directory_ = Teuchos::rcp( new Directory<LocalOrdinal,GlobalOrdinal,Node>(Teuchos::rcp(this,false)) );
00718       }
00719     }
00720   }
00721 
00722   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00723   LookupStatus Map<LocalOrdinal,GlobalOrdinal,Node>::getRemoteIndexList(
00724                     const Teuchos::ArrayView<const GlobalOrdinal> & GIDList, 
00725                     const Teuchos::ArrayView<int> & imageIDList, 
00726                     const Teuchos::ArrayView<LocalOrdinal> & LIDList) const {
00727     if (distributed_ == false) {
00728       TEST_FOR_EXCEPTION(GIDList.size() > 0, std::runtime_error,
00729         Teuchos::typeName(*this) << "::getRemoteIndexList() cannot be called for local maps.");
00730       return AllIDsPresent;
00731     }
00732     TEST_FOR_EXCEPTION(GIDList.size() > 0 && getGlobalNumElements() == 0, std::runtime_error,
00733         Teuchos::typeName(*this) << "::getRemoteIndexList(): getRemoteIndexList() cannot be called, zero entries in Map.");
00734     return directory_->getDirectoryEntries(GIDList, imageIDList, LIDList);
00735   }
00736 
00737   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00738   LookupStatus Map<LocalOrdinal,GlobalOrdinal,Node>::getRemoteIndexList(
00739                     const Teuchos::ArrayView<const GlobalOrdinal> & GIDList, 
00740                     const Teuchos::ArrayView<int> & imageIDList) const {
00741     if (distributed_ == false) {
00742       TEST_FOR_EXCEPTION(GIDList.size() > 0, std::runtime_error,
00743         Teuchos::typeName(*this) << "::getRemoteIndexList() cannot be called for local maps.");
00744       return AllIDsPresent;
00745     }
00746     TEST_FOR_EXCEPTION(GIDList.size() > 0 && getGlobalNumElements() == 0, std::runtime_error,
00747         Teuchos::typeName(*this) << "::getRemoteIndexList(): getRemoteIndexList() cannot be called, zero entries in Map.");
00748     return directory_->getDirectoryEntries(GIDList, imageIDList);
00749   }
00750 
00751   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00752   const Teuchos::RCP<const Teuchos::Comm<int> > &
00753   Map<LocalOrdinal,GlobalOrdinal,Node>::getComm() const {
00754     return comm_;
00755   }
00756 
00757   template <class LocalOrdinal, class GlobalOrdinal, class Node>
00758   const Teuchos::RCP<Node> &
00759   Map<LocalOrdinal,GlobalOrdinal,Node>::getNode() const {
00760     return node_;
00761   }
00762 
00763   template <class LocalOrdinal,class GlobalOrdinal, class Node>
00764   bool Map<LocalOrdinal,GlobalOrdinal,Node>::checkIsDist() const {
00765     using Teuchos::outArg;
00766     bool global = false;
00767     if(comm_->getSize() > 1) {
00768       char localRep = 0;
00769       if (numGlobalElements_ == Teuchos::as<global_size_t>(numLocalElements_)) {
00770         localRep = 1;
00771       }
00772       char allLocalRep;
00773       Teuchos::reduceAll<int>(*comm_,Teuchos::REDUCE_MIN,localRep,outArg(allLocalRep));
00774       if (allLocalRep != 1) {
00775         global = true;
00776       }
00777     }
00778     return global;
00779   }
00780 
00781 } // Tpetra namespace
00782 
00783 template <class LocalOrdinal, class GlobalOrdinal>
00784 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType> >
00785 Tpetra::createLocalMap(size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm) {
00786   return createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType>(numElements, comm, Kokkos::DefaultNode::getDefaultNode());
00787 }
00788 
00789 template <class LocalOrdinal, class GlobalOrdinal>
00790 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType> >
00791 Tpetra::createUniformContigMap(global_size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm) {
00792   return createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType>(numElements, comm, Kokkos::DefaultNode::getDefaultNode());
00793 }
00794 
00795 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00796 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
00797 Tpetra::createUniformContigMapWithNode(global_size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP<Node> &node) {
00798   Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map;
00799   map = Teuchos::rcp( new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(numElements,                                    // num elements, global and local
00800                                                               Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),  // index base is zero
00801                                                               comm, GloballyDistributed, node) );
00802   return map.getConst();
00803 }
00804 
00805 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00806 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
00807 Tpetra::createLocalMapWithNode(size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< Node > &node) {
00808   Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map;
00809   map = Teuchos::rcp( new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>((Tpetra::global_size_t)numElements,                     // num elements, global and local
00810                                                               Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),  // index base is zero
00811                                                               comm, LocallyReplicated, node) );
00812   return map.getConst();
00813 }
00814 
00815 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00816 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
00817 Tpetra::createContigMapWithNode(Tpetra::global_size_t numElements, size_t localNumElements, 
00818                                 const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< Node > &node) {
00819   Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map;
00820   map = Teuchos::rcp( new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(numElements,localNumElements,
00821                                                               Teuchos::OrdinalTraits<GlobalOrdinal>::zero(),  // index base is zero
00822                                                               comm, node) );
00823   return map.getConst();
00824 }
00825 
00826 template <class LocalOrdinal, class GlobalOrdinal>
00827 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType> >
00828 Tpetra::createContigMap(Tpetra::global_size_t numElements, size_t localNumElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm) {
00829   return Tpetra::createContigMapWithNode<LocalOrdinal,GlobalOrdinal,Kokkos::DefaultNode::DefaultNodeType>(numElements, localNumElements, comm, Kokkos::DefaultNode::getDefaultNode() );
00830 }
00831 
00832 
00833 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00834 Teuchos::RCP< const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
00835 Tpetra::createWeightedContigMapWithNode(int myWeight, Tpetra::global_size_t numElements, 
00836                                         const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< Node > &node) {
00837   Teuchos::RCP< Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > map;
00838   int sumOfWeights, elemsLeft, localNumElements;
00839   const int numImages = comm->getSize(), 
00840             myImageID = comm->getRank();
00841   Teuchos::reduceAll<int>(*comm,Teuchos::REDUCE_SUM,myWeight,Teuchos::outArg(sumOfWeights));
00842   const double myShare = ((double)myWeight) / ((double)sumOfWeights);
00843   localNumElements = (int)std::floor( myShare * ((double)numElements) );
00844   // std::cout << "numElements: " << numElements << "  myWeight: " << myWeight << "  sumOfWeights: " << sumOfWeights << "  myShare: " << myShare << std::endl;
00845   Teuchos::reduceAll<int>(*comm,Teuchos::REDUCE_SUM,localNumElements,Teuchos::outArg(elemsLeft));
00846   elemsLeft = numElements - elemsLeft;
00847   // std::cout << "(before) localNumElements: " << localNumElements << "  elemsLeft: " << elemsLeft << std::endl;
00848   // i think this is true. just test it for now.
00849   TEST_FOR_EXCEPT(elemsLeft < -numImages || numImages < elemsLeft);
00850   if (elemsLeft < 0) {
00851     // last elemsLeft nodes lose an element
00852     if (myImageID >= numImages-elemsLeft) --localNumElements;
00853   }
00854   else if (elemsLeft > 0) {
00855     // first elemsLeft nodes gain an element
00856     if (myImageID < elemsLeft) ++localNumElements;
00857   }
00858   // std::cout << "(after) localNumElements: " << localNumElements << std::endl;
00859   return createContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(numElements,localNumElements,comm,node);
00860 }
00861 
00863 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00864 bool operator== (const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &map1, const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &map2) {
00865   return map1.isSameAs(map2);
00866 }
00867 
00869 template <class LocalOrdinal, class GlobalOrdinal, class Node>
00870 bool operator!= (const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &map1, const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> &map2) {
00871   return !map1.isSameAs(map2);
00872 }
00873 
00874 
00875 //
00876 // Explicit instantiation macro
00877 //
00878 // Must be expanded from within the Tpetra namespace!
00879 //
00880 
00882 #define TPETRA_MAP_INSTANT(LO,GO,NODE) \
00883   \
00884   template class Map< LO , GO , NODE >; \
00885   \
00886   template Teuchos::RCP< const Map<LO,GO,NODE> > \
00887   createLocalMapWithNode<LO,GO,NODE>(size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< NODE > &node); \
00888   \
00889   template Teuchos::RCP< const Map<LO,GO,NODE> > \
00890   createContigMapWithNode<LO,GO,NODE>(global_size_t numElements, size_t localNumElements, \
00891                                       const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< NODE > &node); \
00892   \
00893   template Teuchos::RCP< const Map<LO,GO,NODE> > \
00894   createUniformContigMapWithNode<LO,GO,NODE>(global_size_t numElements, \
00895                                              const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< NODE > &node); \
00896   \
00897   template Teuchos::RCP< const Map<LO,GO,NODE> > \
00898   createWeightedContigMapWithNode<LO,GO,NODE>(int thisNodeWeight, global_size_t numElements, \
00899                                               const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< NODE > &node); \
00900 
00901 
00902 #endif // TPETRA_MAP_DEF_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