MueLu_UncoupledAggregationFactory_kokkos_def.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 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
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef MUELU_UNCOUPLEDAGGREGATIONFACTORY_KOKKOS_DEF_HPP_
#define MUELU_UNCOUPLEDAGGREGATIONFACTORY_KOKKOS_DEF_HPP_

#include <climits>

#include <Xpetra_Map.hpp>
#include <Xpetra_Vector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>

#include "MueLu_UncoupledAggregationFactory_kokkos_decl.hpp"

#include "MueLu_OnePtAggregationAlgorithm_kokkos.hpp"
#include "MueLu_PreserveDirichletAggregationAlgorithm_kokkos.hpp"

#include "MueLu_AggregationPhase1Algorithm_kokkos.hpp"
#include "MueLu_AggregationPhase2aAlgorithm_kokkos.hpp"
#include "MueLu_AggregationPhase2bAlgorithm_kokkos.hpp"
#include "MueLu_AggregationPhase3Algorithm_kokkos.hpp"

#include "MueLu_Level.hpp"
#include "MueLu_LWGraph_kokkos.hpp"
#include "MueLu_Aggregates.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"

#include "KokkosGraph_Distance2ColorHandle.hpp"
#include "KokkosGraph_Distance2Color.hpp"
#include "KokkosGraph_MIS2.hpp"

namespace MueLu {

template <class LocalOrdinal, class GlobalOrdinal, class Node>
UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::UncoupledAggregationFactory_kokkos()
  : bDefinitionPhase_(true) {}

template <class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());

  // Aggregation parameters (used in aggregation algorithms)
  // TODO introduce local member function for each aggregation algorithm such that each aggregation algorithm can define its own parameters

  typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
  SET_VALID_ENTRY("aggregation: max agg size");
  SET_VALID_ENTRY("aggregation: min agg size");
  SET_VALID_ENTRY("aggregation: max selected neighbors");
  SET_VALID_ENTRY("aggregation: ordering");
  validParamList->getEntry("aggregation: ordering").setValidator(rcp(new validatorType(Teuchos::tuple<std::string>("natural", "graph", "random"), "aggregation: ordering")));
  SET_VALID_ENTRY("aggregation: deterministic");
  SET_VALID_ENTRY("aggregation: coloring algorithm");
  SET_VALID_ENTRY("aggregation: enable phase 1");
  SET_VALID_ENTRY("aggregation: enable phase 2a");
  SET_VALID_ENTRY("aggregation: enable phase 2b");
  SET_VALID_ENTRY("aggregation: enable phase 3");
  SET_VALID_ENTRY("aggregation: match ML phase1");
  SET_VALID_ENTRY("aggregation: match ML phase2a");
  SET_VALID_ENTRY("aggregation: match ML phase2b");
  SET_VALID_ENTRY("aggregation: phase3 avoid singletons");
  SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
  SET_VALID_ENTRY("aggregation: preserve Dirichlet points");
  SET_VALID_ENTRY("aggregation: allow user-specified singletons");
  SET_VALID_ENTRY("aggregation: phase 1 algorithm");
#undef SET_VALID_ENTRY

  // general variables needed in AggregationFactory
  validParamList->set<RCP<const FactoryBase> >("Graph", null, "Generating factory of the graph");
  validParamList->set<RCP<const FactoryBase> >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'");

  // special variables necessary for OnePtAggregationAlgorithm
  validParamList->set<std::string>("OnePt aggregate map name", "", "Name of input map for single node aggregates. (default='')");
  validParamList->set<std::string>("OnePt aggregate map factory", "", "Generating factory of (DOF) map for single node aggregates.");
  // validParamList->set< RCP<const FactoryBase> >("OnePt aggregate map factory",    NoFactory::getRCP(), "Generating factory of (DOF) map for single node aggregates.");

  return validParamList;
}

template <class LocalOrdinal, class GlobalOrdinal, class Node>
void UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
  Input(currentLevel, "Graph");
  Input(currentLevel, "DofsPerNode");

  const ParameterList& pL = GetParameterList();

  // request special data necessary for OnePtAggregationAlgorithm
  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
  if (mapOnePtName.length() > 0) {
    std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
    if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
      currentLevel.DeclareInput(mapOnePtName, NoFactory::get());
    } else {
      RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
      currentLevel.DeclareInput(mapOnePtName, mapOnePtFact.get());
    }
  }
}

