Tpetra_Details_makeOptimizedColMap.hpp

Go to the documentation of this file.

/*
// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @HEADER
*/


#include "Tpetra_Map.hpp"
#include "Tpetra_Import.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_Details_Behavior.hpp"
#include "Teuchos_FancyOStream.hpp"
#include <memory>

namespace Tpetra {
namespace Details {

  template<class MapType>
  class OptColMap {
  public:
    using local_ordinal_type = typename MapType::local_ordinal_type;
    using global_ordinal_type = typename MapType::global_ordinal_type;
    using node_type = typename MapType::node_type;
    using map_type = ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>;
    using import_type = ::Tpetra::Import<local_ordinal_type, global_ordinal_type, node_type>;

    static Teuchos::RCP<const map_type>
    makeOptColMap (std::ostream& errStream,
                   bool& lclErr,
                   const map_type& domMap,
                   const map_type& colMap,
                   const import_type* /* oldImport */)
    {
      using ::Tpetra::Details::Behavior;
      using Teuchos::Array;
      using Teuchos::ArrayView;
      using Teuchos::FancyOStream;
      using Teuchos::getFancyOStream;
      using Teuchos::RCP;
      using Teuchos::rcp;
      using Teuchos::rcpFromRef;
      using std::endl;
      using LO = local_ordinal_type;
      using GO = global_ordinal_type;
      const char prefix[] = "Tpetra::Details::makeOptimizedColMap: ";

      RCP<const Teuchos::Comm<int> > comm = colMap.getComm ();
      std::ostream& err = errStream;

      const bool verbose = Behavior::verbose ("Tpetra::Details::makeOptimizedColMap");

      RCP<FancyOStream> outPtr = getFancyOStream (rcpFromRef (std::cerr));
      TEUCHOS_TEST_FOR_EXCEPTION
        (outPtr.is_null (), std::logic_error,
         "outPtr is null; this should never happen!");
      FancyOStream& out = *outPtr;
      Teuchos::OSTab tab1 (out);

      std::unique_ptr<std::string> verboseHeader;
      if (verbose) {
        std::ostringstream os;
        const int myRank = comm->getRank ();
        os << "Proc " << myRank << ": ";
        verboseHeader = std::unique_ptr<std::string> (new std::string (os.str ()));
      }
      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "Tpetra::Details::makeOptimizedColMap" << endl;
        out << os.str ();
      }

      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "Domain Map GIDs: [";
        const LO domMapLclNumInds = static_cast<LO> (domMap.getNodeNumElements ());
        for (LO lid = 0; lid < domMapLclNumInds; ++lid) {
          const GO gid = domMap.getGlobalElement (lid);
          os << gid;
          if (lid + LO (1) < domMapLclNumInds) {
            os << ", ";
          }
        }
        os << "]" << endl;
        out << os.str ();
      }

      const LO colMapLclNumInds = static_cast<LO> (colMap.getNodeNumElements ());

      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "Column Map GIDs: [";
        for (LO lid = 0; lid < colMapLclNumInds; ++lid) {
          const GO gid = colMap.getGlobalElement (lid);
          os << gid;
          if (lid + LO (1) < colMapLclNumInds) {
            os << ", ";
          }
        }
        os << "]" << endl;
        out << os.str ();
      }

