MueLu_HybridAggregationFactory_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_HYBRIDAGGREGATIONFACTORY_DEF_HPP_
#define MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP_

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

#include "MueLu_HybridAggregationFactory_decl.hpp"

// Uncoupled Agg
#include "MueLu_InterfaceAggregationAlgorithm.hpp"
#include "MueLu_OnePtAggregationAlgorithm.hpp"
#include "MueLu_PreserveDirichletAggregationAlgorithm.hpp"
#include "MueLu_IsolatedNodeAggregationAlgorithm.hpp"

#include "MueLu_AggregationPhase1Algorithm.hpp"
#include "MueLu_AggregationPhase2aAlgorithm.hpp"
#include "MueLu_AggregationPhase2bAlgorithm.hpp"
#include "MueLu_AggregationPhase3Algorithm.hpp"

// Structured Agg
#include "MueLu_AggregationStructuredAlgorithm.hpp"
#include "MueLu_UncoupledIndexManager.hpp"
//#include "MueLu_LocalLexicographicIndexManager.hpp"
//#include "MueLu_GlobalLexicographicIndexManager.hpp"

// Shared
#include "MueLu_Level.hpp"
#include "MueLu_GraphBase.hpp"
#include "MueLu_Aggregates.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_AmalgamationInfo.hpp"


namespace MueLu {

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

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

    typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
    // From UncoupledAggregationFactory
    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: 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: preserve Dirichlet points");
    SET_VALID_ENTRY("aggregation: allow user-specified singletons");
    SET_VALID_ENTRY("aggregation: error on nodes with no on-rank neighbors");
    SET_VALID_ENTRY("aggregation: phase2a include root");
    SET_VALID_ENTRY("aggregation: phase3 avoid singletons");

    // From StructuredAggregationFactory
    SET_VALID_ENTRY("aggregation: coarsening rate");
    SET_VALID_ENTRY("aggregation: coarsening order");
    SET_VALID_ENTRY("aggregation: number of spatial dimensions");

    // From HybridAggregationFactory
    SET_VALID_ENTRY("aggregation: use interface aggregation");
#undef  SET_VALID_ENTRY

    /* From UncoupledAggregation */
    // 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.");

    // InterfaceAggregation parameters
    validParamList->set<std::string>             ("Interface aggregate map name",                "",
                                                  "Name of input map for interface aggregates. (default='')");
    validParamList->set<std::string>             ("Interface aggregate map factory",             "",
                                                  "Generating factory of (DOF) map for interface aggregates.");
    validParamList->set<RCP<const FactoryBase> > ("interfacesDimensions",             Teuchos::null,
                                                  "Describes the dimensions of all the interfaces on this rank.");
    validParamList->set<RCP<const FactoryBase> > ("nodeOnInterface",                  Teuchos::null,
                                                  "List the LIDs of the nodes on any interface.");

    /* From StructuredAggregation */
    // general variables needed in AggregationFactory
    validParamList->set<RCP<const FactoryBase> >("numDimensions",           Teuchos::null,
                                                 "Number of spatial dimension provided by CoordinatesTransferFactory.");
    validParamList->set<RCP<const FactoryBase> >("lNodesPerDim",                     Teuchos::null,
                                                 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");


    // Hybrid Aggregation Params
    validParamList->set<RCP<const FactoryBase> > ("aggregationRegionType",            Teuchos::null,
                                                  "Type of aggregation to use on the region (\"structured\" or \"uncoupled\")");

    return validParamList;
  }

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

    ParameterList pL = GetParameterList();



    /* StructuredAggregation */

