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