template <class LocalOrdinal, class GlobalOrdinal, class Node>
void UncoupledAggregationFactory_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
    Build(Level& currentLevel) const {
  using execution_space    = typename LWGraph_kokkos::execution_space;
  using memory_space       = typename LWGraph_kokkos::memory_space;
  using local_ordinal_type = typename LWGraph_kokkos::local_ordinal_type;
  FactoryMonitor m(*this, "Build", currentLevel);

  ParameterList pL  = GetParameterList();
  bDefinitionPhase_ = false;  // definition phase is finished, now all aggregation algorithm information is fixed

  if (pL.get<int>("aggregation: max agg size") == -1)
    pL.set("aggregation: max agg size", INT_MAX);

  // define aggregation algorithms
  RCP<const FactoryBase> graphFact = GetFactory("Graph");

  // TODO Can we keep different aggregation algorithms over more Build calls?
  algos_.clear();
  algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm_kokkos(graphFact)));
  if (pL.get<bool>("aggregation: allow user-specified singletons") == true) algos_.push_back(rcp(new OnePtAggregationAlgorithm_kokkos(graphFact)));
  if (pL.get<bool>("aggregation: enable phase 1") == true) algos_.push_back(rcp(new AggregationPhase1Algorithm_kokkos(graphFact)));
  if (pL.get<bool>("aggregation: enable phase 2a") == true) algos_.push_back(rcp(new AggregationPhase2aAlgorithm_kokkos(graphFact)));
  if (pL.get<bool>("aggregation: enable phase 2b") == true) algos_.push_back(rcp(new AggregationPhase2bAlgorithm_kokkos(graphFact)));
  if (pL.get<bool>("aggregation: enable phase 3") == true) algos_.push_back(rcp(new AggregationPhase3Algorithm_kokkos(graphFact)));

  // Sanity Checking: match ML behavior is not supported in UncoupledAggregation_Kokkos in Phase 1 or Phase 2b, but is in 2a
  TEUCHOS_TEST_FOR_EXCEPTION(pL.get<bool>("aggregation: match ML phase1"), std::invalid_argument, "Option: 'aggregation: match ML phase1' is not supported in the Kokkos version of uncoupled aggregation");
  TEUCHOS_TEST_FOR_EXCEPTION(pL.get<bool>("aggregation: match ML phase2b"), std::invalid_argument, "Option: 'aggregation: match ML phase2b' is not supported in the Kokkos version of uncoupled aggregation");

  std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name");
  RCP<Map> OnePtMap        = Teuchos::null;
  if (mapOnePtName.length()) {
    std::string mapOnePtFactName = pL.get<std::string>("OnePt aggregate map factory");
    if (mapOnePtFactName == "" || mapOnePtFactName == "NoFactory") {
      OnePtMap = currentLevel.Get<RCP<Map> >(mapOnePtName, NoFactory::get());
    } else {
      RCP<const FactoryBase> mapOnePtFact = GetFactory(mapOnePtFactName);
      OnePtMap                            = currentLevel.Get<RCP<Map> >(mapOnePtName, mapOnePtFact.get());
    }
  }

  RCP<const LWGraph_kokkos> graph = Get<RCP<LWGraph_kokkos> >(currentLevel, "Graph");

  // Build
  RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
  aggregates->setObjectLabel("UC");

  const LO numRows = graph->GetNodeNumVertices();

  // construct aggStat information
  Kokkos::View<unsigned*, typename LWGraph_kokkos::device_type> aggStat(Kokkos::ViewAllocateWithoutInitializing("aggregation status"),
                                                                        numRows);
  Kokkos::deep_copy(aggStat, READY);

  // LBV on Sept 06 2019: re-commenting out the dirichlet boundary map
  // even if the map is correctly extracted from the graph, aggStat is
  // now a Kokkos::View<unsinged*, memory_space> and filling it will
  // require a parallel_for or to copy it to the Host which is not really
  // good from a performance point of view.
  // If dirichletBoundaryMap was an actual Xpetra::Map, one could call
  // getLocalMap to have a Kokkos::View on the appropriate memory_space
  // instead of an ArrayRCP.
  {
    auto dirichletBoundaryMap = graph->GetBoundaryNodeMap();
    Kokkos::parallel_for(
        "MueLu - UncoupledAggregation: tagging boundary nodes in aggStat",
        Kokkos::RangePolicy<local_ordinal_type, execution_space>(0, numRows),
        KOKKOS_LAMBDA(const local_ordinal_type nodeIdx) {
          if (dirichletBoundaryMap(nodeIdx) == true) {
            aggStat(nodeIdx) = BOUNDARY;
          }
        });
  }

  LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
  GO indexBase    = graph->GetDomainMap()->getIndexBase();

  /* FIXME: This chunk of code is still executing on the host */
  if (OnePtMap != Teuchos::null) {
    typename Kokkos::View<unsigned*, typename LWGraph_kokkos::device_type>::HostMirror aggStatHost = Kokkos::create_mirror_view(aggStat);
    Kokkos::deep_copy(aggStatHost, aggStat);

    for (LO i = 0; i < numRows; i++) {
      // reconstruct global row id (FIXME only works for contiguous maps)
      GO grid = (graph->GetDomainMap()->getGlobalElement(i) - indexBase) * nDofsPerNode + indexBase;

      for (LO kr = 0; kr < nDofsPerNode; kr++)
        if (OnePtMap->isNodeGlobalElement(grid + kr))
          aggStatHost(i) = ONEPT;
    }

    Kokkos::deep_copy(aggStat, aggStatHost);
  }

  const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
  GO numGlobalRows                          = 0;
  if (IsPrint(Statistics1))
    MueLu_sumAll(comm, as<GO>(numRows), numGlobalRows);

  LO numNonAggregatedNodes = numRows;
  std::string aggAlgo      = pL.get<std::string>("aggregation: coloring algorithm");
  if (aggAlgo == "mis2 coarsening" || aggAlgo == "mis2 aggregation") {
    SubFactoryMonitor sfm(*this, "Algo \"MIS2\"", currentLevel);
    using graph_t      = typename LWGraph_kokkos::local_graph_type;
    using device_t     = typename graph_t::device_type;
    using exec_space   = typename device_t::execution_space;
    using rowmap_t     = typename graph_t::row_map_type;
    using colinds_t    = typename graph_t::entries_type;
    using lno_t        = typename colinds_t::non_const_value_type;
    rowmap_t aRowptrs  = graph->getRowPtrs();
    colinds_t aColinds = graph->getEntries();
    lno_t numAggs      = 0;
    typename colinds_t::non_const_type labels;

    if (aggAlgo == "mis2 coarsening") {
      if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: MIS-2 coarsening" << std::endl;
      labels = KokkosGraph::graph_mis2_coarsen<device_t, rowmap_t, colinds_t>(aRowptrs, aColinds, numAggs);
    } else if (aggAlgo == "mis2 aggregation") {
      if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: MIS-2 aggregation" << std::endl;
      labels = KokkosGraph::graph_mis2_aggregate<device_t, rowmap_t, colinds_t>(aRowptrs, aColinds, numAggs);
    }
    auto vertex2AggId = aggregates->GetVertex2AggId()->getDeviceLocalView(Xpetra::Access::ReadWrite);
    auto procWinner   = aggregates->GetProcWinner()->getDeviceLocalView(Xpetra::Access::OverwriteAll);
    int rank          = comm->getRank();
    Kokkos::parallel_for(
        Kokkos::RangePolicy<exec_space>(0, numRows),
        KOKKOS_LAMBDA(lno_t i) {
          procWinner(i, 0) = rank;
          if (aggStat(i) == READY) {
            aggStat(i)         = AGGREGATED;
            vertex2AggId(i, 0) = labels(i);
          }
        });
    numNonAggregatedNodes = 0;
    aggregates->SetNumAggregates(numAggs);
  } else {
    {
      SubFactoryMonitor sfm(*this, "Algo \"Graph Coloring\"", currentLevel);

      // LBV on Sept 06 2019: the note below is a little worrisome,
      // can we guarantee that MueLu is never used on a non-symmetric
      // graph?
      // note: just using colinds_view in place of scalar_view_t type
      // (it won't be used at all by symbolic SPGEMM)
      using graph_t      = typename LWGraph_kokkos::local_graph_type;
      using KernelHandle = KokkosKernels::Experimental::
          KokkosKernelsHandle<typename graph_t::row_map_type::value_type,
                              typename graph_t::entries_type::value_type,
                              typename graph_t::entries_type::value_type,
                              typename graph_t::device_type::execution_space,
                              typename graph_t::device_type::memory_space,
                              typename graph_t::device_type::memory_space>;
      KernelHandle kh;
      // leave gc algorithm choice as the default
      kh.create_distance2_graph_coloring_handle();

      // get the distance-2 graph coloring handle
      auto coloringHandle = kh.get_distance2_graph_coloring_handle();

      // Set the distance-2 graph coloring algorithm to use.
      // Options:
      //     COLORING_D2_DEFAULT        - Let the kernel handle pick the variation
      //     COLORING_D2_SERIAL         - Use the legacy serial-only implementation
      //     COLORING_D2_VB             - Use the parallel vertex based direct method
      //     COLORING_D2_VB_BIT         - Same as VB but using the bitvector forbidden array
      //     COLORING_D2_VB_BIT_EF      - Add experimental edge-filtering to VB_BIT
      //     COLORING_D2_NB_BIT         - Net-based coloring (generally the fastest)
      if (pL.get<bool>("aggregation: deterministic") == true) {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_SERIAL);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: serial" << std::endl;
      } else if (aggAlgo == "serial") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_SERIAL);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: serial" << std::endl;
      } else if (aggAlgo == "default") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_DEFAULT);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: default" << std::endl;
      } else if (aggAlgo == "vertex based") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_VB);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: vertex based" << std::endl;
      } else if (aggAlgo == "vertex based bit set") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_VB_BIT);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: vertex based bit set" << std::endl;
      } else if (aggAlgo == "edge filtering") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_VB_BIT_EF);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: edge filtering" << std::endl;
      } else if (aggAlgo == "net based bit set") {
        coloringHandle->set_algorithm(KokkosGraph::COLORING_D2_NB_BIT);
        if (IsPrint(Statistics1)) GetOStream(Statistics1) << "  algorithm: net based bit set" << std::endl;
      } else {
        TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unrecognized distance 2 coloring algorithm, valid options are: serial, default, matrix squared, vertex based, vertex based bit set, edge filtering")
      }

      // Create device views for graph rowptrs/colinds
      typename graph_t::row_map_type aRowptrs = graph->getRowPtrs();
      typename graph_t::entries_type aColinds = graph->getEntries();

      // run d2 graph coloring
      // graph is symmetric so row map/entries and col map/entries are the same
      {
        SubFactoryMonitor sfm2(*this, "Algo \"Graph Coloring\": KokkosGraph Call", currentLevel);  // CMS HACK
        KokkosGraph::Experimental::graph_color_distance2(&kh, numRows, aRowptrs, aColinds);
      }

      // extract the colors and store them in the aggregates
      aggregates->SetGraphColors(coloringHandle->get_vertex_colors());
      aggregates->SetGraphNumColors(static_cast<LO>(coloringHandle->get_num_colors()));

      // clean up coloring handle
      kh.destroy_distance2_graph_coloring_handle();
    }

    if (IsPrint(Statistics1)) {
      GetOStream(Statistics1) << "  num colors: " << aggregates->GetGraphNumColors() << std::endl;
    }
    GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
    for (size_t a = 0; a < algos_.size(); a++) {
      std::string phase = algos_[a]->description();
      SubFactoryMonitor sfm2(*this, "Algo \"" + phase + "\"", currentLevel);

      int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
      algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
      algos_[a]->SetProcRankVerbose(oldRank);

      if (IsPrint(Statistics1)) {
        GO numLocalAggregated = numRows - numNonAggregatedNodes, numGlobalAggregated = 0;
        GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
        MueLu_sumAll(comm, numLocalAggregated, numGlobalAggregated);
        MueLu_sumAll(comm, numLocalAggs, numGlobalAggs);

        double aggPercent = 100 * as<double>(numGlobalAggregated) / as<double>(numGlobalRows);
        if (aggPercent > 99.99 && aggPercent < 100.00) {
          // Due to round off (for instance, for 140465733/140466897), we could
          // get 100.00% display even if there are some remaining nodes. This
          // is bad from the users point of view. It is much better to change
          // it to display 99.99%.
          aggPercent = 99.99;
        }
        GetOStream(Statistics1) << "  aggregated : " << (numGlobalAggregated - numGlobalAggregatedPrev) << " (phase), " << std::fixed
                                << std::setprecision(2) << numGlobalAggregated << "/" << numGlobalRows << " [" << aggPercent << "%] (total)\n"
                                << "  remaining  : " << numGlobalRows - numGlobalAggregated << "\n"
                                << "  aggregates : " << numGlobalAggs - numGlobalAggsPrev << " (phase), " << numGlobalAggs << " (total)" << std::endl;
        numGlobalAggregatedPrev = numGlobalAggregated;
        numGlobalAggsPrev       = numGlobalAggs;
      }
    }
  }

  TEUCHOS_TEST_FOR_EXCEPTION(numNonAggregatedNodes, Exceptions::RuntimeError, "MueLu::UncoupledAggregationFactory::Build: Leftover nodes found! Error!");

  aggregates->AggregatesCrossProcessors(false);
  aggregates->ComputeAggregateSizes(true /*forceRecompute*/);

  Set(currentLevel, "Aggregates", aggregates);
}

}  // namespace MueLu

#endif /* MUELU_UNCOUPLEDAGGREGATIONFACTORY_DEF_HPP_ */