      // Count remote GIDs.
      LO numOwnedGids = 0;
      LO numRemoteGids = 0;
      for (LO colMapLid = 0; colMapLid < colMapLclNumInds; ++colMapLid) {
        const GO colMapGid = colMap.getGlobalElement (colMapLid);
        if (domMap.isNodeGlobalElement (colMapGid)) {
          ++numOwnedGids;
        }
        else {
          ++numRemoteGids;
        }
      }

      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "- numOwnedGids: " << numOwnedGids << endl
           << *verboseHeader << "- numRemoteGids: " << numRemoteGids << endl;
        out << os.str ();
      }

      // Put all colMap GIDs on the calling process in a single array.
      // Owned GIDs go in front, and remote GIDs at the end.
      Array<GO> allGids (numOwnedGids + numRemoteGids);
      ArrayView<GO> ownedGids = allGids.view (0, numOwnedGids);
      ArrayView<GO> remoteGids = allGids.view (numOwnedGids, numRemoteGids);

      // Fill ownedGids and remoteGids (and therefore allGids).  We use
      // two loops, one to count (above) and one to fill (here), in
      // order to avoid dynamic memory allocation during the loop (in
      // this case, lots of calls to push_back()).  That will simplify
      // use of Kokkos to parallelize these loops later.
      LO ownedPos = 0;
      LO remotePos = 0;
      for (LO colMapLid = 0; colMapLid < colMapLclNumInds; ++colMapLid) {
        const GO colMapGid = colMap.getGlobalElement (colMapLid);
        if (domMap.isNodeGlobalElement (colMapGid)) {
          ownedGids[ownedPos++] = colMapGid;
        }
        else {
          remoteGids[remotePos++] = colMapGid;
        }
      }

      // If, for some reason, the running count doesn't match the
      // orignal count, fill in any remaining GID spots with an
      // obviously invalid value.  We don't want to stop yet, because
      // other processes might not have noticed this error; Map
      // construction is a collective, so we can't stop now.
      if (ownedPos != numOwnedGids) {
        lclErr = true;
        err << prefix << "On Process " << comm->getRank () << ", ownedPos = "
            << ownedPos << " != numOwnedGids = " << numOwnedGids << endl;
        for (LO colMapLid = ownedPos; colMapLid < numOwnedGids; ++colMapLid) {
          ownedGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
        }
      }
      if (remotePos != numRemoteGids) {
        lclErr = true;
        err << prefix << "On Process " << comm->getRank () << ", remotePos = "
            << remotePos << " != numRemoteGids = " << numRemoteGids << endl;
        for (LO colMapLid = remotePos; colMapLid < numRemoteGids; ++colMapLid) {
          remoteGids[colMapLid] = Teuchos::OrdinalTraits<GO>::invalid ();
        }
      }

      // Figure out what processes own what GIDs in the domain Map.
      // Initialize the output array of remote PIDs with the "invalid
      // process rank" -1, to help us test whether getRemoteIndexList
      // did its job.
      Array<int> remotePids (numRemoteGids, -1);
      const LookupStatus lookupStatus =
        domMap.getRemoteIndexList (remoteGids, remotePids ());

      // If any process returns IDNotPresent, then at least one of the
      // remote indices was not present in the domain Map.  This means
      // that the Import object cannot be constructed, because of
      // incongruity between the column Map and domain Map.  This
      // means that either the column Map or domain Map, or both, is
      // incorrect.
      const bool getRemoteIndexListFailed = (lookupStatus == IDNotPresent);
      if (getRemoteIndexListFailed) {
        lclErr = true;
        err << prefix << "On Process " << comm->getRank () << ", some indices "
          "in the input colMap (the original column Map) are not in domMap (the "
          "domain Map).  Either these indices or the domain Map is invalid.  "
          "Likely cause: For a nonsquare matrix, you must give the domain and "
          "range Maps as input to fillComplete." << endl;
      }

      // Check that getRemoteIndexList actually worked, by making sure
      // that none of the remote PIDs are -1.
      for (LO k = 0; k < numRemoteGids; ++k) {
        bool foundInvalidPid = false;
        if (remotePids[k] == -1) {
          foundInvalidPid = true;
          break;
        }
        if (foundInvalidPid) {
          lclErr = true;
          err << prefix << "On Process " << comm->getRank () << ", "
            "getRemoteIndexList returned -1 for the process ranks of "
            "one or more GIDs on this process." << endl;
        }
      }

      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "- Before sort2:" << endl
           << *verboseHeader << "-- ownedGids: " << Teuchos::toString (ownedGids) << endl
           << *verboseHeader << "-- remoteGids: " << Teuchos::toString (remoteGids) << endl
           << *verboseHeader << "-- allGids: " << Teuchos::toString (allGids ()) << endl;
        out << os.str ();
      }
      using Tpetra::sort2;
      sort2 (remotePids.begin (), remotePids.end (), remoteGids.begin ());
      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "- After sort2:" << endl
           << *verboseHeader << "-- ownedGids: " << Teuchos::toString (ownedGids) << endl
           << *verboseHeader << "-- remoteGids: " << Teuchos::toString (remoteGids) << endl
           << *verboseHeader << "-- allGids: " << Teuchos::toString (allGids ()) << endl;
        out << os.str ();
      }

      auto optColMap = rcp (new map_type (colMap.getGlobalNumElements (),
                                          allGids (),
                                          colMap.getIndexBase (),
                                          comm));
      if (verbose) {
        std::ostringstream os;
        os << *verboseHeader << "Tpetra::Details::makeOptimizedColMap: Done" << endl;
        out << os.str ();
      }
      return optColMap;
    }

    static std::pair<Teuchos::RCP<const map_type>,
                     Teuchos::RCP<import_type> >
    makeOptColMapAndImport (std::ostream& errStream,
                            bool& lclErr,
                            const map_type& domMap,
                            const map_type& colMap,
                            const import_type* oldImport)
    {
      using Teuchos::RCP;
      using Teuchos::rcp;

      // mfh 15 May 2018: For now, just call makeOptColMap, and use
      // the conventional two-Map (source and target) Import
      // constructor.
      RCP<const map_type> newColMap =
        makeOptColMap (errStream, lclErr, domMap, colMap, oldImport);
      RCP<import_type> imp (new import_type (rcp (new map_type (domMap)), newColMap));

      // FIXME (mfh 06 Jul 2014) This constructor throws a runtime
      // error, so I'm not using it for now.
      //
      // imp = rcp (new import_type (domMap, newColMap, remoteGids,
      //                             remotePids (), remoteLids (),
      //                             Teuchos::null, Teuchos::null));
      return std::make_pair (newColMap, imp);
    }
  };

  template<class MapType>
  Teuchos::RCP<const MapType>
  makeOptimizedColMap (std::ostream& errStream,
                       bool& lclErr,
                       const MapType& domMap,
                       const MapType& colMap,
                       const Tpetra::Import<
                         typename MapType::local_ordinal_type,
                         typename MapType::global_ordinal_type,
                         typename MapType::node_type>* oldImport = nullptr)
  {
    using map_type = ::Tpetra::Map<
      typename MapType::local_ordinal_type,
      typename MapType::global_ordinal_type,
      typename MapType::node_type>;
    using impl_type = OptColMap<map_type>;
    auto mapPtr = impl_type::makeOptColMap (errStream, lclErr,
                                            domMap, colMap, oldImport);
    return mapPtr;
  }

  template<class MapType>
  std::pair<Teuchos::RCP<const MapType>,
            Teuchos::RCP<typename OptColMap<MapType>::import_type> >
  makeOptimizedColMapAndImport (std::ostream& errStream,
                                bool& lclErr,
                                const MapType& domMap,
                                const MapType& colMap,
                                const typename OptColMap<MapType>::import_type* oldImport = nullptr)
  {
    using local_ordinal_type = typename MapType::local_ordinal_type;
    using global_ordinal_type = typename MapType::global_ordinal_type;
    using node_type = typename MapType::node_type;
    using map_type = ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>;
    using impl_type = OptColMap<map_type>;

    auto mapAndImp = impl_type::makeOptColMapAndImport (errStream, lclErr, domMap, colMap, oldImport);
    return std::make_pair (mapAndImp.first, mapAndImp.second);
  }

} // namespace Details
} // namespace Tpetra