diff --git a/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp b/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp index ad179b1fb5d5..9bb82a44d103 100644 --- a/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp +++ b/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp @@ -59,6 +59,9 @@ #include "MueLu_SingleLevelFactoryBase.hpp" #include "MueLu_BrickAggregationFactory_fwd.hpp" +#include "MueLu_GraphBase.hpp" +#include "MueLu_Graph_fwd.hpp" +#include "MueLu_LWGraph_fwd.hpp" #include "MueLu_Level_fwd.hpp" #include "MueLu_Aggregates_fwd.hpp" #include "MueLu_Exceptions.hpp" @@ -66,11 +69,10 @@ /*! @class BrickAggregationFactory - @brief Aggregation method for generating "brick" aggregates. + @brief Aggregation method for generating "brick" aggregates. It also does "hotdogs" and "pancakes." - This factory can generate: - - In 2D, 2x2 and 3x3 aggregates - - In 3D, 2x2x2 or 3x3x3 aggregates + This factory can generate aggregates of size 1, 2 or 3 in each dimension, in any combination. + */ namespace MueLu { @@ -136,6 +138,9 @@ namespace MueLu { void Setup(const RCP >& comm, const RCP::magnitudeType,LO,GO,NO> >& coords, const RCP& map) const; RCP Construct1DMap(const RCP >& comm, const ArrayRCP::magnitudeType>& x) const; + void BuildGraph(Level& currentLevel, const RCP& A) const; + + bool isDirichlet(LocalOrdinal LID) const; bool isRoot (LocalOrdinal LID) const; GlobalOrdinal getRoot (LocalOrdinal LID) const; diff --git a/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp b/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp index 9f8c1b74a4ae..1b9aafcf3ae8 100644 --- a/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp +++ b/packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp @@ -61,13 +61,15 @@ #include #include - #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 { @@ -86,9 +88,6 @@ namespace MueLu { validParamList->set< RCP >("A", Teuchos::null, "Generating factory for matrix"); validParamList->set< RCP >("Coordinates", Teuchos::null, "Generating factory for coordinates"); - validParamList->set< RCP >("Graph" , Teuchos::null, "Generating factory of the graph"); - validParamList->set< RCP >("DofsPerNode", null, "Generating factory for variable \'DofsPerNode\', usually the same as for \'Graph\'"); - return validParamList; } @@ -192,6 +191,9 @@ namespace MueLu { << (nDim_ > 1 ? "x " + toString(by_) : "") << (nDim_ > 2 ? "x " + toString(bz_) : "") << std::endl; + // Build the graph + BuildGraph(currentLevel,A); + // Construct aggregates RCP aggregates = rcp(new Aggregates(colMap)); aggregates->setObjectLabel("Brick"); @@ -480,6 +482,106 @@ namespace MueLu { } + template + void BrickAggregationFactory::BuildGraph(Level& currentLevel, const RCP& 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 graph = rcp(new Graph(A->getCrsGraph(), "graph of A")); + ArrayRCP boundaryNodes = Teuchos::arcp_const_cast(MueLu::Utilities::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 > 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 rows (A->getLocalNumRows()+1); + ArrayRCP 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; rowgetColMap()->getLocalElement(A->getRowMap()->getGlobalElement(row)); + getIJK(row2,ir,jr,kr); + + for(size_t cidx=rowptr[row]; cidx graph = rcp(new LWGraph(rows, columns, A->getRowMap(), A->getColMap(), "thresholded graph of A")); + + + ArrayRCP boundaryNodes = Teuchos::arcp_const_cast(MueLu::Utilities::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 > 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_ */ diff --git a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp index 6d8586cf953d..8f553bf879b0 100644 --- a/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp +++ b/packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp @@ -1157,6 +1157,10 @@ namespace MueLu { MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: brick z Dirichlet", bool, aggParams); aggFactory->SetParameterList(aggParams); + // Unlike other factories, BrickAggregationFactory makes the Graph/DofsPerNode itself + manager.SetFactory("Graph", aggFactory); + manager.SetFactory("DofsPerNode", aggFactory); + manager.SetFactory("Filtering", aggFactory); if (levelID > 1) { // We check for levelID > 0, as in the interpreter aggFactory for // levelID really corresponds to level 0. Managers are clunky, as they