MueLu: Adding filtering support for BrickAggregation

packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_decl.hpp

@@ -59,18 +59,20 @@
 #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"
 #include "MueLu_Utilities_fwd.hpp"
 
 /*!
   @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<const Teuchos::Comm<int> >& comm, const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> >& coords, const RCP<const Map>& map) const;
     RCP<container> Construct1DMap(const RCP<const Teuchos::Comm<int> >& comm, const ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& x) const;
 
+    void BuildGraph(Level& currentLevel, const RCP<Matrix>& A) const;
+
+
     bool        isDirichlet(LocalOrdinal LID) const;
     bool           isRoot  (LocalOrdinal LID) const;
     GlobalOrdinal getRoot  (LocalOrdinal LID) const;

packages/muelu/src/Graph/BrickAggregation/MueLu_BrickAggregationFactory_def.hpp

@@ -61,13 +61,15 @@
 #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 {
 
@@ -86,9 +88,6 @@ namespace MueLu {
 
     validParamList->set< RCP<const FactoryBase> >("A",           Teuchos::null, "Generating factory for matrix");
     validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coordinates");
-    validParamList->set< RCP<const FactoryBase> >("Graph"      , Teuchos::null, "Generating factory of the graph");
-    validParamList->set< RCP<const FactoryBase> >("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> aggregates = rcp(new Aggregates(colMap));
     aggregates->setObjectLabel("Brick");
@@ -480,6 +482,106 @@ namespace MueLu {
   }
 
 
+  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_ */

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

packages/muelu/test/interface/default/Output/aggregation5_epetra.gold

@@ -8,16 +8,14 @@ smoother ->
 Level 1
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)
@@ -42,16 +40,14 @@ smoother ->
 Level 2
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)

packages/muelu/test/interface/default/Output/aggregation5_np4_epetra.gold

@@ -8,16 +8,14 @@ smoother ->
 Level 1
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)
@@ -42,16 +40,14 @@ smoother ->
 Level 2
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)

packages/muelu/test/interface/default/Output/aggregation5_np4_tpetra.gold

@@ -8,16 +8,14 @@ smoother ->
 Level 1
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)
@@ -42,16 +40,14 @@ smoother ->
 Level 2
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)

packages/muelu/test/interface/default/Output/aggregation5_tpetra.gold

@@ -8,16 +8,14 @@ smoother ->
 Level 1
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)
@@ -42,16 +40,14 @@ smoother ->
 Level 2
 Prolongator smoothing (MueLu::SaPFactory)
 Matrix filtering (MueLu::FilteredAFactory)
-Build (MueLu::CoalesceDropFactory)
-Build (MueLu::AmalgamationFactory)
-[empty list]
-algorithm = "classical" classical algorithm = "default": threshold = 0, blocksize = 1
-lightweight wrap = 1
-Filtered matrix is not being constructed as no filtering is being done
-Build (MueLu::TentativePFactory)
 Build (MueLu::BrickAggregationFactory)
 Using brick size: 3
+Generating Graph (trivial) (MueLu::BrickAggregationFactory)
 aggregation: brick x size = 3
+Filtered matrix is not being constructed as no filtering is being done
+Build (MueLu::TentativePFactory)
+Build (MueLu::AmalgamationFactory)
+[empty list]
 Nullspace factory (MueLu::NullspaceFactory)
 Fine level nullspace = Nullspace
 Build (MueLu::CoarseMapFactory)