    // Request the local number of nodes per dimensions
    if(currentLevel.GetLevelID() == 0) {
      if(currentLevel.IsAvailable("aggregationRegionType", NoFactory::get())) {
        currentLevel.DeclareInput("aggregationRegionType", NoFactory::get(), this);
      } else {
        TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("aggregationRegionType",NoFactory::get()),
                                   Exceptions::RuntimeError,
                                   "Aggregation region type was not provided by the user!");
      }
      if(currentLevel.IsAvailable("numDimensions", NoFactory::get())) {
        currentLevel.DeclareInput("numDimensions", NoFactory::get(), this);
      } else {
        TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("numDimensions", NoFactory::get()),
                                   Exceptions::RuntimeError,
                                   "numDimensions was not provided by the user on level0!");
      }
      if(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
        currentLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
      } else {
        TEUCHOS_TEST_FOR_EXCEPTION(currentLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
                                   Exceptions::RuntimeError,
                                   "lNodesPerDim was not provided by the user on level0!");
      }
    } else {
      Input(currentLevel, "aggregationRegionType");
      Input(currentLevel, "numDimensions");
      Input(currentLevel, "lNodesPerDim");
    }



    /* UncoupledAggregation */
    Input(currentLevel, "DofsPerNode");

    // request special data necessary for InterfaceAggregation
    if (pL.get<bool>("aggregation: use interface aggregation") == true){
      if(currentLevel.GetLevelID() == 0) {
        if(currentLevel.IsAvailable("interfacesDimensions", NoFactory::get())) {
          currentLevel.DeclareInput("interfacesDimensions", NoFactory::get(), this);
        } else {
          TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("interfacesDimensions", NoFactory::get()),
                                       Exceptions::RuntimeError,
                                       "interfacesDimensions was not provided by the user on level0!");
        }
        if(currentLevel.IsAvailable("nodeOnInterface", NoFactory::get())) {
          currentLevel.DeclareInput("nodeOnInterface", NoFactory::get(), this);
        } else {
          TEUCHOS_TEST_FOR_EXCEPTION(!currentLevel.IsAvailable("nodeOnInterface", NoFactory::get()),
                                       Exceptions::RuntimeError,
                                       "nodeOnInterface was not provided by the user on level0!");
        }
      } else {
        Input(currentLevel, "interfacesDimensions");
        Input(currentLevel, "nodeOnInterface");
      }
    }

    // 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());
      }
    }
  } // DeclareInput()

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void HybridAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::
  Build(Level &currentLevel) const {
    FactoryMonitor m(*this, "Build", currentLevel);

    RCP<Teuchos::FancyOStream> out;
    if(const char* dbg = std::getenv("MUELU_HYBRIDAGGREGATION_DEBUG")) {
      out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
      out->setShowAllFrontMatter(false).setShowProcRank(true);
    } else {
      out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
    }

    *out << "Entering hybrid aggregation" << std::endl;

    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");

    // General problem informations are gathered from data stored in the problem matix.
    RCP<const GraphBase> graph  = Get< RCP<GraphBase> >(currentLevel, "Graph");
    RCP<const Map> fineMap      = graph->GetDomainMap();
    const int myRank            = fineMap->getComm()->getRank();
    const int numRanks          = fineMap->getComm()->getSize();

    out->setProcRankAndSize(graph->GetImportMap()->getComm()->getRank(),
                            graph->GetImportMap()->getComm()->getSize());

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

    // construct aggStat information
    const LO numRows = graph->GetNodeNumVertices();
    std::vector<unsigned> aggStat(numRows, READY);

    // Get aggregation type for region
    std::string regionType;
    if(currentLevel.GetLevelID() == 0) {
      // On level 0, data is provided by applications and has no associated factory.
      regionType = currentLevel.Get<std::string>("aggregationRegionType", NoFactory::get());
    } else {
      // On level > 0, data is provided directly by generating factories.
      regionType = Get< std::string >(currentLevel, "aggregationRegionType");
    }

    int numDimensions = 0;
    if(currentLevel.GetLevelID() == 0) {
      // On level 0, data is provided by applications and has no associated factory.
      numDimensions = currentLevel.Get<int>("numDimensions", NoFactory::get());
    } else {
      // On level > 0, data is provided directly by generating factories.
      numDimensions = Get<int>(currentLevel, "numDimensions");
    }

    // Get the coarsening rate (potentially used for both structured and uncoupled aggregation if interface)
    std::string coarseningRate = pL.get<std::string>("aggregation: coarsening rate");
    Teuchos::Array<LO> coarseRate;
    try {
      coarseRate = Teuchos::fromStringToArray<LO>(coarseningRate);
    } catch(const Teuchos::InvalidArrayStringRepresentation& e) {
      GetOStream(Errors,-1) << " *** \"aggregation: coarsening rate\" must be a string convertible into an array! *** "
                            << std::endl;
      throw e;
    }
    TEUCHOS_TEST_FOR_EXCEPTION((coarseRate.size() > 1) && (coarseRate.size() < numDimensions),
                               Exceptions::RuntimeError,
                               "\"aggregation: coarsening rate\" must have at least as many"
                               " components as the number of spatial dimensions in the problem.");

    algos_.clear();
    LO numNonAggregatedNodes = numRows;
    if (regionType == "structured") {
      // Add AggregationStructuredAlgorithm
      algos_.push_back(rcp(new AggregationStructuredAlgorithm(graphFact)));

      // Since we want to operate on nodes and not dof, we need to modify the rowMap in order to
      // obtain a nodeMap.
      const int interpolationOrder = pL.get<int>("aggregation: coarsening order");
      Array<LO> lFineNodesPerDir(3);
      if(currentLevel.GetLevelID() == 0) {
        // On level 0, data is provided by applications and has no associated factory.
        lFineNodesPerDir = currentLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
      } else {
        // On level > 0, data is provided directly by generating factories.
        lFineNodesPerDir = Get<Array<LO> >(currentLevel, "lNodesPerDim");
      }

      // Set lFineNodesPerDir to 1 for directions beyond numDimensions
      for(int dim = numDimensions; dim < 3; ++dim) {
        lFineNodesPerDir[dim] = 1;
      }

      // Now that we have extracted info from the level, create the IndexManager
      RCP<MueLu::IndexManager<LO,GO,NO> > geoData;
      geoData = rcp(new MueLu::UncoupledIndexManager<LO,GO,NO>(fineMap->getComm(),
                                                               false,
                                                               numDimensions,
                                                               interpolationOrder,
                                                               myRank,
                                                               numRanks,
                                                               Array<GO>(3, -1),
                                                               lFineNodesPerDir,
                                                               coarseRate));

      TEUCHOS_TEST_FOR_EXCEPTION(fineMap->getNodeNumElements()
                                 != static_cast<size_t>(geoData->getNumLocalFineNodes()),
                                 Exceptions::RuntimeError,
                                 "The local number of elements in the graph's map is not equal to "
                                 "the number of nodes given by: lNodesPerDim!");

      aggregates->SetIndexManager(geoData);
      aggregates->SetNumAggregates(geoData->getNumLocalCoarseNodes());

      Set(currentLevel, "lCoarseNodesPerDim", geoData->getLocalCoarseNodesPerDir());

    } // end structured aggregation setup

    if (regionType == "uncoupled"){
      // Add unstructred aggregation phases
      algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm(graphFact)));
      if (pL.get<bool>("aggregation: use interface aggregation")       == true)   algos_.push_back(rcp(new InterfaceAggregationAlgorithm         (graphFact)));
      if (pL.get<bool>("aggregation: allow user-specified singletons") == true)   algos_.push_back(rcp(new OnePtAggregationAlgorithm             (graphFact)));
      if (pL.get<bool>("aggregation: enable phase 1" )                 == true)   algos_.push_back(rcp(new AggregationPhase1Algorithm            (graphFact)));
      if (pL.get<bool>("aggregation: enable phase 2a")                 == true)   algos_.push_back(rcp(new AggregationPhase2aAlgorithm           (graphFact)));
      if (pL.get<bool>("aggregation: enable phase 2b")                 == true)   algos_.push_back(rcp(new AggregationPhase2bAlgorithm           (graphFact)));
      if (pL.get<bool>("aggregation: enable phase 3" )                 == true)   algos_.push_back(rcp(new AggregationPhase3Algorithm            (graphFact)));

      *out << " Build interface aggregates" << std::endl;
      // interface
      if (pL.get<bool>("aggregation: use interface aggregation") == true) {
        BuildInterfaceAggregates(currentLevel, aggregates, aggStat, numNonAggregatedNodes,
                                 coarseRate);
      }

      *out << "Treat Dirichlet BC" << std::endl;
      // Dirichlet boundary
      ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
      if (dirichletBoundaryMap != Teuchos::null)
        for (LO i = 0; i < numRows; i++)
          if (dirichletBoundaryMap[i] == true)
            aggStat[i] = BOUNDARY;

      // OnePt 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());
        }
      }

      LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
      GO indexBase = graph->GetDomainMap()->getIndexBase();
      if (OnePtMap != Teuchos::null) {
        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))
              aggStat[i] = ONEPT;
        }
      }

      // Create a fake lCoarseNodesPerDir for CoordinatesTranferFactory
      Array<LO> lCoarseNodesPerDir(3,-1);
      Set(currentLevel, "lCoarseNodesPerDim", lCoarseNodesPerDir);
    } // end uncoupled aggregation setup

    aggregates->AggregatesCrossProcessors(false); // No coupled aggregation

    *out << "Run all the algorithms on the local rank" << std::endl;
    for (size_t a = 0; a < algos_.size(); a++) {
      std::string phase = algos_[a]->description();
      SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
      *out << regionType <<" | Executing phase " << a << std::endl;

      int oldRank = algos_[a]->SetProcRankVerbose(this->GetProcRankVerbose());
      algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
      algos_[a]->SetProcRankVerbose(oldRank);
      *out << regionType <<" | Done Executing phase " << a << std::endl;
    }

    *out << "Compute statistics on aggregates" << std::endl;
    aggregates->ComputeAggregateSizes(true/*forceRecompute*/);

    Set(currentLevel, "Aggregates",                  aggregates);
    Set(currentLevel, "numDimensions",               numDimensions);
    Set(currentLevel, "aggregationRegionTypeCoarse", regionType);

    GetOStream(Statistics1) << aggregates->description() << std::endl;
    *out << "HybridAggregation done!" << std::endl;
  }

  template <class LocalOrdinal, class GlobalOrdinal, class Node>
  void HybridAggregationFactory<LocalOrdinal, GlobalOrdinal, Node>::
  BuildInterfaceAggregates(Level& currentLevel, RCP<Aggregates> aggregates,
                           std::vector<unsigned>& aggStat, LO& numNonAggregatedNodes,
                           Array<LO> coarseRate) const {
    FactoryMonitor m(*this, "BuildInterfaceAggregates", currentLevel);

    RCP<Teuchos::FancyOStream> out;
    if(const char* dbg = std::getenv("MUELU_HYBRIDAGGREGATION_DEBUG")) {
      out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
      out->setShowAllFrontMatter(false).setShowProcRank(true);
    } else {
      out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
    }

    // Extract and format input data for algo
    if(coarseRate.size() == 1) {coarseRate.resize(3, coarseRate[0]);}
    ArrayRCP<LO> vertex2AggId      = aggregates->GetVertex2AggId()->getDataNonConst(0);
    ArrayRCP<LO> procWinner        = aggregates->GetProcWinner()  ->getDataNonConst(0);
    Array<LO> interfacesDimensions = Get<Array<LO> >(currentLevel, "interfacesDimensions");
    Array<LO> nodesOnInterfaces    = Get<Array<LO> >(currentLevel, "nodeOnInterface");
    const int numInterfaces        = interfacesDimensions.size() / 3;
    const int myRank               = aggregates->GetMap()->getComm()->getRank();

    // Create coarse level container to gather data on the fly
    Array<LO> coarseInterfacesDimensions(interfacesDimensions.size());
    Array<LO> nodesOnCoarseInterfaces;
    { // Scoping the temporary variables...
      LO endRate, totalNumCoarseNodes = 0, numCoarseNodes;
      for(int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
        numCoarseNodes = 1;
        for(int dim = 0; dim < 3; ++dim) {
          endRate = (interfacesDimensions[3*interfaceIdx + dim] - 1) % coarseRate[dim];
          if(interfacesDimensions[3*interfaceIdx + dim] == 1) {
            coarseInterfacesDimensions[3*interfaceIdx + dim] = 1;
          } else {
            coarseInterfacesDimensions[3*interfaceIdx + dim]
              = (interfacesDimensions[3*interfaceIdx+dim]-1) / coarseRate[dim] + 2;
            if(endRate==0){ coarseInterfacesDimensions[3*interfaceIdx + dim]--;}
          }
          numCoarseNodes *= coarseInterfacesDimensions[3*interfaceIdx + dim];
        }
        totalNumCoarseNodes += numCoarseNodes;
      }
      nodesOnCoarseInterfaces.resize(totalNumCoarseNodes, -1);
    }

    Array<LO> endRate(3);
    LO interfaceOffset = 0, aggregateCount = 0, coarseNodeCount = 0;
    for(int interfaceIdx = 0; interfaceIdx < numInterfaces; ++interfaceIdx) {
      ArrayView<LO> fineNodesPerDim   = interfacesDimensions(3*interfaceIdx, 3);
      ArrayView<LO> coarseNodesPerDim = coarseInterfacesDimensions(3*interfaceIdx, 3);
      LO numInterfaceNodes = 1, numCoarseNodes = 1;
      for(int dim = 0; dim < 3; ++dim) {
        numInterfaceNodes *= fineNodesPerDim[dim];
        numCoarseNodes    *= coarseNodesPerDim[dim];
        endRate[dim]       = fineNodesPerDim[dim] % coarseRate[dim];
      }
      ArrayView<LO> interfaceNodes = nodesOnInterfaces(interfaceOffset, numInterfaceNodes);

      interfaceOffset += numInterfaceNodes;

      LO rem, rate, fineNodeIdx;
      Array<LO> nodeIJK(3), coarseIJK(3), rootIJK(3);
      // First find treat coarse nodes as they generate the aggregate IDs
      // and they might be repeated on multiple interfaces (think corners and edges).
      for(LO coarseNodeIdx = 0; coarseNodeIdx < numCoarseNodes; ++coarseNodeIdx) {
        coarseIJK[2] = coarseNodeIdx / (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
        rem          = coarseNodeIdx % (coarseNodesPerDim[0]*coarseNodesPerDim[1]);
        coarseIJK[1] = rem / coarseNodesPerDim[0];
        coarseIJK[0] = rem % coarseNodesPerDim[0];

        for(LO dim = 0; dim < 3; ++dim) {
          if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
            nodeIJK[dim] = fineNodesPerDim[dim] - 1;
          } else {
            nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
          }
        }
        fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];

        if(aggStat[interfaceNodes[fineNodeIdx]] == READY) {
          vertex2AggId[interfaceNodes[fineNodeIdx]] = aggregateCount;
          procWinner[interfaceNodes[fineNodeIdx]]   = myRank;
          aggStat[interfaceNodes[fineNodeIdx]]      = AGGREGATED;
          ++aggregateCount;
          --numNonAggregatedNodes;
        }
        nodesOnCoarseInterfaces[coarseNodeCount] = vertex2AggId[interfaceNodes[fineNodeIdx]];
        ++coarseNodeCount;
      }

      // Now loop over all the node on the interface
      // skip the coarse nodes as they are already aggregated
      // and find the appropriate aggregate ID for the fine nodes.
      for(LO nodeIdx = 0; nodeIdx < numInterfaceNodes; ++nodeIdx) {

        // If the node is already aggregated skip it!
        if(aggStat[interfaceNodes[nodeIdx]] == AGGREGATED) {continue;}

        nodeIJK[2] = nodeIdx / (fineNodesPerDim[0]*fineNodesPerDim[1]);
        rem        = nodeIdx % (fineNodesPerDim[0]*fineNodesPerDim[1]);
        nodeIJK[1] = rem / fineNodesPerDim[0];
        nodeIJK[0] = rem % fineNodesPerDim[0];

        for(int dim = 0; dim < 3; ++dim) {
          coarseIJK[dim] = nodeIJK[dim] / coarseRate[dim];
          rem            = nodeIJK[dim] % coarseRate[dim];
          if(nodeIJK[dim] < fineNodesPerDim[dim] - endRate[dim]) {
            rate = coarseRate[dim];
          } else {
            rate = endRate[dim];
          }
          if(rem > (rate / 2)) {++coarseIJK[dim];}
          if(coarseNodesPerDim[dim] - coarseIJK[dim] > fineNodesPerDim[dim]-nodeIJK[dim]){
            ++coarseIJK[dim];
          }
        }

        for(LO dim = 0; dim < 3; ++dim) {
          if(coarseIJK[dim] == coarseNodesPerDim[dim] - 1) {
            nodeIJK[dim] = fineNodesPerDim[dim] - 1;
          } else {
            nodeIJK[dim] = coarseIJK[dim]*coarseRate[dim];
          }
        }
        fineNodeIdx = (nodeIJK[2]*fineNodesPerDim[1] + nodeIJK[1])*fineNodesPerDim[0] + nodeIJK[0];

        vertex2AggId[interfaceNodes[nodeIdx]] = vertex2AggId[interfaceNodes[fineNodeIdx]];
        procWinner[interfaceNodes[nodeIdx]]   = myRank;
        aggStat[interfaceNodes[nodeIdx]]      = AGGREGATED;
        --numNonAggregatedNodes;
      } // Loop over interface nodes
    } // Loop over the interfaces

    // Update aggregates information before subsequent aggregation algorithms
    aggregates->SetNumAggregates(aggregateCount);

    // Set coarse data for next level
    Set(currentLevel, "coarseInterfacesDimensions", coarseInterfacesDimensions);
    Set(currentLevel, "nodeOnCoarseInterface",      nodesOnCoarseInterfaces);

  } // BuildInterfaceAggregates()

} //namespace MueLu


#endif /* MUELU_HYBRIDAGGREGATIONFACTORY_DEF_HPP */