MueLu_BrickAggregationFactory_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_BRICKAGGREGATIONFACTORY_DEF_HPP_
#define MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_

#include "MueLu_BrickAggregationFactory_decl.hpp"
#ifdef HAVE_MPI
#include <Teuchos_DefaultMpiComm.hpp>
#include <Teuchos_CommHelpers.hpp>
#endif
#include <Teuchos_OrdinalTraits.hpp>

#include <Xpetra_Import.hpp>
#include <Xpetra_ImportFactory.hpp>
#include <Xpetra_Map.hpp>
#include <Xpetra_MapFactory.hpp>
#include <Xpetra_Matrix.hpp>
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>

#include "MueLu_Aggregates.hpp"
#include "MueLu_Level.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_GraphBase.hpp"
#include "MueLu_Graph.hpp"
#include "MueLu_LWGraph.hpp"


namespace MueLu {

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

#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
    SET_VALID_ENTRY("aggregation: brick x size");
    SET_VALID_ENTRY("aggregation: brick y size");
    SET_VALID_ENTRY("aggregation: brick z size");
    SET_VALID_ENTRY("aggregation: brick x Dirichlet");
    SET_VALID_ENTRY("aggregation: brick y Dirichlet");
    SET_VALID_ENTRY("aggregation: brick z Dirichlet");
#undef  SET_VALID_ENTRY

