Merge Pull Request #12025 from trilinos/Trilinos/master_merge_2023070…

…7_175815

Automatically Merged using Trilinos Master Merge AutoTester
PR Title: Trilinos Master Merge PR Generator: Auto PR created to promote from master_merge_20230707_175815 branch to master
PR Author: trilinos-autotester

packages/kokkos-kernels/ode/CMakeLists.txt

@@ -11,5 +11,8 @@ ENDIF()
 
 
 # Adding unit-tests
-KOKKOSKERNELS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}/ode)
+# Note BMK: Since ODE has no auto-generated ETI files, this directory does not exist in a build without unit tests.
+# This causes configure errors when building an app against a Trilinos install, and the unit test build dir doesn't contain any headers that need to be found.
+# But uncomment the next line if ETI headers are ever added.
+# KOKKOSKERNELS_INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR}/ode)
 KOKKOSKERNELS_INCLUDE_DIRECTORIES(REQUIRED_DURING_INSTALLATION_TESTING ${CMAKE_CURRENT_SOURCE_DIR}/ode)

packages/ml/src/Utils/ml_MultiLevelPreconditioner.cpp

@@ -2242,8 +2242,10 @@ ComputePreconditioner(const bool CheckPreconditioner)
   if(AMGSolver_ != ML_CLASSICAL_FAMILY) 
     ML_CHK_ERR(SetupCoordinates());
 
-  if (List_.get("RAP: sort columns",0))                                     //
+  if (List_.get("RAP: sort columns",0)) {                                    
+    if (ml_nodes_) ml_nodes_->sortColumnsAfterRAP = 1;
     ml_->sortColumnsAfterRAP = 1;
+  }
   // ========================================================================//
   //               Setting Repartitioning                                    //
   // ========================================================================//

packages/muelu/doc/UsersGuide/masterList.xml

@@ -521,6 +521,16 @@
       <name-ML>aggregation: threshold</name-ML>
     </parameter>
 
+  <parameter>
+      <name>aggregation: use ml scaling of drop tol</name>
+      <type>bool</type>
+      <default>false</default>
+      <description>Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.</description>
+      <comment-ML>parameter not existing in ML</comment-ML>
+    </parameter>
+
+
+
     <parameter>
       <name>aggregation: min agg size</name>
       <type>int</type>

packages/muelu/doc/UsersGuide/options_aggregation.tex

@@ -13,6 +13,8 @@
 
 \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.}
 
+\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.}
+
 \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.}
 
 \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).}

packages/muelu/doc/UsersGuide/paramlist.tex

@@ -66,6 +66,8 @@
 
 \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.}
 
+\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.}
+
 \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.}
 
 \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).}

packages/muelu/doc/UsersGuide/paramlist_hidden.tex

@@ -95,6 +95,8 @@
 
 \cbb{aggregation: drop tol}{double}{0.0}{Connectivity dropping threshold for a graph used in aggregation.}
 
+\cbb{aggregation: use ml scaling of drop tol}{bool}{false}{Enables ML-style scaling of drop tol, where the drop tol halves with each successive level.}
+
 \cbb{aggregation: min agg size}{int}{2}{Minimum size of an aggregate.}
 
 \cbb{aggregation: max agg size}{int}{-1}{Maximum size of an aggregate (-1 means unlimited).}

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

