MueLu: More alternative coarsening

packages/muelu/doc/UsersGuide/masterList.xml

@@ -624,6 +624,17 @@
       <comment-ML>parameter not existing in ML</comment-ML>
     </parameter>
 
+
+    <parameter>
+      <name>aggregation: coloring: use color graph</name>
+      <type>bool</type>
+      <default>false</default>
+      <visible>false</visible>
+      <description>Have CoalesceDropFactory generate a seperate graph for SOC and for coloring</description>
+      <comment-ML>parameter not existing in ML</comment-ML>
+    </parameter>
+
+
     <parameter>
       <name>aggregation: enable phase 1</name>
       <type>bool</type>

packages/muelu/doc/UsersGuide/paramlist_hidden.tex

@@ -122,6 +122,8 @@
 
 \cbb{aggregation: coloring algorithm}{string}{serial}{Choice of distance 2 coloring algorithm used by Uncoupled Aggregation. Currently set by default to COLORING\_D2\_SERIAL.}
 
+\cbb{aggregation: coloring: use color graph}{bool}{false}{Have CoalesceDropFactory generate a seperate graph for SOC and for coloring}
+
 \cbb{aggregation: enable phase 1}{bool}{true}{Turn on/off phase 1 of aggregation}
 
 \cbb{aggregation: enable phase 2a}{bool}{true}{Turn on/off phase 2a of aggregation}

packages/muelu/src/Graph/Containers/MueLu_LWGraph_def.hpp

@@ -55,7 +55,7 @@ namespace MueLu {
   //void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel = Default) const {
   template <class LocalOrdinal, class GlobalOrdinal, class Node>
   void LWGraph<LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
-    MUELU_DESCRIBE;
+    //    MUELU_DESCRIBE;
 
     if (verbLevel & Parameters0) {
       //out0 << "Prec. type: " << type_ << std::endl;
@@ -69,9 +69,11 @@ namespace MueLu {
     }
 
     if (verbLevel & Debug) {
+      RCP<const Map> col_map = importMap_.is_null() ? domainMap_ : importMap_;
+
       for (LO i = 0; i < rows_.size()-1; i++) {
         for (LO j = rows_[i]; j < rows_[i+1]; j++)
-          out0 << i << " " << columns_[j]<<std::endl;
+          out<< domainMap_->getGlobalElement(i) << " " << col_map->getGlobalElement(columns_[j])<<std::endl;
       }
     }
   }

packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_decl.hpp

@@ -178,6 +178,9 @@ namespace MueLu {
     // When we want to decouple a block diagonal system (returns Teuchos::null if generate_matrix is false)
     Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > BlockDiagonalize(Level & currentLevel,const RCP<Matrix> & A, bool generate_matrix) const;
 
+    // When we want to decouple a block diagonal system via a *graph*
+    void BlockDiagonalizeGraph(const RCP<GraphBase> & inputGraph, const RCP<LocalOrdinalVector> & ghostedBlockNumber, RCP<GraphBase> & outputGraph) const;
+
   }; //class CoalesceDropFactory
 
 } //namespace MueLu

packages/muelu/src/Graph/MatrixTransformation/MueLu_CoalesceDropFactory_def.hpp