    validParamList->set< RCP<const FactoryBase> >("A",           Teuchos::null, "Generating factory for matrix");
    validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coordinates");
    return validParamList;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
    Input(currentLevel, "A");
    Input(currentLevel, "Coordinates");
  }

  // The current implementation cannot deal with bricks larger than 3x3(x3) in
  // parallel. The reason is that aggregation infrastructure in place has
  // major drawbacks.
  //
  // Aggregates class is constructed with a help of a provided map, either
  // taken from a graph, or provided directly. This map is usually taken to be
  // a column map of a matrix. The reason for that is that if we have an
  // overlapped aggregation, we want the processor owning aggregates to store
  // agg id for all nodes in this aggregate. If we used row map, there would
  // be no way for the processor to know whether there are some other nodes on
  // a different processor which belong to its aggregate. On the other hand,
  // using column map allows both vertex2AggId and procWinner arrays in
  // Aggregates class to store some extra data, such as whether nodes belonging
  // to a different processor belong to this processor aggregate.
  //
  // The drawback of this is that it stores only overlap=1 data. For aggressive
  // coarsening, such a brick aggregation with a large single dimension of
  // brick, it could happen that we need to know depth two or more extra nodes
  // in the other processor subdomain.
  //
  // Another issue is that we may have some implicit connection between
  // aggregate map and maps of A used in the construction of a tentative
  // prolongator.
  //
  // Another issue is that it seems that some info is unused or not required.
  // Specifically, it seems that if a node belongs to an aggregate on a
  // different processor, we don't actually need to set vertex2AggId and
  // procWinner, despite the following comment in
  // Aggregates decl:
  //      vertex2AggId[k] gives a local id
  //      corresponding to the aggregate to which
  //      local id k has been assigned.  While k
  //      is the local id on my processor (MyPID)
  //      vertex2AggId[k] is the local id on the
  //      processor which actually owns the
  //      aggregate. This owning processor has id
  //      given by procWinner[k].
  // It is possible that that info is only used during arbitration in
  // CoupledAggregationFactory.
  //
  // The steps that we need to do to resolve this issue:
  //   - Break the link between maps in TentativePFactory, allowing any maps in Aggregates
  //   - Allow Aggregates to construct their own maps, if necessary, OR
  //   - construct aggregates based on row map
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& currentLevel) const {
    FactoryMonitor m(*this, "Build", currentLevel);

    typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> MultiVector_d;

    const ParameterList& pL = GetParameterList();
    RCP<MultiVector_d> coords = Get<RCP<MultiVector_d> >(currentLevel, "Coordinates");
    RCP<Matrix>        A      = Get< RCP<Matrix> >      (currentLevel, "A");
    RCP<const Map>     rowMap = A->getRowMap();
    RCP<const Map>     colMap = A->getColMap();
    GO GO_INVALID = Teuchos::OrdinalTraits<GO>::invalid();

    RCP<const Teuchos::Comm<int> > comm = rowMap->getComm();
    int numProcs = comm->getSize();
    int myRank   = comm->getRank();

    int numPoints = colMap->getLocalNumElements();

    bx_ = pL.get<int>("aggregation: brick x size");
    by_ = pL.get<int>("aggregation: brick y size");
    bz_ = pL.get<int>("aggregation: brick z size");

    dirichletX_ = pL.get<bool>("aggregation: brick x Dirichlet");
    dirichletY_ = pL.get<bool>("aggregation: brick y Dirichlet");
    dirichletZ_ = pL.get<bool>("aggregation: brick z Dirichlet");
    if(dirichletX_) GetOStream(Runtime0) << "Dirichlet boundaries in the x direction"<<std::endl;
    if(dirichletY_) GetOStream(Runtime0) << "Dirichlet boundaries in the y direction"<<std::endl;
    if(dirichletZ_) GetOStream(Runtime0) << "Dirichlet boundaries in the z direction"<<std::endl;

    if (numProcs > 1) {
      // TODO: deal with block size > 1  (see comments above)
      //TEUCHOS_TEST_FOR_EXCEPTION(bx_ > 3 || by_ > 3 || bz_ > 3, Exceptions::RuntimeError, "Currently cannot deal with brick size > 3");
    }

    RCP<MultiVector_d> overlappedCoords = coords;
    RCP<const Import> importer = ImportFactory::Build(coords->getMap(), colMap);
    if (!importer.is_null()) {
      overlappedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(colMap, coords->getNumVectors());
      overlappedCoords->doImport(*coords, *importer, Xpetra::INSERT);
    }

    // Setup misc structures
    // Logically, we construct enough data to query topological information of a rectangular grid
    Setup(comm, overlappedCoords, colMap);

    GetOStream(Runtime0) << "Using brick size: " << bx_
                                                 << (nDim_ > 1 ? "x " + toString(by_) : "")
                                                 << (nDim_ > 2 ? "x " + toString(bz_) : "") << std::endl;

    // Build the graph
    BuildGraph(currentLevel,A);

    // Construct aggregates
    RCP<Aggregates> aggregates = rcp(new Aggregates(colMap));
    aggregates->setObjectLabel("Brick");

    ArrayRCP<LO> vertex2AggId = aggregates->GetVertex2AggId()->getDataNonConst(0);
    ArrayRCP<LO> procWinner   = aggregates->GetProcWinner()  ->getDataNonConst(0);

    // In the first pass, we set a mapping from a vertex to aggregate global id. We deal with a structured
    // rectangular mesh, therefore we know the structure of aggregates. For each vertex we can tell exactly
    // which aggregate it belongs to.
    // If we determine that the aggregate does not belong to us (i.e. the root vertex does not belong to this
    // processor, or is outside and we lost "" arbitration), we record the global aggregate id in order to
    // fetch the local info from the processor owning the aggregate. This is required for aggregates, as it
    // uses the local aggregate ids of the owning processor.
    std::set<GO> myAggGIDs, remoteAggGIDs;
    for (LO LID = 0; LID < numPoints; LID++) {
      GO aggGID = getAggGID(LID);
      //      printf("[%d] (%d,%d,%d) => agg %d\n",LID,(int)(*xMap_)[x_[LID]],nDim_ > 1 ? (int)(*yMap_)[y_[LID]] : -1,nDim_ > 2 ? (int)(*zMap_)[z_[LID]] : -1,(int)aggGID);
      if(aggGID == GO_INVALID) continue;
      //      printf("[%d] getRoot = %d\n",(int)LID,(int)getRoot(LID));

      if ((revMap_.find(getRoot(LID)) != revMap_.end()) && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
        // Root of the brick aggregate containing GID (<- LID) belongs to us
        vertex2AggId[LID] = aggGID;
        myAggGIDs.insert(aggGID);

        if (isRoot(LID))
          aggregates->SetIsRoot(LID);
  //  printf("[%d] initial vertex2AggId = %d\n",(int)LID,(int)vertex2AggId[LID]);
      } else {
        remoteAggGIDs.insert(aggGID);
      }
    }
    size_t numAggregates = myAggGIDs    .size();
    size_t numRemote     = remoteAggGIDs.size();
    aggregates->SetNumAggregates(numAggregates);

    std::map<GO,LO>  AggG2L;  // Map: Agg GID -> Agg LID (possibly on a different processor)
    std::map<GO,int> AggG2R;  // Map: Agg GID -> processor rank owning aggregate

    Array<GO> myAggGIDsArray(numAggregates), remoteAggGIDsArray(numRemote);

    // Fill in the maps for aggregates that we own
    size_t ind = 0;
    for (typename std::set<GO>::const_iterator it = myAggGIDs.begin(); it != myAggGIDs.end(); it++) {
      AggG2L[*it] = ind;
      AggG2R[*it] = myRank;

      myAggGIDsArray[ind++] = *it;
    }

    // The map is a convenient way to fetch remote local indices from global indices.
    RCP<Map> aggMap = MapFactory::Build(rowMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
                                        myAggGIDsArray, 0, comm);

    ind = 0;
    for (typename std::set<GO>::const_iterator it = remoteAggGIDs.begin(); it != remoteAggGIDs.end(); it++)
      remoteAggGIDsArray[ind++] = *it;

    // Fetch the required aggregate local ids and ranks
    Array<int> remoteProcIDs(numRemote);
    Array<LO>  remoteLIDs   (numRemote);
    aggMap->getRemoteIndexList(remoteAggGIDsArray, remoteProcIDs, remoteLIDs);

    // Fill in the maps for aggregates that we don't own but which have some of our vertices
    for (size_t i = 0; i < numRemote; i++) {
      AggG2L[remoteAggGIDsArray[i]] = remoteLIDs   [i];
      AggG2R[remoteAggGIDsArray[i]] = remoteProcIDs[i];
    }

    // Remap aggregate GIDs to LIDs and set up owning processors
    for (LO LID = 0; LID < numPoints; LID++) {
      if (revMap_.find(getRoot(LID)) != revMap_.end() && rowMap->isNodeGlobalElement(colMap->getGlobalElement(revMap_[getRoot(LID)]))) {
        GO aggGID = vertex2AggId[LID];
        if(aggGID != MUELU_UNAGGREGATED) {
          vertex2AggId[LID] = AggG2L[aggGID];
          procWinner  [LID] = AggG2R[aggGID];
        }
      }
    }


    GO numGlobalRemote;
    MueLu_sumAll(comm, as<GO>(numRemote), numGlobalRemote);
    aggregates->AggregatesCrossProcessors(numGlobalRemote);

    Set(currentLevel, "Aggregates", aggregates);

    GetOStream(Statistics1) << aggregates->description() << std::endl;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  Setup(const RCP<const Teuchos::Comm<int> >& comm, const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> >& coords, const RCP<const Map>& /* map */) const {
    nDim_ = coords->getNumVectors();

    x_    = coords->getData(0);
    xMap_ = Construct1DMap(comm, x_);
    nx_   = xMap_->size();

    ny_   = 1;
    if (nDim_ > 1) {
      y_    = coords->getData(1);
      yMap_ = Construct1DMap(comm, y_);
      ny_   = yMap_->size();
    }

    nz_   = 1;
    if (nDim_ > 2) {
      z_    = coords->getData(2);
      zMap_ = Construct1DMap(comm, z_);
      nz_   = zMap_->size();
    }

    for (size_t ind = 0; ind < coords->getLocalLength(); ind++) {
      GO i = (*xMap_)[(coords->getData(0))[ind]], j = 0, k = 0;
      if (nDim_ > 1)
        j = (*yMap_)[(coords->getData(1))[ind]];
      if (nDim_ > 2)
        k = (*zMap_)[(coords->getData(2))[ind]];

      revMap_[k*ny_*nx_ + j*nx_ + i] = ind;
    }

   
    // Get the number of aggregates in each direction, correcting for Dirichlet
    int xboost = dirichletX_ ? 1 : 0;
    int yboost = dirichletY_ ? 1 : 0;
    int zboost = dirichletZ_ ? 1 : 0;
    naggx_ = (nx_-2*xboost)/bx_ + ((nx_-2*xboost) % bx_ ? 1 : 0);

    if(nDim_ > 1)
      naggy_ = (ny_-2*yboost)/by_ + ( (ny_-2*yboost) % by_ ? 1 : 0);
    else
      naggy_ = 1;

    if(nDim_ > 2)
      naggz_ = (nz_-2*zboost)/bz_ + ( (nz_-2*zboost) % bz_ ? 1 : 0);
    else
      naggz_ = 1;

  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  RCP<typename BrickAggregationFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::container>
  BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  Construct1DMap (const RCP<const Teuchos::Comm<int> >& comm,
                  const ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x) const
  {
    int n = x.size();

    // Step 1: Create a local vector with unique coordinate points
    RCP<container> gMap = rcp(new container);
    for (int i = 0; i < n; i++)
      (*gMap)[x[i]] = 0;

#ifdef HAVE_MPI
    // Step 2: exchange coordinates
    // NOTE: we assume the coordinates are double, or double compatible
    // That means that for complex case, we assume that all imaginary parts are zeros
    int numProcs = comm->getSize();
    if (numProcs > 1) {
      RCP<const Teuchos::MpiComm<int> > dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm->duplicate());

      MPI_Comm rawComm = (*dupMpiComm->getRawMpiComm())();

      int           sendCnt = gMap->size(), cnt = 0, recvSize;
      Array<int>    recvCnt(numProcs), Displs(numProcs);
      Array<double> sendBuf, recvBuf;

      sendBuf.resize(sendCnt);
      for (typename container::const_iterator cit = gMap->begin(); cit != gMap->end(); cit++)
        sendBuf[cnt++] = Teuchos::as<double>(STS::real(cit->first));

      MPI_Allgather(&sendCnt, 1, MPI_INT, recvCnt.getRawPtr(), 1, MPI_INT, rawComm);
      Displs[0] = 0;
      for (int i = 0; i < numProcs-1; i++)
        Displs[i+1] = Displs[i] + recvCnt[i];
      recvSize = Displs[numProcs-1] + recvCnt[numProcs-1];
      recvBuf.resize(recvSize);
      MPI_Allgatherv(sendBuf.getRawPtr(), sendCnt, MPI_DOUBLE, recvBuf.getRawPtr(), recvCnt.getRawPtr(), Displs.getRawPtr(), MPI_DOUBLE, rawComm);

      for (int i = 0; i < recvSize; i++)
        (*gMap)[as<SC>(recvBuf[i])] = 0;
    }
#endif

    GO cnt = 0;
    for (typename container::iterator it = gMap->begin(); it != gMap->end(); it++)
      it->second = cnt++;

    return gMap;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  bool BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isRoot(LocalOrdinal LID) const {
    int i,j,k;
    getIJK(LID,i,j,k);

    return (k*ny_*nx_ + j*nx_ + i) == getRoot(LID);
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  bool BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isDirichlet(LocalOrdinal LID) const {
    bool boundary = false;
    int i,j,k;
    getIJK(LID,i,j,k);
    if(             dirichletX_  &&  (i == 0 || i == nx_-1) )
      boundary = true;
    if(nDim_ > 1 && dirichletY_  &&  (j == 0 || j == ny_-1) )
      boundary = true;
    if(nDim_ > 2 && dirichletZ_  &&  (k == 0 || k == nz_-1) )
      boundary = true;
    
    return boundary;
  }


  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  GlobalOrdinal BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getRoot(LocalOrdinal LID) const {
    if(isDirichlet(LID))
      return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();

    int aggI,aggJ,aggK;
    getAggIJK(LID,aggI,aggJ,aggK);
    int xboost = dirichletX_ ? 1 : 0;
    int yboost = dirichletY_ ? 1 : 0;
    int zboost = dirichletZ_ ? 1 : 0;
   
    int i =             xboost + aggI*bx_ + (bx_-1)/2;
    int j = (nDim_>1) ? yboost + aggJ*by_ + (by_-1)/2 : 0;
    int k = (nDim_>2) ? zboost + aggK*bz_ + (bz_-1)/2 : 0;

    return k*ny_*nx_ + j*nx_ + i;
  }

 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void  BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getIJK(LocalOrdinal LID, int &i, int &j, int &k) const {
   i = (*xMap_)[x_[LID]];
   j = (nDim_>1) ? (*yMap_)[y_[LID]] : 0;
   k =  (nDim_>2) ? (*zMap_)[z_[LID]] : 0;
  }


  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void  BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getAggIJK(LocalOrdinal LID, int &i, int &j, int &k) const {
    int xboost = dirichletX_ ? 1 : 0;
    int yboost = dirichletY_ ? 1 : 0;
    int zboost = dirichletZ_ ? 1 : 0;
    int pointI, pointJ, pointK;
    getIJK(LID,pointI,pointJ,pointK);
    i = (pointI-xboost)/bx_;

    if (nDim_ > 1) j = (pointJ-yboost)/by_;
    else j = 0;

    if (nDim_ > 2) k = (pointK-zboost)/bz_;
    else k = 0;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  GlobalOrdinal BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getAggGID(LocalOrdinal LID) const {
    bool boundary = false;

    int i, j, k;
    getIJK(LID,i,j,k);
    int ii , jj, kk;
    getAggIJK(LID,ii,jj,kk);

    if(              dirichletX_ && (i == 0 || i == nx_ - 1)) boundary = true;
    if (nDim_ > 1 && dirichletY_ && (j == 0 || j == ny_ - 1)) boundary = true;
    if (nDim_ > 2 && dirichletZ_ && (k == 0 || k == nz_ - 1)) boundary = true;

    /*
    if(boundary)
      printf("[%d] coord = (%d,%d,%d) {%d,%d,%d} agg = (%d,%d,%d) {%d,%d,%d} => agg %s\n",LID,i,j,k,nx_,ny_,nz_,ii,jj,kk,naggx_,naggy_,naggz_,"BOUNDARY");
    else
      printf("[%d] coord = (%d,%d,%d) {%d,%d,%d} agg = (%d,%d,%d) {%d,%d,%d} => agg %d\n",LID,i,j,k,nx_,ny_,nz_,ii,jj,kk,naggx_,naggy_,naggz_,kk*naggy_*naggx_ + jj*naggx_ + ii);
    */

    if (boundary)
      return Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
    else
      return Teuchos::as<GlobalOrdinal>(kk*naggy_*naggx_) + Teuchos::as<GlobalOrdinal>(jj*naggx_) + ii;
    
  }


  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void BrickAggregationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildGraph(Level& currentLevel, const RCP<Matrix>& A) const {
    // TODO: Currently only works w/ 1 DOF per node
    double dirichletThreshold = 0.0;

    if(bx_ > 1 && (nDim_ <= 1 || by_ > 1) && (nDim_ <=2 || bz_>1) ) {
      FactoryMonitor m(*this, "Generating Graph (trivial)", currentLevel);
      /*** Case 1: Use the matrix is the graph ***/
      // Bricks are of non-trivial size in all active dimensions
      RCP<GraphBase> graph = rcp(new Graph(A->getCrsGraph(), "graph of A"));
      ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
      graph->SetBoundaryNodeMap(boundaryNodes);
      
      if (GetVerbLevel() & Statistics1) {
        GO numLocalBoundaryNodes  = 0;
        GO numGlobalBoundaryNodes = 0;
        for (LO i = 0; i < boundaryNodes.size(); ++i)
          if (boundaryNodes[i])
            numLocalBoundaryNodes++;
        RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
        MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
        GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
      }
      Set(currentLevel, "DofsPerNode", 1);
      Set(currentLevel, "Graph", graph);
      Set(currentLevel, "Filtering",false);
    }
    else {
      FactoryMonitor m(*this, "Generating Graph", currentLevel);
      /*** Case 2: Dropping required ***/
      // There is at least one active dimension in which we are not coarsening.
      // Those connections need to be dropped
      bool drop_x = (bx_ == 1);
      bool drop_y = (nDim_> 1 && by_ == 1);
      bool drop_z = (nDim_> 2 && bz_ == 1);
      
      ArrayRCP<LO> rows   (A->getLocalNumRows()+1);
      ArrayRCP<LO> columns(A->getLocalNumEntries());
      
      size_t N = A->getRowMap()->getLocalNumElements();

      // FIXME: Do this on the host because indexing functions are host functions
      auto G = A->getLocalMatrixHost().graph;
      auto rowptr = G.row_map;
      auto colind = G.entries;

      int ct=0;
      rows[0] = 0;
      for(size_t row=0; row<N; row++) {
        // NOTE: Assumes that the first part of the colmap is the rowmap
        int ir,jr,kr;
        LO row2 = A->getColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row));
        getIJK(row2,ir,jr,kr);

        for(size_t cidx=rowptr[row]; cidx<rowptr[row+1]; cidx++) {
          int ic,jc,kc;
          LO col = colind[cidx];
          getIJK(col,ic,jc,kc);

          if( (row2 !=col) && ((drop_x && ir != ic) || (drop_y && jr != jc) || (drop_z && kr != kc) )) {
            // Drop it
            //            printf("[%4d] DROP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
          }
          else {
            // Keep it
            //            printf("[%4d] KEEP row = (%d,%d,%d) col = (%d,%d,%d)\n",(int)row,ir,jr,kr,ic,jc,kc);
            columns[ct] = col;
            ct++;
          }
        }
        rows[row+1] = ct;
      }//end for

      RCP<GraphBase> graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A"));


      ArrayRCP<bool > boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
      graph->SetBoundaryNodeMap(boundaryNodes);
      
      if (GetVerbLevel() & Statistics1) {
        GO numLocalBoundaryNodes  = 0;
        GO numGlobalBoundaryNodes = 0;
        for (LO i = 0; i < boundaryNodes.size(); ++i)
          if (boundaryNodes[i])
            numLocalBoundaryNodes++;
        RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
        MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
        GetOStream(Statistics1) << "Detected " << numGlobalBoundaryNodes << " Dirichlet nodes" << std::endl;
      }
      Set(currentLevel, "DofsPerNode", 1);
      Set(currentLevel, "Graph", graph);
      Set(currentLevel, "Filtering",true);
    }//end else


  }//end BuildGraph

  


} //namespace MueLu

#endif /* MUELU_BRICKAGGREGATIONFACTORY_DEF_HPP_ */