@@ -122,6 +122,7 @@ namespace MueLu {
 
 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
     SET_VALID_ENTRY("aggregation: drop tol");
+    SET_VALID_ENTRY("aggregation: use ml scaling of drop tol");
     SET_VALID_ENTRY("aggregation: Dirichlet threshold");
     SET_VALID_ENTRY("aggregation: greedy Dirichlet");
     SET_VALID_ENTRY("aggregation: row sum drop tol");
@@ -205,8 +206,10 @@ namespace MueLu {
     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
+    // in the block diagonalization). 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;
 
@@ -326,7 +329,14 @@ namespace MueLu {
       TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null   && algo != "classical", Exceptions::RuntimeError, "Dropping function must not be provided for \"" << algo << "\" algorithm");
       TEUCHOS_TEST_FOR_EXCEPTION(algo != "classical" && algo != "distance laplacian" && algo != "signed classical", Exceptions::RuntimeError, "\"algorithm\" must be one of (classical|distance laplacian|signed classical)");
 
-      SC threshold = as<SC>(pL.get<double>("aggregation: drop tol"));
+      SC threshold;
+    // If we're doing the ML-style halving of the drop tol at each level, we do that here.
+    if (pL.get<bool>("aggregation: use ml scaling of drop tol"))
+      threshold = pL.get<double>("aggregation: drop tol") / pow(2.0,currentLevel.GetLevelID());
+    else
+      threshold =  as<SC>(pL.get<double>("aggregation: drop tol"));
+
+
       std::string distanceLaplacianAlgoStr = pL.get<std::string>("aggregation: distance laplacian algo");
       std::string classicalAlgoStr = pL.get<std::string>("aggregation: classical algo");
       real_type realThreshold = STS::magnitude(threshold);// CMS: Rename this to "magnitude threshold" sometime

packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp

@@ -399,6 +399,9 @@ namespace MueLu {
       // make sure that MueLu's phase2a matches ML's
       mueluss << "<Parameter name=\"aggregation: match ML phase2a\"      type=\"bool\"     value=\"true\"/>" << std::endl;
 
+      // make sure that MueLu's drop tol matches ML's
+      mueluss << "<Parameter name=\"aggregation: use ml scaling of drop tol\"      type=\"bool\"     value=\"true\"/>" << std::endl;
+
       // special handling for energy minimization
       // TAW: this is not optimal for symmetric problems but at least works.
       //      for symmetric problems the "energy minimization" parameter should not exist anyway...

packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp

@@ -1090,6 +1090,8 @@ namespace MueLu {
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: row sum drop tol",        double, dropParams);
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: block diagonal: interleaved blocksize", int, dropParams);
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: drop tol",                     double, dropParams);
+       MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: use ml scaling of drop tol",   bool, dropParams);
+
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: Dirichlet threshold",          double, dropParams);
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: greedy Dirichlet",          bool, dropParams);
        MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "aggregation: distance laplacian algo", std::string, dropParams);

packages/muelu/src/MueCentral/MueLu_MasterList.cpp

@@ -221,6 +221,7 @@ namespace MueLu {
   "<Parameter name=\"aggregation: distance laplacian algo\" type=\"string\" value=\"default\"/>"
   "<Parameter name=\"aggregation: classical algo\" type=\"string\" value=\"default\"/>"
   "<Parameter name=\"aggregation: drop tol\" type=\"double\" value=\"0.0\"/>"
+  "<Parameter name=\"aggregation: use ml scaling of drop tol\" type=\"bool\" value=\"false\"/>"
   "<Parameter name=\"aggregation: min agg size\" type=\"int\" value=\"2\"/>"
   "<Parameter name=\"aggregation: max agg size\" type=\"int\" value=\"-1\"/>"
   "<Parameter name=\"aggregation: compute aggregate qualities\" type=\"bool\" value=\"false\"/>"
@@ -661,6 +662,8 @@ namespace MueLu {
 
          ("aggregation: threshold","aggregation: drop tol")
 
+         ("aggregation: use ml scaling of drop tol","aggregation: use ml scaling of drop tol")
+
          ("aggregation: min agg size","aggregation: min agg size")
 
          ("aggregation: max agg size","aggregation: max agg size")

packages/muelu/src/Operators/MueLu_Maxwell1_def.hpp

@@ -112,6 +112,10 @@ namespace MueLu {
       newList.sublist("maxwell1: 11list").set("tentative: constant column sums", false);
       newList.sublist("maxwell1: 11list").set("tentative: calculate qr", false);
 
+      newList.sublist("maxwell1: 11list").set("aggregation: use ml scaling of drop tol",true);
+      newList.sublist("maxwell1: 22list").set("aggregation: use ml scaling of drop tol",true);
+
+
       if(list.isParameter("aggregation: damping factor") && list.get<double>("aggregation: damping factor") == 0.0)
         newList.sublist("maxwell1: 11list").set("multigrid algorithm", "unsmoothed reitzinger");
       else

packages/muelu/src/Smoothers/BlockedSmoothers/MueLu_BlockedGaussSeidelSmoother_def.hpp

@@ -201,85 +201,73 @@ namespace MueLu {
 
     // make sure that both rcpX and rcpB are BlockedMultiVector objects
     bool bCopyResultX = false;
+    bool bReorderX = false;
+    bool bReorderB = false;
     RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
     MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
     RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
     RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
 
-    if(bX.is_null() == true) {
-      RCP<MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
-      rcpX.swap(test);
-      bCopyResultX = true;
+    // check the type of operator
+    RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA =
+            Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
+
+    if(rbA.is_null() == false) {
+      // A is a ReorderedBlockedCrsMatrix and thus uses nested maps, retrieve BlockedCrsMatrix and use plain blocked
+      // map for the construction of the blocked multivectors
+
+      // check type of X vector
+      if (bX.is_null() == true) {
+        RCP<MultiVector> vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedDomainMap(), rcpX));
+        rcpX.swap(vectorWithBlockedMap);
+        bCopyResultX = true;
+        bReorderX = true;
+      }
+
+      // check type of B vector
+      if (bB.is_null() == true) {
+        RCP<const MultiVector> vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(rbA->getBlockedCrsMatrix()->getBlockedRangeMap(), rcpB));
+        rcpB.swap(vectorWithBlockedMap);
+        bReorderB = true;
+      }
     }
+    else {
+      // A is a BlockedCrsMatrix and uses a plain blocked map
+      if (bX.is_null() == true) {
+        RCP<MultiVector> vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(), rcpX));
+        rcpX.swap(vectorWithBlockedMap);
+        bCopyResultX = true;
+      }
 
-    if(bB.is_null() == true) {
-      RCP<const MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
-      rcpB.swap(test);
+      if(bB.is_null() == true) {
+        RCP<const MultiVector> vectorWithBlockedMap = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
+        rcpB.swap(vectorWithBlockedMap);
+      }
     }
 
     // we now can guarantee that X and B are blocked multi vectors
     bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
     bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
 
-    // check the type of operator
-    RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
+    // Finally we can do a reordering of the blocked multivectors if A is a ReorderedBlockedCrsMatrix
     if(rbA.is_null() == false) {
-      // A is a ReorderedBlockedCrsMatrix
+
       Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
 
       // check type of X vector
-      if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
+      if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) {
         // X is a blocked multi vector but incompatible to the reordered blocked operator A
-        Teuchos::RCP<MultiVector> test =
-            buildReorderedBlockedMultiVector(brm, bX);
-        rcpX.swap(test);
+        Teuchos::RCP<MultiVector> reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bX);
+        rcpX.swap(reorderedBlockedVector);
       }
-      if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
+
+      if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) {
         // B is a blocked multi vector but incompatible to the reordered blocked operator A
-        Teuchos::RCP<const MultiVector> test =
-            buildReorderedBlockedMultiVector(brm, bB);
-        rcpB.swap(test);
+        Teuchos::RCP<const MultiVector> reorderedBlockedVector = buildReorderedBlockedMultiVector(brm, bB);
+        rcpB.swap(reorderedBlockedVector);
       }
     }
 
-
-    // check the type of operator
-    /*RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
-    MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
-    RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
-    if(rbA.is_null() == false) {
-      // A is a ReorderedBlockedCrsMatrix
-      Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
-
-      // check type of vectors
-      RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
-      RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
-
-      // check type of X vector
-      if(bX.is_null() == false && bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
-        RCP<ReorderedBlockedMultiVector> rbX = Teuchos::rcp_dynamic_cast<ReorderedBlockedMultiVector>(bX);
-        if(rbX.is_null() == true) {
-          // X is a blocked multi vector but not reordered
-          // However, A is a reordered blocked operator
-          // We have to make sure, that A and X use compatible maps
-          Teuchos::RCP<MultiVector> test =
-              buildReorderedBlockedMultiVector(brm, bX);
-          rcpX.swap(test);
-        }
-      }
-      if(bB.is_null() == false && bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
-        RCP<const ReorderedBlockedMultiVector> rbB = Teuchos::rcp_dynamic_cast<const ReorderedBlockedMultiVector>(bB);
-        if(rbB.is_null() == true) {
-          // B is a blocked multi vector but not reordered
-          // However, A is a reordered blocked operator
-          // We have to make sure, that A and X use compatible maps
-          Teuchos::RCP<const MultiVector> test =
-              buildReorderedBlockedMultiVector(brm, bB);
-          rcpB.swap(test);
-        }
-      }
-    }*/
-
     // Throughout the rest of the algorithm rcpX and rcpB are used for solution vector and RHS
 
     RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());