@@ -125,7 +125,7 @@ namespace MueLu {
 
     {
       typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
-      validParamList->getEntry("aggregation: drop scheme").setValidator(rcp(new validatorType(Teuchos::tuple<std::string>("classical", "distance laplacian","signed classical","block diagonal","block diagonal classical","block diagonal distance laplacian","block diagonal signed classical"), "aggregation: drop scheme")));
+      validParamList->getEntry("aggregation: drop scheme").setValidator(rcp(new validatorType(Teuchos::tuple<std::string>("classical", "distance laplacian","signed classical","block diagonal","block diagonal classical","block diagonal distance laplacian","block diagonal signed classical","block diagonal colored signed classical"), "aggregation: drop scheme")));
 
     }
     SET_VALID_ENTRY("aggregation: distance laplacian algo");
@@ -155,8 +155,7 @@ namespace MueLu {
       if (algo == "distance laplacian" || algo == "block diagonal distance laplacian") {
         Input(currentLevel, "Coordinates");    
       }
-      if (algo == "block diagonal classical" || algo == "block diagonal distance laplacian" 
-          || algo == "block diagonal" || algo == "block diagonal signed classical")  {
+      if (algo.find("block diagonal") != std::string::npos || algo.find("signed classical") != std::string::npos)  {
         Input(currentLevel, "BlockNumber");
       }
     }     
@@ -176,7 +175,6 @@ namespace MueLu {
     if (predrop_ != Teuchos::null)
       GetOStream(Parameters0) << predrop_->description();
 
-
     RCP<Matrix> realA = Get< RCP<Matrix> >(currentLevel, "A");
     RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel, "UnAmalgamationInfo");
     const ParameterList  & pL = GetParameterList();
@@ -191,27 +189,47 @@ namespace MueLu {
     bool use_block_algorithm=false;
     LO interleaved_blocksize = as<LO>(pL.get<int>("aggregation: block diagonal: interleaved blocksize"));
     bool useSignedClassical = false;
+    bool generateColoringGraph = false;
 
     // NOTE:  If we're doing blockDiagonal, we'll not want to do rowSum twice (we'll do it
     // in the block diagonalizaiton). So we'll clobber the rowSumTol with -1.0 in this case
     typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<double>("aggregation: row sum drop tol"));
+    RCP<LocalOrdinalVector> ghostedBlockNumber;
+    ArrayRCP<const LO> g_block_id;
 
     if(algo == "distance laplacian" ) { 
       // Grab the coordinates for distance laplacian
       Coords = Get< RCP<RealValuedMultiVector > >(currentLevel, "Coordinates");
       A = realA;
     }
-    else if(algo == "signed classical") {
+    else if(algo == "signed classical" || algo == "block diagonal colored signed classical" || algo == "block diagonal signed classical") {
       useSignedClassical = true;
+      //      if(realA->GetFixedBlockSize() > 1) {
+        RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel, "BlockNumber");
+        // Ghost the column block numbers if we need to
+        RCP<const Import> importer = realA->getCrsGraph()->getImporter();
+        if(!importer.is_null()) {
+          SubFactoryMonitor m1(*this, "Block Number import", currentLevel);      
+          ghostedBlockNumber= Xpetra::VectorFactory<LO,LO,GO,NO>::Build(importer->getTargetMap());
+          ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);          
+        } 
+        else {
+          ghostedBlockNumber = BlockNumber;
+        }
+        g_block_id = ghostedBlockNumber->getData(0);
+        //      }
+      if(algo == "block diagonal colored signed classical") 
+        generateColoringGraph=true;
       algo = "classical";
       A = realA;
+
     }
     else if(algo == "block diagonal") {
       // Handle the "block diagonal" filtering and then leave
       BlockDiagonalize(currentLevel,realA,false);
       return;
     }
-    else if (algo == "block diagonal classical" || algo == "block diagonal distance laplacian" || algo == "block diagonal signed classical")  {
+    else if (algo == "block diagonal classical" || algo == "block diagonal distance laplacian")  {
       // Handle the "block diagonal" filtering, and then continue onward
       use_block_algorithm = true;
       RCP<Matrix> filteredMatrix = BlockDiagonalize(currentLevel,realA,true);
@@ -239,10 +257,6 @@ namespace MueLu {
       else if(algo == "block diagonal classical") {
         algo = "classical";
       }
-      else if(algo == "block diagonal signed classical") {
-        algo = "classical";
-        useSignedClassical = true;
-      }
       // All cases
       A = filteredMatrix;
       rowSumTol = -1.0;
@@ -423,20 +437,34 @@ namespace MueLu {
           ArrayRCP<const SC> ghostedDiagVals;
           ArrayRCP<const MT> negMaxOffDiagonal;
           if(useSignedClassical) {
-            negMaxOffDiagonal = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixMaxMinusOffDiagonal(*A);
+            if(ghostedBlockNumber.is_null()) {
+              if (GetVerbLevel() & Statistics1)
+                GetOStream(Statistics1) << "Calculating max point off-diagonal"<<std::endl;
+              negMaxOffDiagonal = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixMaxMinusOffDiagonal(*A);
+            }
+            else {
+              if (GetVerbLevel() & Statistics1)
+                GetOStream(Statistics1) << "Calculating max block off-diagonal"<<std::endl;
+              negMaxOffDiagonal = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixMaxMinusOffDiagonal(*A,*ghostedBlockNumber);
+            }
           }
           else {
             ghostedDiag = MueLu::Utilities<SC,LO,GO,NO>::GetMatrixOverlappedDiagonal(*A);
             ghostedDiagVals = ghostedDiag->getData(0);
           }
           ArrayRCP<bool> boundaryNodes = Teuchos::arcp_const_cast<bool>(MueLu::Utilities<SC,LO,GO,NO>::DetectDirichletRows(*A, dirichletThreshold));
-          if (rowSumTol > 0.)
-            Utilities::ApplyRowSumCriterion(*A, rowSumTol, boundaryNodes);          
+          if (rowSumTol > 0.) {
+            if(ghostedBlockNumber.is_null())
+              Utilities::ApplyRowSumCriterion(*A, rowSumTol, boundaryNodes);          
+            else
+              Utilities::ApplyRowSumCriterion(*A, *ghostedBlockNumber, rowSumTol, boundaryNodes);          
+          }
 
           LO realnnz = 0;
           rows[0] = 0;
           for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
             size_t nnz = A->getNumEntriesInLocalRow(row);
+            bool rowIsDirichlet = boundaryNodes[row];
             ArrayView<const LO> indices;
             ArrayView<const SC> vals;
             A->getLocalRowView(row, indices, vals);
@@ -454,8 +482,11 @@ namespace MueLu {
                   LO col = indices[colID];               
                   MT max_neg_aik = realThreshold * STS::real(negMaxOffDiagonal[row]);
                   MT neg_aij    = - STS::real(vals[colID]);
-                  //printf(" - a_ij = %6.4e >? %6.4e * %6.4e = alpha max(-aik)\n",neg_aij,threshold, negMaxOffDiagonal[row]);
-                  if (neg_aij > max_neg_aik || row == col) {
+                  /*                  if(row==1326) printf("A(%d,%d) = %6.4e, block = (%d,%d) neg_aij = %6.4e max_neg_aik = %6.4e\n",row,col,vals[colID],
+                                       g_block_id.is_null() ? -1 :  g_block_id[row],
+                                       g_block_id.is_null() ? -1 :  g_block_id[col],
+                                       neg_aij, max_neg_aik);*/
+                  if ((!rowIsDirichlet && (g_block_id.is_null() || g_block_id[row] == g_block_id[col]) && neg_aij > max_neg_aik) || row == col) {
                     columns[realnnz++] = col;
                     rownnz++;
                   } else
@@ -470,7 +501,7 @@ namespace MueLu {
                   MT  aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);  // eps^2*|a_ii|*|a_jj|
                   MT aij    = STS::magnitude(vals[colID]*vals[colID]);                                          // |a_ij|^2
 
-                  if (aij > aiiajj || row == col) {
+                  if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
                   columns[realnnz++] = col;
                   rownnz++;
                   } else
@@ -488,6 +519,7 @@ namespace MueLu {
               const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
               const real_type one  = Teuchos::ScalarTraits<real_type>::one();
               LO rownnz = 0;
+              // NOTE: This probably needs to be fixed for rowsum
 
               // find magnitudes
               for (LO colID = 0; colID < (LO)nnz; colID++) {
@@ -605,6 +637,29 @@ namespace MueLu {
           Set(currentLevel, "Graph",       graph);
           Set(currentLevel, "DofsPerNode", 1);
 
+          // If we're doing signed classical, we might want to block-diagonalize *after* the dropping
+          if(generateColoringGraph) {
+            RCP<GraphBase> colorGraph;
+            BlockDiagonalizeGraph(graph,ghostedBlockNumber,colorGraph);
+            Set(currentLevel, "Coloring Graph",colorGraph);
+            //#define CMS_DUMP
+#ifdef CMS_DUMP
+           {
+             int rank = graph->GetDomainMap()->getComm()->getRank();
+             {
+               std::ofstream ofs(std::string("m_color_graph_") + std::to_string(currentLevel.GetLevelID())+std::string("_") + std::to_string(rank) + std::string(".dat"),std::ofstream::out);
+               RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(ofs));
+               colorGraph->print(*fancy,Debug);
+             }
+             {
+               std::ofstream ofs(std::string("m_regular_graph_") + std::to_string(currentLevel.GetLevelID())+std::string("_") + std::to_string(rank) + std::string(".dat"),std::ofstream::out);
+               RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(ofs));
+               graph->print(*fancy,Debug);
+             }
+
+           }
+#endif
+          }//end generateColoringGraph
         } else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
           // Case 3:  Multiple DOF/node problem without dropping
           const RCP<const Map> rowMap = A->getRowMap();
@@ -1621,7 +1676,7 @@ namespace MueLu {
 
     RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel, "BlockNumber");
     RCP<LocalOrdinalVector> ghostedBlockNumber;
-    GetOStream(Statistics1) << "Using BlockDiagonal Graph (with provided blocking)"<<std::endl;      
+    GetOStream(Statistics1) << "Using BlockDiagonal Graph before dropping (with provided blocking)"<<std::endl;      
 
     // Ghost the column block numbers if we need to
     RCP<const Import> importer = A->getCrsGraph()->getImporter();
@@ -1730,6 +1785,71 @@ namespace MueLu {
     return crs_matrix_wrap;
   }
 
+
+ template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
+ void  CoalesceDropFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BlockDiagonalizeGraph(const RCP<GraphBase> & inputGraph, const RCP<LocalOrdinalVector> & ghostedBlockNumber, RCP<GraphBase> & outputGraph) const {
+typedef Teuchos::ScalarTraits<SC> STS;
+
+ TEUCHOS_TEST_FOR_EXCEPTION(ghostedBlockNumber.is_null(), Exceptions::RuntimeError, "BlockDiagonalizeGraph(): ghostedBlockNumber is null.");
+//    const ParameterList  & pL = GetParameterList();
+
+    GetOStream(Statistics1) << "Using BlockDiagonal Graph after Dropping (with provided blocking)"<<std::endl;      
+
+    // Accessors for block numbers
+    Teuchos::ArrayRCP<const LO> row_block_number = ghostedBlockNumber->getData(0);
+    Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
+
+    // allocate space for the local graph
+    ArrayRCP<size_t> rows_mat;   
+    ArrayRCP<LO> rows_graph,columns;
+
+    rows_graph.resize(inputGraph->GetNodeNumVertices()+1);
+    columns.resize(inputGraph->GetNodeNumEdges());
+
+    LO realnnz = 0;
+    GO numDropped = 0, numTotal = 0;
+    for (LO row = 0; row < Teuchos::as<LO>(inputGraph->GetDomainMap()->getNodeNumElements()); ++row) {
+      LO row_block = row_block_number[row];
+      ArrayView<const LO> indices = inputGraph->getNeighborVertices(row);
+
+      LO rownnz = 0;
+      for (LO colID = 0; colID < Teuchos::as<LO>(indices.size()); colID++) {
+        LO col = indices[colID];
+        LO col_block = col_block_number[col];
+
+        if(row_block == col_block) {
+          columns[realnnz++] = col;
+          rownnz++;
+        } else
+          numDropped++;
+      }
+      rows_graph[row+1] = realnnz;
+    }
+
+
+    columns.resize(realnnz);
+    numTotal = inputGraph->GetNodeNumEdges();
+
+    if (GetVerbLevel() & Statistics1) {
+      RCP<const Teuchos::Comm<int> > comm = inputGraph->GetDomainMap()->getComm();
+      GO numGlobalTotal, numGlobalDropped;
+      MueLu_sumAll(comm, numTotal,   numGlobalTotal);
+      MueLu_sumAll(comm, numDropped, numGlobalDropped);
+      GetOStream(Statistics1) << "Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped << "/" << numGlobalTotal;
+      if (numGlobalTotal != 0)
+        GetOStream(Statistics1) << " (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) << "%)";
+      GetOStream(Statistics1) << std::endl;
+    }
+
+    outputGraph =  rcp(new LWGraph(rows_graph, columns, inputGraph->GetDomainMap(), inputGraph->GetImportMap(), "block-diagonalized graph of A"));
+    outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
+
+
+
+ }
+
+
+
 } //namespace MueLu
 
 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP

packages/muelu/src/Interface/MueLu_HierarchyManager.hpp

@@ -174,7 +174,7 @@ namespace MueLu {
       H.SetMaxCoarseSize(maxCoarseSize_);
       VerboseObject::SetDefaultVerbLevel(verbosity_);
       if (graphOutputLevel_ >= 0)
-        H.EnableGraphDumping("dep_graph.dot", graphOutputLevel_);
+        H.EnableGraphDumping("dep_graph", graphOutputLevel_);
 
       if (VerboseObject::IsPrint(Statistics2)) {
         RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(Op);