Tpetra Matrix/Vector Services Version of the Day
Tpetra_DirectoryImpl_def.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef __Tpetra_DirectoryImpl_def_hpp
00043 #define __Tpetra_DirectoryImpl_def_hpp
00044 
00045 #ifdef DOXYGEN_USE_ONLY
00046 #  include <Tpetra_DirectoryImpl_decl.hpp>
00047 #endif
00048 #include <Tpetra_Distributor.hpp>
00049 #include <Tpetra_Map.hpp>
00050 #include <Tpetra_TieBreak.hpp>
00051 
00052 #include <Tpetra_Details_FixedHashTable.hpp>
00053 #include <Tpetra_HashTable.hpp>
00054 
00055 
00056 // FIXME (mfh 16 Apr 2013) GIANT HACK BELOW
00057 #ifdef HAVE_MPI
00058 #  include "mpi.h"
00059 #endif // HAVE_MPI
00060 // FIXME (mfh 16 Apr 2013) GIANT HACK ABOVE
00061 
00062 
00063 namespace Tpetra {
00064   namespace Details {
00065     template<class LO, class GO, class NT>
00066     Directory<LO, GO, NT>::
00067     Directory () {}
00068 
00069     template<class LO, class GO, class NT>
00070     LookupStatus
00071     Directory<LO, GO, NT>::
00072     getEntries (const map_type& map,
00073                 const Teuchos::ArrayView<const GO> &globalIDs,
00074                 const Teuchos::ArrayView<int> &nodeIDs,
00075                 const Teuchos::ArrayView<LO> &localIDs,
00076                 const bool computeLIDs) const
00077     {
00078       // Ensure that globalIDs, nodeIDs, and localIDs (if applicable)
00079       // all have the same size, before modifying any output arguments.
00080       TEUCHOS_TEST_FOR_EXCEPTION(nodeIDs.size() != globalIDs.size(),
00081         std::invalid_argument, Teuchos::typeName(*this) << "::getEntries(): "
00082         "Output arrays do not have the right sizes.  nodeIDs.size() = "
00083         << nodeIDs.size() << " != globalIDs.size() = " << globalIDs.size()
00084         << ".");
00085       TEUCHOS_TEST_FOR_EXCEPTION(
00086         computeLIDs && localIDs.size() != globalIDs.size(),
00087         std::invalid_argument, Teuchos::typeName(*this) << "::getEntries(): "
00088         "Output array do not have the right sizes.  localIDs.size() = "
00089         << localIDs.size() << " != globalIDs.size() = " << globalIDs.size()
00090         << ".");
00091 
00092       // Initially, fill nodeIDs and localIDs (if applicable) with
00093       // invalid values.  The "invalid" process ID is -1 (this means
00094       // the same thing as MPI_ANY_SOURCE to Teuchos, so it's an
00095       // "invalid" process ID); the invalid local ID comes from
00096       // OrdinalTraits.
00097       std::fill (nodeIDs.begin(), nodeIDs.end(), -1);
00098       if (computeLIDs) {
00099         std::fill (localIDs.begin(), localIDs.end(),
00100                    Teuchos::OrdinalTraits<LO>::invalid ());
00101       }
00102       // Actually do the work.
00103       return this->getEntriesImpl (map, globalIDs, nodeIDs, localIDs, computeLIDs);
00104     }
00105 
00106 
00107     template<class LO, class GO, class NT>
00108     ReplicatedDirectory<LO, GO, NT>::
00109     ReplicatedDirectory (const map_type& map) :
00110       numProcs_ (map.getComm ()->getSize ())
00111     {}
00112 
00113 
00114     template<class LO, class GO, class NT>
00115     ReplicatedDirectory<LO, GO, NT>::
00116     ReplicatedDirectory () :
00117       numProcs_ (0) // to be set later
00118     {}
00119 
00120 
00121     template<class LO, class GO, class NT>
00122     bool
00123     ReplicatedDirectory<LO, GO, NT>::
00124     isOneToOne (const Teuchos::Comm<int>& comm) const
00125     {
00126       // A locally replicated Map is one-to-one only if there is no
00127       // replication, that is, only if the Map's communicator only has
00128       // one process.
00129       return (numProcs_ == 1);
00130     }
00131 
00132 
00133     template<class LO, class GO, class NT>
00134     std::string
00135     ReplicatedDirectory<LO, GO, NT>::description () const
00136     {
00137       std::ostringstream os;
00138       os << "ReplicatedDirectory"
00139          << "<" << Teuchos::TypeNameTraits<LO>::name ()
00140          << ", " << Teuchos::TypeNameTraits<GO>::name ()
00141          << ", " << Teuchos::TypeNameTraits<NT>::name () << ">";
00142       return os.str ();
00143     }
00144 
00145 
00146     template<class LO, class GO, class NT>
00147     ContiguousUniformDirectory<LO, GO, NT>::
00148     ContiguousUniformDirectory (const map_type& map)
00149     {
00150       TEUCHOS_TEST_FOR_EXCEPTION(! map.isContiguous (), std::invalid_argument,
00151         Teuchos::typeName (*this) << " constructor: Map is not contiguous.");
00152       TEUCHOS_TEST_FOR_EXCEPTION(! map.isUniform (), std::invalid_argument,
00153         Teuchos::typeName (*this) << " constructor: Map is not uniform.");
00154     }
00155 
00156 
00157     template<class LO, class GO, class NT>
00158     std::string
00159     ContiguousUniformDirectory<LO, GO, NT>::description () const
00160     {
00161       std::ostringstream os;
00162       os << "ContiguousUniformDirectory"
00163          << "<" << Teuchos::TypeNameTraits<LO>::name ()
00164          << ", " << Teuchos::TypeNameTraits<GO>::name ()
00165          << ", " << Teuchos::TypeNameTraits<NT>::name () << ">";
00166       return os.str ();
00167     }
00168 
00169 
00170     template<class LO, class GO, class NT>
00171     LookupStatus
00172     ContiguousUniformDirectory<LO, GO, NT>::
00173     getEntriesImpl (const map_type& map,
00174                     const Teuchos::ArrayView<const GO> &globalIDs,
00175                     const Teuchos::ArrayView<int> &nodeIDs,
00176                     const Teuchos::ArrayView<LO> &localIDs,
00177                     const bool computeLIDs) const
00178     {
00179       using Teuchos::as;
00180       using Teuchos::Comm;
00181       using Teuchos::RCP;
00182       typedef typename Teuchos::ArrayView<const GO>::size_type size_type;
00183       const LO invalidLid = Teuchos::OrdinalTraits<LO>::invalid ();
00184       LookupStatus res = AllIDsPresent;
00185 
00186       RCP<const Comm<int> > comm = map.getComm ();
00187       const GO g_min = map.getMinAllGlobalIndex ();
00188 
00189       // Let N_G be the global number of elements in the Map,
00190       // and P be the number of processes in its communicator.
00191       // Then, N_G = P * N_L + R = R*(N_L + 1) + (P - R)*N_L.
00192       //
00193       // The first R processes own N_L+1 elements.
00194       // The remaining P-R processes own N_L elements.
00195       //
00196       // Let g be the current GID, g_min be the global minimum GID,
00197       // and g_0 = g - g_min.  If g is a valid GID in this Map, then
00198       // g_0 is in [0, N_G - 1].
00199       //
00200       // If g is a valid GID in this Map and g_0 < R*(N_L + 1), then
00201       // the rank of the process that owns g is floor(g_0 / (N_L +
00202       // 1)), and its corresponding local index on that process is g_0
00203       // mod (N_L + 1).
00204       //
00205       // Let g_R = g_0 - R*(N_L + 1).  If g is a valid GID in this Map
00206       // and g_0 >= R*(N_L + 1), then the rank of the process that
00207       // owns g is then R + floor(g_R / N_L), and its corresponding
00208       // local index on that process is g_R mod N_L.
00209 
00210       const size_type N_G = as<size_type> (map.getGlobalNumElements ());
00211       const size_type P = as<size_type> (comm->getSize ());
00212       const size_type N_L  = N_G / P;
00213       const size_type R = N_G - N_L * P; // N_G mod P
00214       const size_type N_R = R * (N_L + 1);
00215 
00216 #ifdef HAVE_TPETRA_DEBUG
00217       TEUCHOS_TEST_FOR_EXCEPTION(
00218         N_G != P*N_L + R, std::logic_error,
00219         "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
00220         "N_G = " << N_G << " != P*N_L + R = " << P << "*" << N_L << " + " << R
00221         << " = " << P*N_L + R << ".  "
00222         "Please report this bug to the Tpetra developers.");
00223 #endif // HAVE_TPETRA_DEBUG
00224 
00225       const size_type numGids = globalIDs.size (); // for const loop bound
00226       if (computeLIDs) {
00227         for (size_type k = 0; k < numGids; ++k) {
00228           const GO g_0 = globalIDs[k] - g_min;
00229           // The first test is a little strange just in case GO is
00230           // unsigned.  Compilers raise a warning on tests like "x <
00231           // 0" if x is unsigned, but don't usually raise a warning if
00232           // the expression is a bit more complicated than that.
00233           if (g_0 + 1 < 1 || g_0 >= N_G) {
00234             nodeIDs[k] = -1;
00235             localIDs[k] = invalidLid;
00236             res = IDNotPresent;
00237           }
00238           else if (g_0 < as<GO> (N_R)) {
00239             // The GID comes from the initial sequence of R processes.
00240             nodeIDs[k] = as<int> (g_0 / (N_L + 1));
00241             localIDs[k] = g_0 % (N_L + 1);
00242           }
00243           else if (g_0 >= as<GO> (N_R)) {
00244             // The GID comes from the remaining P-R processes.
00245             const GO g_R = g_0 - N_R;
00246             nodeIDs[k] = as<int> (R + g_R / N_L);
00247             localIDs[k] = as<int> (g_R % N_L);
00248           }
00249 #ifdef HAVE_TPETRA_DEBUG
00250           else {
00251             TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00252               "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
00253               "should never get here.  "
00254               "Please report this bug to the Tpetra developers.");
00255           }
00256 #endif // HAVE_TPETRA_DEBUG
00257         }
00258       }
00259       else { // don't compute local indices
00260         for (size_type k = 0; k < numGids; ++k) {
00261           const GO g_0 = globalIDs[k] - g_min;
00262           // The first test is a little strange just in case GO is
00263           // unsigned.  Compilers raise a warning on tests like "x <
00264           // 0" if x is unsigned, but don't usually raise a warning if
00265           // the expression is a bit more complicated than that.
00266           if (g_0 + 1 < 1 || g_0 >= N_G) {
00267             nodeIDs[k] = -1;
00268             res = IDNotPresent;
00269           }
00270           else if (g_0 < as<GO> (N_R)) {
00271             // The GID comes from the initial sequence of R processes.
00272             nodeIDs[k] = as<int> (g_0 / (N_L + 1));
00273           }
00274           else if (g_0 >= as<GO> (N_R)) {
00275             // The GID comes from the remaining P-R processes.
00276             const GO g_R = g_0 - N_R;
00277             nodeIDs[k] = as<int> (R + g_R / N_L);
00278           }
00279 #ifdef HAVE_TPETRA_DEBUG
00280           else {
00281             TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00282               "Tpetra::ContiguousUniformDirectory::getEntriesImpl: "
00283               "should never get here.  "
00284               "Please report this bug to the Tpetra developers.");
00285           }
00286 #endif // HAVE_TPETRA_DEBUG
00287         }
00288       }
00289       return res;
00290     }
00291 
00292     template<class LO, class GO, class NT>
00293     DistributedContiguousDirectory<LO, GO, NT>::
00294     DistributedContiguousDirectory (const map_type& map)
00295     {
00296       using Teuchos::gatherAll;
00297       using Teuchos::RCP;
00298 
00299       RCP<const Teuchos::Comm<int> > comm = map.getComm ();
00300 
00301       TEUCHOS_TEST_FOR_EXCEPTION(! map.isDistributed (), std::invalid_argument,
00302         Teuchos::typeName (*this) << " constructor: Map is not distributed.");
00303       TEUCHOS_TEST_FOR_EXCEPTION(! map.isContiguous (), std::invalid_argument,
00304         Teuchos::typeName (*this) << " constructor: Map is not contiguous.");
00305 
00306       const int numProcs = comm->getSize ();
00307 
00308       // Make room for the min global ID on each process, plus one
00309       // entry at the end for the "max cap."
00310       allMinGIDs_ = arcp<GO> (numProcs + 1);
00311       // Get my process' min global ID.
00312       GO minMyGID = map.getMinGlobalIndex ();
00313       // Gather all of the min global IDs into the first numProcs
00314       // entries of allMinGIDs_.
00315 
00316       // FIXME (mfh 16 Apr 2013) GIANT HACK BELOW
00317       //
00318       // The purpose of this giant hack is that gatherAll appears to
00319       // interpret the "receive count" argument differently than
00320       // MPI_Allgather does.  Matt Bettencourt reports Valgrind issues
00321       // (memcpy with overlapping data) with MpiComm<int>::gatherAll,
00322       // which could relate either to this, or to OpenMPI.
00323 #ifdef HAVE_MPI
00324       MPI_Datatype rawMpiType = MPI_INT;
00325       bool useRawMpi = true;
00326       if (typeid (GO) == typeid (int)) {
00327         rawMpiType = MPI_INT;
00328       } else if (typeid (GO) == typeid (long)) {
00329         rawMpiType = MPI_LONG;
00330       } else {
00331         useRawMpi = false;
00332       }
00333       if (useRawMpi) {
00334         using Teuchos::rcp_dynamic_cast;
00335         using Teuchos::MpiComm;
00336         RCP<const MpiComm<int> > mpiComm =
00337           rcp_dynamic_cast<const MpiComm<int> > (comm);
00338         // It could be a SerialComm instead, even in an MPI build, so
00339         // be sure to check.
00340         if (! comm.is_null ()) {
00341           MPI_Comm rawMpiComm = * (mpiComm->getRawMpiComm ());
00342           const int err =
00343             MPI_Allgather (&minMyGID, 1, rawMpiType,
00344                            allMinGIDs_.getRawPtr (), 1, rawMpiType,
00345                            rawMpiComm);
00346           TEUCHOS_TEST_FOR_EXCEPTION(err != MPI_SUCCESS, std::runtime_error,
00347             "Tpetra::DistributedContiguousDirectory: MPI_Allgather failed");
00348         } else {
00349           gatherAll<int, GO> (*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr ());
00350         }
00351       } else {
00352         gatherAll<int, GO> (*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr ());
00353       }
00354 #else // NOT HAVE_MPI
00355       gatherAll<int, GO> (*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr ());
00356 #endif // HAVE_MPI
00357       // FIXME (mfh 16 Apr 2013) GIANT HACK ABOVE
00358 
00359       //gatherAll<int, GO> (*comm, 1, &minMyGID, numProcs, allMinGIDs_.getRawPtr ());
00360 
00361       // Put the max cap at the end.  Adding one lets us write loops
00362       // over the global IDs with the usual strict less-than bound.
00363       allMinGIDs_[numProcs] = map.getMaxAllGlobalIndex ()
00364         + Teuchos::OrdinalTraits<GO>::one ();
00365     }
00366 
00367     template<class LO, class GO, class NT>
00368     std::string
00369     DistributedContiguousDirectory<LO, GO, NT>::description () const
00370     {
00371       std::ostringstream os;
00372       os << "DistributedContiguousDirectory"
00373          << "<" << Teuchos::TypeNameTraits<LO>::name ()
00374          << ", " << Teuchos::TypeNameTraits<GO>::name ()
00375          << ", " << Teuchos::TypeNameTraits<NT>::name () << ">";
00376       return os.str ();
00377     }
00378 
00379     template<class LO, class GO, class NT>
00380     LookupStatus
00381     ReplicatedDirectory<LO, GO, NT>::
00382     getEntriesImpl (const map_type& map,
00383                     const Teuchos::ArrayView<const GO> &globalIDs,
00384                     const Teuchos::ArrayView<int> &nodeIDs,
00385                     const Teuchos::ArrayView<LO> &localIDs,
00386                     const bool computeLIDs) const
00387     {
00388       using Teuchos::Array;
00389       using Teuchos::ArrayRCP;
00390       using Teuchos::ArrayView;
00391       using Teuchos::as;
00392       using Teuchos::Comm;
00393       using Teuchos::RCP;
00394 
00395       LookupStatus res = AllIDsPresent;
00396       RCP<const Teuchos::Comm<int> > comm = map.getComm ();
00397       const int myRank = comm->getRank ();
00398 
00399       // Map is on one process or is locally replicated.
00400       typename ArrayView<int>::iterator procIter = nodeIDs.begin();
00401       typename ArrayView<LO>::iterator lidIter = localIDs.begin();
00402       typename ArrayView<const GO>::iterator gidIter;
00403       for (gidIter = globalIDs.begin(); gidIter != globalIDs.end(); ++gidIter) {
00404         if (map.isNodeGlobalElement (*gidIter)) {
00405           *procIter++ = myRank;
00406           if (computeLIDs) {
00407             *lidIter++ = map.getLocalElement (*gidIter);
00408           }
00409         }
00410         else {
00411           // Advance the pointers, leaving these values set to invalid
00412           procIter++;
00413           if (computeLIDs) {
00414             lidIter++;
00415           }
00416           res = IDNotPresent;
00417         }
00418       }
00419       return res;
00420     }
00421 
00422     template<class LO, class GO, class NT>
00423     LookupStatus
00424     DistributedContiguousDirectory<LO, GO, NT>::
00425     getEntriesImpl (const map_type& map,
00426                     const Teuchos::ArrayView<const GO> &globalIDs,
00427                     const Teuchos::ArrayView<int> &nodeIDs,
00428                     const Teuchos::ArrayView<LO> &localIDs,
00429                     const bool computeLIDs) const
00430     {
00431       using Teuchos::Array;
00432       using Teuchos::ArrayRCP;
00433       using Teuchos::ArrayView;
00434       using Teuchos::as;
00435       using Teuchos::Comm;
00436       using Teuchos::RCP;
00437 
00438       RCP<const Teuchos::Comm<int> > comm = map.getComm ();
00439       const int numProcs = comm->getSize ();
00440       const global_size_t nOverP = map.getGlobalNumElements () / numProcs;
00441       const LO LINVALID = Teuchos::OrdinalTraits<LO>::invalid();
00442       LookupStatus res = AllIDsPresent;
00443 
00444       // Map is distributed but contiguous.
00445       typename ArrayView<int>::iterator procIter = nodeIDs.begin();
00446       typename ArrayView<LO>::iterator lidIter = localIDs.begin();
00447       typename ArrayView<const GO>::iterator gidIter;
00448       for (gidIter = globalIDs.begin(); gidIter != globalIDs.end(); ++gidIter) {
00449         LO LID = LINVALID; // Assume not found until proven otherwise
00450         int image = -1;
00451         GO GID = *gidIter;
00452         // Guess uniform distribution and start a little above it
00453         // TODO: replace by a binary search
00454         int curRank;
00455         { // We go through all this trouble to avoid overflow and
00456           // signed / unsigned casting mistakes (that were made in
00457           // previous versions of this code).
00458           const GO one = as<GO> (1);
00459           const GO two = as<GO> (2);
00460           const GO nOverP_GID = as<GO> (nOverP);
00461           const GO lowerBound = GID / std::max(nOverP_GID, one) + two;
00462           curRank = as<int>(std::min(lowerBound, as<GO>(numProcs - 1)));
00463         }
00464         bool found = false;
00465         while (curRank >= 0 && curRank < numProcs) {
00466           if (allMinGIDs_[curRank] <= GID) {
00467             if (GID < allMinGIDs_[curRank + 1]) {
00468               found = true;
00469               break;
00470             }
00471             else {
00472               curRank++;
00473             }
00474           }
00475           else {
00476             curRank--;
00477           }
00478         }
00479         if (found) {
00480           image = curRank;
00481           LID = as<LO> (GID - allMinGIDs_[image]);
00482         }
00483         else {
00484           res = IDNotPresent;
00485         }
00486         *procIter++ = image;
00487         if (computeLIDs) {
00488           *lidIter++ = LID;
00489         }
00490       }
00491       return res;
00492     }
00493 
00494     template<class LO, class GO, class NT>
00495     DistributedNoncontiguousDirectory<LO, GO, NT>::
00496     DistributedNoncontiguousDirectory (const map_type& map) :
00497       oneToOneResult_ (ONE_TO_ONE_NOT_CALLED_YET), // to be revised below
00498       locallyOneToOne_ (true), // to be revised below
00499       useHashTables_ (false) // to be revised below
00500     {
00501       initialize (map, Teuchos::null);
00502     }
00503 
00504     template<class LO, class GO, class NT>
00505     DistributedNoncontiguousDirectory<LO, GO, NT>::
00506     DistributedNoncontiguousDirectory (const map_type& map,
00507                                        const tie_break_type& tie_break) :
00508       oneToOneResult_ (ONE_TO_ONE_NOT_CALLED_YET), // to be revised below
00509       locallyOneToOne_ (true), // to be revised below
00510       useHashTables_ (false) // to be revised below
00511     {
00512       initialize (map, Teuchos::ptrFromRef (tie_break));
00513     }
00514 
00515     template<class LO, class GO, class NT>
00516     void
00517     DistributedNoncontiguousDirectory<LO, GO, NT>::
00518     initialize (const map_type& map,
00519                 Teuchos::Ptr<const tie_break_type> tie_break)
00520     {
00521       using Teuchos::Array;
00522       using Teuchos::ArrayView;
00523       using Teuchos::as;
00524       using Teuchos::RCP;
00525       using Teuchos::rcp;
00526       using Teuchos::typeName;
00527       using Teuchos::TypeNameTraits;
00528       using std::cerr;
00529       using std::endl;
00530       typedef Array<int>::size_type size_type;
00531 
00532       // This class' implementation of getEntriesImpl() currently
00533       // encodes the following assumptions:
00534       //
00535       // 1. global_size_t >= GO
00536       // 2. global_size_t >= int
00537       // 3. global_size_t >= LO
00538       //
00539       // We check these assumptions here.
00540       TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(GO),
00541         std::logic_error, typeName (*this) << ": sizeof(Tpetra::"
00542         "global_size_t) = " << sizeof(global_size_t) << " < sizeof(Global"
00543         "Ordinal = " << TypeNameTraits<LO>::name () << ") = " << sizeof(GO)
00544         << ".");
00545       TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(int),
00546         std::logic_error, typeName (*this) << ": sizeof(Tpetra::"
00547         "global_size_t) = " << sizeof(global_size_t) << " < sizeof(int) = "
00548         << sizeof(int) << ".");
00549       TEUCHOS_TEST_FOR_EXCEPTION(sizeof(global_size_t) < sizeof(LO),
00550         std::logic_error, typeName (*this) << ": sizeof(Tpetra::"
00551         "global_size_t) = " << sizeof(global_size_t) << " < sizeof(Local"
00552         "Ordinal = " << TypeNameTraits<LO>::name () << ") = " << sizeof(LO)
00553         << ".");
00554 
00555       RCP<const Teuchos::Comm<int> > comm = map.getComm ();
00556       const LO LINVALID = Teuchos::OrdinalTraits<LO>::invalid ();
00557       const GO minAllGID = map.getMinAllGlobalIndex ();
00558       const GO maxAllGID = map.getMaxAllGlobalIndex ();
00559 
00560       // The "Directory Map" (see below) will have a range of elements
00561       // from the minimum to the maximum GID of the user Map, and a
00562       // minimum GID of minAllGID from the user Map.  It doesn't
00563       // actually have to store all those entries, though do beware of
00564       // calling getNodeElementList on it (see Bug 5822).
00565       const global_size_t numGlobalEntries = maxAllGID - minAllGID + 1;
00566 
00567       // We can't afford to replicate the whole directory on each
00568       // process, so create the "Directory Map", a uniform contiguous
00569       // Map that describes how we will distribute the directory over
00570       // processes.
00571       //
00572       // FIXME (mfh 08 May 2012) Here we're setting minAllGID to be
00573       // the index base.  The index base should be separate from the
00574       // minimum GID.
00575       directoryMap_ = rcp (new map_type (numGlobalEntries, minAllGID, comm,
00576                                          GloballyDistributed, map.getNode ()));
00577       // The number of Directory elements that my process owns.
00578       const size_t dir_numMyEntries = directoryMap_->getNodeNumElements ();
00579 
00580       // Fix for Bug 5822: If the input Map is "sparse," that is if
00581       // the difference between the global min and global max GID is
00582       // much larger than the global number of elements in the input
00583       // Map, then it's possible that the Directory Map might have
00584       // many more entries than the input Map on this process.  This
00585       // can cause memory scalability issues.  In that case, we switch
00586       // from the array-based implementation of Directory storage to
00587       // the hash table - based implementation.  We don't use hash
00588       // tables all the time, because they are slower in the common
00589       // case of a nonsparse Map.
00590       //
00591       // NOTE: This is a per-process decision.  Some processes may use
00592       // array-based storage, whereas others may use hash table -
00593       // based storage.
00594 
00595       // A hash table takes a constant factor more space, more or
00596       // less, than an array.  Thus, it's not worthwhile, even in
00597       // terms of memory usage, always to use a hash table.
00598       // Furthermore, array lookups are faster than hash table
00599       // lookups, so it may be worthwhile to use an array even if it
00600       // takes more space.  The "sparsity threshold" governs when to
00601       // switch to a hash table - based implementation.
00602       const size_t inverseSparsityThreshold = 10;
00603       useHashTables_ =
00604         dir_numMyEntries >= inverseSparsityThreshold * map.getNodeNumElements ();
00605 
00606       // Get list of process IDs that own the directory entries for the
00607       // Map GIDs.  These will be the targets of the sends that the
00608       // Distributor will do.
00609       const int myRank = comm->getRank ();
00610       const size_t numMyEntries = map.getNodeNumElements ();
00611       Array<int> sendImageIDs (numMyEntries);
00612       ArrayView<const GO> myGlobalEntries = map.getNodeElementList ();
00613       // An ID not present in this lookup indicates that it lies outside
00614       // of the range [minAllGID,maxAllGID] (from map_).  this means
00615       // something is wrong with map_, our fault.
00616       const LookupStatus lookupStatus =
00617         directoryMap_->getRemoteIndexList (myGlobalEntries, sendImageIDs);
00618       TEUCHOS_TEST_FOR_EXCEPTION(
00619         lookupStatus == IDNotPresent, std::logic_error, Teuchos::typeName(*this)
00620         << " constructor: the Directory Map could not find out where one or "
00621         "more of my Map's indices should go.  The input to getRemoteIndexList "
00622         "is " << Teuchos::toString (myGlobalEntries) << ", and the output is "
00623         << Teuchos::toString (sendImageIDs ()) << ".  The input Map itself has "
00624         "the following entries on the calling process " <<
00625         map.getComm ()->getRank () << ": " <<
00626         Teuchos::toString (map.getNodeElementList ()) << ", and has "
00627         << map.getGlobalNumElements () << " total global indices in ["
00628         << map.getMinAllGlobalIndex () << "," << map.getMaxAllGlobalIndex ()
00629         << "].  The Directory Map has "
00630         << directoryMap_->getGlobalNumElements () << " total global indices in "
00631         "[" << directoryMap_->getMinAllGlobalIndex () << "," <<
00632         directoryMap_->getMaxAllGlobalIndex () << "], and the calling process "
00633         "has GIDs [" << directoryMap_->getMinGlobalIndex () << "," <<
00634         directoryMap_->getMaxGlobalIndex () << "].  "
00635         "This probably means there is a bug in Map or Directory.  "
00636         "Please report this bug to the Tpetra developers.");
00637 
00638       // Initialize the distributor using the list of process IDs to
00639       // which to send.  We'll use the distributor to send out triples
00640       // of (GID, process ID, LID).  We're sending the entries to the
00641       // processes that the Directory Map says should own them, which is
00642       // why we called directoryMap_->getRemoteIndexList() above.
00643       Distributor distor (comm);
00644       const size_t numReceives = distor.createFromSends (sendImageIDs);
00645 
00646       // NOTE (mfh 21 Mar 2012) The following code assumes that
00647       // sizeof(GO) >= sizeof(int) and sizeof(GO) >= sizeof(LO).
00648       //
00649       // Create and fill buffer of (GID, PID, LID) triples to send
00650       // out.  We pack the (GID, PID, LID) triples into a single Array
00651       // of GO, casting the PID from int to GO and the LID from LO to
00652       // GO as we do so.
00653       //
00654       // FIXME (mfh 23 Mar 2014) This assumes that sizeof(LO) <=
00655       // sizeof(GO) and sizeof(int) <= sizeof(GO).  The former is
00656       // required, and the latter is generally the case, but we should
00657       // still check for this.
00658       const int packetSize = 3; // We're sending triples, so packet size is 3.
00659       Array<GO> exportEntries (packetSize * numMyEntries); // data to send out
00660       {
00661         size_type exportIndex = 0;
00662         for (size_type i = 0; i < static_cast<size_type> (numMyEntries); ++i) {
00663           exportEntries[exportIndex++] = myGlobalEntries[i];
00664           exportEntries[exportIndex++] = as<GO> (myRank);
00665           exportEntries[exportIndex++] = as<GO> (i);
00666         }
00667       }
00668       // Buffer of data to receive.  The Distributor figured out for
00669       // us how many packets we're receiving, when we called its
00670       // createFromSends() method to set up the distribution plan.
00671       Array<GO> importElements (packetSize * distor.getTotalReceiveLength ());
00672 
00673       // Distribute the triples of (GID, process ID, LID).
00674       distor.doPostsAndWaits (exportEntries ().getConst (), packetSize, importElements ());
00675 
00676       // Unpack the redistributed data.  Both implementations of
00677       // Directory storage map from an LID in the Directory Map (which
00678       // is the LID of the GID to store) to either a PID or an LID in
00679       // the input Map.  Each "packet" (contiguous chunk of
00680       // importElements) contains a triple: (GID, PID, LID).
00681       if (useHashTables_) {
00682         // Create the hash tables.  We know exactly how many elements
00683         // to expect in each hash table.  FixedHashTable's constructor
00684         // currently requires all the keys and values at once, so we
00685         // have to extract them in temporary arrays.  It may be
00686         // possible to rewrite FixedHashTable to use a "start fill" /
00687         // "end fill" approach that avoids the temporary arrays, but
00688         // we won't try that for now.
00689 
00690         // The constructors of Array and ArrayRCP that take a number
00691         // of elements all initialize the arrays.  Instead, allocate
00692         // raw arrays, then hand them off to ArrayRCP, to avoid the
00693         // initial unnecessary initialization without losing the
00694         // benefit of exception safety (and bounds checking, in a
00695         // debug build).
00696         LO* tableKeysRaw = NULL;
00697         LO* tableLidsRaw = NULL;
00698         int* tablePidsRaw = NULL;
00699         try {
00700           tableKeysRaw = new LO [numReceives];
00701           tableLidsRaw = new LO [numReceives];
00702           tablePidsRaw = new int [numReceives];
00703         } catch (...) {
00704           if (tableKeysRaw != NULL) {
00705             delete [] tableKeysRaw;
00706           }
00707           if (tableLidsRaw != NULL) {
00708             delete [] tableLidsRaw;
00709           }
00710           if (tablePidsRaw != NULL) {
00711             delete [] tablePidsRaw;
00712           }
00713           throw;
00714         }
00715         ArrayRCP<LO> tableKeys (tableKeysRaw, 0, numReceives, true);
00716         ArrayRCP<LO> tableLids (tableLidsRaw, 0, numReceives, true);
00717         ArrayRCP<int> tablePids (tablePidsRaw, 0, numReceives, true);
00718 
00719         if (tie_break.is_null ()) {
00720           // Fill the temporary arrays of keys and values.
00721           size_type importIndex = 0;
00722           for (size_type i = 0; i < static_cast<size_type> (numReceives); ++i) {
00723             const GO curGID = importElements[importIndex++];
00724             const LO curLID = directoryMap_->getLocalElement (curGID);
00725             TEUCHOS_TEST_FOR_EXCEPTION(
00726               curLID == LINVALID, std::logic_error,
00727               Teuchos::typeName(*this) << " constructor: Incoming global index "
00728               << curGID << " does not have a corresponding local index in the "
00729               "Directory Map.  Please report this bug to the Tpetra developers.");
00730             tableKeys[i] = curLID;
00731             tablePids[i] = importElements[importIndex++];
00732             tableLids[i] = importElements[importIndex++];
00733           }
00734           // Set up the hash tables.  The hash tables' constructor
00735           // detects whether there are duplicates, so that we can set
00736           // locallyOneToOne_.
00737           lidToPidTable_ =
00738             rcp (new Details::FixedHashTable<LO, int> (tableKeys (),
00739                                                        tablePids ()));
00740           locallyOneToOne_ = ! (lidToPidTable_->hasDuplicateKeys ());
00741           lidToLidTable_ =
00742             rcp (new Details::FixedHashTable<LO, LO> (tableKeys (),
00743                                                       tableLids ()));
00744         }
00745         else { // tie_break is NOT null
00746 
00747           // For each directory Map LID received, collect all the
00748           // corresponding (PID,LID) pairs.  If the input Map is not
00749           // one-to-one, corresponding directory Map LIDs will have
00750           // more than one pair.  In that case, we will use the
00751           // TieBreak object to pick exactly one pair.
00752           typedef std::map<LO, std::vector<std::pair<int, LO> > > pair_table_type;
00753           pair_table_type ownedPidLidPairs;
00754 
00755           // For each directory Map LID received, collect the zero or
00756           // more input Map (PID,LID) pairs into ownedPidLidPairs.
00757           size_type importIndex = 0;
00758           for (size_type i = 0; i < static_cast<size_type> (numReceives); ++i) {
00759             const GO curGID = importElements[importIndex++];
00760             const LO dirMapLid = directoryMap_->getLocalElement (curGID);
00761             TEUCHOS_TEST_FOR_EXCEPTION(
00762               dirMapLid == LINVALID, std::logic_error,
00763               Teuchos::typeName(*this) << " constructor: Incoming global index "
00764               << curGID << " does not have a corresponding local index in the "
00765               "Directory Map.  Please report this bug to the Tpetra developers.");
00766             tableKeys[i] = dirMapLid;
00767             const int PID = importElements[importIndex++];
00768             const int LID = importElements[importIndex++];
00769 
00770             // These may change below.  We fill them in just to ensure
00771             // that they won't have invalid values.
00772             tablePids[i] = PID;
00773             tableLids[i] = LID;
00774 
00775             // For every directory Map LID, we have to remember all
00776             // (PID, LID) pairs.  The TieBreak object will arbitrate
00777             // between them in the loop below.
00778             ownedPidLidPairs[dirMapLid].push_back (std::make_pair (PID, LID));
00779           }
00780 
00781           // Use TieBreak to arbitrate between (PID,LID) pairs
00782           // corresponding to each directory Map LID.
00783           //
00784           // FIXME (mfh 23 Mar 2014) How do I know that i is the same
00785           // as the directory Map LID?
00786           const size_type numPairs =
00787             static_cast<size_type> (ownedPidLidPairs.size ());
00788           for (size_type i = 0; i < numPairs; ++i) {
00789             const LO dirMapLid = static_cast<LO> (i);
00790             const GO dirMapGid = directoryMap_->getGlobalElement (dirMapLid);
00791             const std::vector<std::pair<int, LO> >& pidLidList =
00792               ownedPidLidPairs[i];
00793             const size_t listLen = pidLidList.size ();
00794             if (listLen > 0) {
00795               if (listLen > 1) {
00796                 locallyOneToOne_ = false;
00797               }
00798               // If there is some (PID,LID) pair for the current input
00799               // Map LID, then it makes sense to invoke the TieBreak
00800               // object to arbitrate between the options.  Even if
00801               // there is only one (PID,LID) pair, we still want to
00802               // give the TieBreak object a chance to do whatever it
00803               // likes to do, in terms of side effects (e.g., track
00804               // (PID,LID) pairs).
00805               const size_type index =
00806                 static_cast<size_type> (tie_break->selectedIndex (dirMapGid,
00807                                                                   pidLidList));
00808               tablePids[i] = pidLidList[index].first;
00809               tableLids[i] = pidLidList[index].second;
00810             }
00811           }
00812 
00813           // Set up the hash tables.
00814           lidToPidTable_ =
00815             rcp (new Details::FixedHashTable<LO, int> (tableKeys (),
00816                                                        tablePids ()));
00817           lidToLidTable_ =
00818             rcp (new Details::FixedHashTable<LO, LO> (tableKeys (),
00819                                                       tableLids ()));
00820         }
00821       }
00822       else {
00823         if (tie_break.is_null ()) {
00824           // Use array-based implementation of Directory storage.
00825           // Allocate these arrays and fill them with invalid values,
00826           // in case the input Map's GID list is sparse (i.e., does
00827           // not populate all GIDs from minAllGID to maxAllGID).
00828           PIDs_ = arcp<int> (dir_numMyEntries);
00829           std::fill (PIDs_.begin (), PIDs_.end (), -1);
00830           LIDs_ = arcp<LO> (dir_numMyEntries);
00831           std::fill (LIDs_.begin (), LIDs_.end (), LINVALID);
00832           // Fill in the arrays with PIDs resp. LIDs.
00833           size_type importIndex = 0;
00834           for (size_type i = 0; i < static_cast<size_type> (numReceives); ++i) {
00835             const GO curGID = importElements[importIndex++];
00836             const LO curLID = directoryMap_->getLocalElement (curGID);
00837             TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
00838               Teuchos::typeName(*this) << " constructor: Incoming global index "
00839               << curGID << " does not have a corresponding local index in the "
00840               "Directory Map.  Please report this bug to the Tpetra developers.");
00841 
00842             // If PIDs_[curLID] is not -1, then curGID is a duplicate
00843             // on the calling process, so the Directory is not locally
00844             // one-to-one.
00845             if (PIDs_[curLID] != -1) {
00846               locallyOneToOne_ = false;
00847             }
00848             PIDs_[curLID] = importElements[importIndex++];
00849             LIDs_[curLID] = importElements[importIndex++];
00850           }
00851         }
00852         else {
00853           PIDs_ = arcp<int> (dir_numMyEntries);
00854           LIDs_ = arcp<LO> (dir_numMyEntries);
00855           std::fill (PIDs_.begin (), PIDs_.end (), -1);
00856 
00857           // All received (PID, LID) pairs go into ownedPidLidPairs.
00858           // This is a map from the directory Map's LID to the (PID,
00859           // LID) pair (where the latter LID comes from the input Map,
00860           // not the directory Map).  If the input Map is not
00861           // one-to-one, corresponding LIDs will have
00862           // ownedPidLidPairs[curLID].size() > 1.  In that case, we
00863           // will use the TieBreak object to pick exactly one pair.
00864           Array<std::vector<std::pair<int, LO> > > ownedPidLidPairs (dir_numMyEntries);
00865           size_type importIndex = 0;
00866           for (size_type i = 0; i < static_cast<size_type> (numReceives); ++i) {
00867             const GO  GID = importElements[importIndex++];
00868             const int PID = importElements[importIndex++];
00869             const LO  LID = importElements[importIndex++];
00870 
00871             const LO dirMapLid = directoryMap_->getLocalElement (GID);
00872             TEUCHOS_TEST_FOR_EXCEPTION(
00873               dirMapLid == LINVALID, std::logic_error,
00874               Teuchos::typeName(*this) << " constructor: Incoming global index "
00875               << GID << " does not have a corresponding local index in the "
00876               "Directory Map.  Please report this bug to the Tpetra developers.");
00877             ownedPidLidPairs[dirMapLid].push_back (std::make_pair (PID, LID));
00878           }
00879 
00880           // Use TieBreak to arbitrate between (PID,LID) pairs
00881           // corresponding to each directory Map LID.
00882           //
00883           // FIXME (mfh 23 Mar 2014) How do I know that i is the same
00884           // as the directory Map LID?
00885           const size_type numPairs =
00886             static_cast<size_type> (ownedPidLidPairs.size ());
00887           for (size_type i = 0; i < numPairs; ++i) {
00888             const LO dirMapLid = static_cast<LO> (i);
00889             const GO dirMapGid = directoryMap_->getGlobalElement (dirMapLid);
00890             const std::vector<std::pair<int, LO> >& pidLidList =
00891               ownedPidLidPairs[i];
00892             const size_t listLen = pidLidList.size ();
00893             if (listLen > 0) {
00894               if (listLen > 1) {
00895                 locallyOneToOne_ = false;
00896               }
00897               // If there is some (PID,LID) pair for the current input
00898               // Map LID, then it makes sense to invoke the TieBreak
00899               // object to arbitrate between the options.  Even if
00900               // there is only one (PID,LID) pair, we still want to
00901               // give the TieBreak object a chance to do whatever it
00902               // likes to do, in terms of side effects (e.g., track
00903               // (PID,LID) pairs).
00904               const size_type index =
00905                 static_cast<size_type> (tie_break->selectedIndex (dirMapGid,
00906                                                                   pidLidList));
00907               PIDs_[i] = pidLidList[index].first;
00908               LIDs_[i] = pidLidList[index].second;
00909             }
00910             // else no GID specified by source map
00911           }
00912         }
00913       }
00914     }
00915 
00916     template<class LO, class GO, class NT>
00917     std::string
00918     DistributedNoncontiguousDirectory<LO, GO, NT>::description () const
00919     {
00920       std::ostringstream os;
00921       os << "DistributedNoncontiguousDirectory"
00922          << "<" << Teuchos::TypeNameTraits<LO>::name ()
00923          << ", " << Teuchos::TypeNameTraits<GO>::name ()
00924          << ", " << Teuchos::TypeNameTraits<NT>::name () << ">";
00925       return os.str ();
00926     }
00927 
00928     template<class LO, class GO, class NT>
00929     LookupStatus
00930     DistributedNoncontiguousDirectory<LO, GO, NT>::
00931     getEntriesImpl (const map_type& map,
00932                     const Teuchos::ArrayView<const GO> &globalIDs,
00933                     const Teuchos::ArrayView<int> &nodeIDs,
00934                     const Teuchos::ArrayView<LO> &localIDs,
00935                     const bool computeLIDs) const
00936     {
00937       using Teuchos::Array;
00938       using Teuchos::ArrayRCP;
00939       using Teuchos::ArrayView;
00940       using Teuchos::as;
00941       using Teuchos::RCP;
00942       using std::cerr;
00943       using std::endl;
00944       typedef typename Array<GO>::size_type size_type;
00945 
00946       RCP<const Teuchos::Comm<int> > comm = map.getComm ();
00947       const size_t numEntries = globalIDs.size ();
00948       const LO LINVALID = Teuchos::OrdinalTraits<LO>::invalid();
00949       LookupStatus res = AllIDsPresent;
00950 
00951       //
00952       // Set up directory structure.
00953       //
00954 
00955       // If we're computing LIDs, we also have to include them in each
00956       // packet, along with the GID and process ID.
00957       const int packetSize = computeLIDs ? 3 : 2;
00958 
00959       // For data distribution, we use: Surprise!  A Distributor!
00960       Distributor distor (comm);
00961 
00962       // Get directory locations for the requested list of entries.
00963       Array<int> dirImages (numEntries);
00964       res = directoryMap_->getRemoteIndexList (globalIDs, dirImages ());
00965       // Check for unfound globalIDs and set corresponding nodeIDs to -1
00966       size_t numMissing = 0;
00967       if (res == IDNotPresent) {
00968         for (size_t i=0; i < numEntries; ++i) {
00969           if (dirImages[i] == -1) {
00970             nodeIDs[i] = -1;
00971             if (computeLIDs) {
00972               localIDs[i] = LINVALID;
00973             }
00974             numMissing++;
00975           }
00976         }
00977       }
00978 
00979       Array<GO> sendGIDs;
00980       Array<int> sendImages;
00981       distor.createFromRecvs (globalIDs, dirImages (), sendGIDs, sendImages);
00982       const size_type numSends = sendGIDs.size ();
00983 
00984       //
00985       // mfh 13 Nov 2012:
00986       //
00987       // The code below temporarily stores LO, GO, and int values in
00988       // an array of global_size_t.  If one of the signed types (LO
00989       // and GO should both be signed) happened to be -1 (or some
00990       // negative number, but -1 is the one that came up today), then
00991       // conversion to global_size_t will result in a huge
00992       // global_size_t value, and thus conversion back may overflow.
00993       // (Teuchos::as doesn't know that we meant it to be an LO or GO
00994       // all along.)
00995       //
00996       // The overflow normally would not be a problem, since it would
00997       // just go back to -1 again.  However, Teuchos::as does range
00998       // checking on conversions in a debug build, so it throws an
00999       // exception (std::range_error) in this case.  Range checking is
01000       // generally useful in debug mode, so we don't want to disable
01001       // this behavior globally.
01002       //
01003       // We solve this problem by forgoing use of Teuchos::as for the
01004       // conversions below from LO, GO, or int to global_size_t, and
01005       // the later conversions back from global_size_t to LO, GO, or
01006       // int.
01007       //
01008       // I've recorded this discussion as Bug 5760.
01009       //
01010 
01011       //    global_size_t >= GO
01012       //    global_size_t >= size_t >= int
01013       //    global_size_t >= size_t >= LO
01014       // Therefore, we can safely store all of these in a global_size_t
01015       Array<global_size_t> exports (packetSize * numSends);
01016       {
01017         // Packet format:
01018         // - If computing LIDs: (GID, PID, LID)
01019         // - Otherwise:         (GID, PID)
01020         //
01021         // "PID" means "process ID" (a.k.a. "node ID," a.k.a. "rank").
01022 
01023         // Current position to which to write in exports array.  If
01024         // sending pairs, we pack the (GID, PID) pair for gid =
01025         // sendGIDs[k] in exports[2*k], exports[2*k+1].  If sending
01026         // triples, we pack the (GID, PID, LID) pair for gid =
01027         // sendGIDs[k] in exports[3*k, 3*k+1, 3*k+2].
01028         size_type exportsIndex = 0;
01029 
01030         if (useHashTables_) {
01031           for (size_type gidIndex = 0; gidIndex < numSends; ++gidIndex) {
01032             const GO curGID = sendGIDs[gidIndex];
01033             // Don't use as() here (see above note).
01034             exports[exportsIndex++] = static_cast<global_size_t> (curGID);
01035             const LO curLID = directoryMap_->getLocalElement (curGID);
01036             TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
01037               Teuchos::typeName (*this) << "::getEntriesImpl(): The Directory "
01038               "Map's global index " << curGID << " does not have a corresponding "
01039               "local index.  Please report this bug to the Tpetra developers.");
01040             // Don't use as() here (see above note).
01041             exports[exportsIndex++] = static_cast<global_size_t> (lidToPidTable_->get (curLID));
01042             if (computeLIDs) {
01043               // Don't use as() here (see above note).
01044               exports[exportsIndex++] = static_cast<global_size_t> (lidToLidTable_->get (curLID));
01045             }
01046           }
01047         } else {
01048           for (size_type gidIndex = 0; gidIndex < numSends; ++gidIndex) {
01049             const GO curGID = sendGIDs[gidIndex];
01050             // Don't use as() here (see above note).
01051             exports[exportsIndex++] = static_cast<global_size_t> (curGID);
01052             const LO curLID = directoryMap_->getLocalElement (curGID);
01053             TEUCHOS_TEST_FOR_EXCEPTION(curLID == LINVALID, std::logic_error,
01054               Teuchos::typeName (*this) << "::getEntriesImpl(): The Directory "
01055               "Map's global index " << curGID << " does not have a corresponding "
01056               "local index.  Please report this bug to the Tpetra developers.");
01057             // Don't use as() here (see above note).
01058             exports[exportsIndex++] = static_cast<global_size_t> (PIDs_[curLID]);
01059             if (computeLIDs) {
01060               // Don't use as() here (see above note).
01061               exports[exportsIndex++] = static_cast<global_size_t> (LIDs_[curLID]);
01062             }
01063           }
01064         }
01065 
01066         TEUCHOS_TEST_FOR_EXCEPTION(
01067           exportsIndex > exports.size (), std::logic_error,
01068           Teuchos::typeName (*this) << "::getEntriesImpl(): On Process " <<
01069           comm->getRank () << ", exportsIndex = " << exportsIndex <<
01070           " > exports.size() = " << exports.size () <<
01071           ".  Please report this bug to the Tpetra developers.");
01072       }
01073 
01074       TEUCHOS_TEST_FOR_EXCEPTION(
01075         numEntries < numMissing, std::logic_error,
01076         Teuchos::typeName (*this) << "::getEntriesImpl(): On Process "
01077         << comm->getRank () << ", numEntries = " << numEntries
01078         << " < numMissing = " << numMissing
01079         << ".  Please report this bug to the Tpetra developers.");
01080 
01081       //
01082       // mfh 13 Nov 2012: See note above on conversions between
01083       // global_size_t and LO, GO, or int.
01084       //
01085       const size_t numRecv = numEntries - numMissing;
01086 
01087       {
01088         const size_t importLen = packetSize * distor.getTotalReceiveLength ();
01089         const size_t requiredImportLen = numRecv * packetSize;
01090         const int myRank = comm->getRank ();
01091         TEUCHOS_TEST_FOR_EXCEPTION(
01092           importLen < requiredImportLen, std::logic_error,
01093           "Tpetra::Details::DistributedNoncontiguousDirectory::getEntriesImpl: "
01094           "On Process " << myRank << ": The 'imports' array must have length "
01095           "at least " << requiredImportLen << ", but its actual length is " <<
01096           importLen << ".  numRecv: " << numRecv << ", packetSize: " <<
01097           packetSize << ", numEntries (# GIDs): " << numEntries <<
01098           ", numMissing: " << numMissing << ": distor.getTotalReceiveLength(): "
01099           << distor.getTotalReceiveLength () << ".  " << std::endl <<
01100           "Distributor description: " << distor.description () << ".  "
01101           << std::endl <<
01102           "Please report this bug to the Tpetra developers.");
01103       }
01104 
01105       Array<global_size_t> imports (packetSize * distor.getTotalReceiveLength ());
01106       // FIXME (mfh 20 Mar 2014) One could overlap the sort2() below
01107       // with communication, by splitting this call into doPosts and
01108       // doWaits.  The code is still correct in this form, however.
01109       distor.doPostsAndWaits (exports ().getConst (), packetSize, imports ());
01110 
01111       Array<GO> sortedIDs (globalIDs); // deep copy (for later sorting)
01112       Array<GO> offset (numEntries); // permutation array (sort2 output)
01113       for (GO ii = 0; ii < static_cast<GO> (numEntries); ++ii) {
01114         offset[ii] = ii;
01115       }
01116       sort2 (sortedIDs.begin(), sortedIDs.begin() + numEntries, offset.begin());
01117 
01118       size_t importsIndex = 0;
01119       //typename Array<global_size_t>::iterator ptr = imports.begin();
01120       typedef typename Array<GO>::iterator IT;
01121 
01122       // we know these conversions are in range, because we loaded this data
01123       for (size_t i = 0; i < numRecv; ++i) {
01124         // Don't use as() here (see above note).
01125         const GO curGID = static_cast<GO> (imports[importsIndex++]);
01126         std::pair<IT, IT> p1 = std::equal_range (sortedIDs.begin(), sortedIDs.end(), curGID);
01127         if (p1.first != p1.second) {
01128           const size_t j = p1.first - sortedIDs.begin();
01129           // Don't use as() here (see above note).
01130           nodeIDs[offset[j]] = static_cast<int> (imports[importsIndex++]);
01131           if (computeLIDs) {
01132             // Don't use as() here (see above note).
01133             localIDs[offset[j]] = static_cast<LO> (imports[importsIndex++]);
01134           }
01135           if (nodeIDs[offset[j]] == -1) {
01136             res = IDNotPresent;
01137           }
01138         }
01139       }
01140 
01141       TEUCHOS_TEST_FOR_EXCEPTION(
01142         static_cast<size_t> (importsIndex) > static_cast<size_t> (imports.size ()),
01143         std::logic_error,
01144         "Tpetra::Details::DistributedNoncontiguousDirectory::getEntriesImpl: "
01145         "On Process " << comm->getRank () << ": importsIndex = " <<
01146         importsIndex << " > imports.size() = " << imports.size () << ".  "
01147         "numRecv: " << numRecv << ", packetSize: " << packetSize << ", "
01148         "numEntries (# GIDs): " << numEntries << ", numMissing: " << numMissing
01149         << ": distor.getTotalReceiveLength(): "
01150         << distor.getTotalReceiveLength () << ".  Please report this bug to "
01151         "the Tpetra developers.");
01152 
01153       return res;
01154     }
01155 
01156 
01157     template<class LO, class GO, class NT>
01158     bool
01159     DistributedNoncontiguousDirectory<LO, GO, NT>::
01160     isOneToOne (const Teuchos::Comm<int>& comm) const
01161     {
01162       if (oneToOneResult_ == ONE_TO_ONE_NOT_CALLED_YET) {
01163         const int lcl121 = isLocallyOneToOne () ? 1 : 0;
01164         int gbl121 = 0;
01165         Teuchos::reduceAll<int, int> (comm, Teuchos::REDUCE_MIN, lcl121,
01166                                       Teuchos::outArg (gbl121));
01167         oneToOneResult_ = (gbl121 == 1) ? ONE_TO_ONE_TRUE : ONE_TO_ONE_FALSE;
01168       }
01169       return (oneToOneResult_ == ONE_TO_ONE_TRUE);
01170     }
01171   } // namespace Details
01172 } // namespace Tpetra
01173 
01174 //
01175 // Explicit instantiation macro
01176 //
01177 // Must be expanded from within the Tpetra::Details namespace!
01178 //
01179 #define TPETRA_DIRECTORY_IMPL_INSTANT(LO,GO,NODE)                     \
01180   template class Directory< LO , GO , NODE >;                         \
01181   template class ReplicatedDirectory< LO , GO , NODE >;               \
01182   template class ContiguousUniformDirectory< LO, GO, NODE >;          \
01183   template class DistributedContiguousDirectory< LO , GO , NODE >;    \
01184   template class DistributedNoncontiguousDirectory< LO , GO , NODE >; \
01185 
01186 
01187 #endif // __Tpetra_DirectoryImpl_def_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines