Tpetra_Map_def.hpp

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

Generated on Tue Jul 13 09:39:06 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7