MueLu: add replace tol and value options

This commit adds two interface options. The first is to control
the tolerance for deciding when to replace a diagonal value
when using the rowsumabs option, the second is to control the
replacement value.

packages/muelu/doc/UsersGuide/masterList.xml

@@ -1105,6 +1105,24 @@ Only used when tentative: calculate qr is set to false.</description>
       <name-ML>not supported by ML</name-ML>
     </parameter>
 
+    <parameter>
+      <name>sa: rowsumabs diagonal replacement tolerance</name>
+      <type>double</type>
+      <default>-1.0</default>
+      <description>If not -1, use this value for deciding whether a diagonal value in prolongator smoothing is too small and should be replaced.</description>
+      <visible>false</visible>
+      <name-ML>not supported by ML</name-ML>
+    </parameter>
+
+    <parameter>
+      <name>sa: rowsumabs diagonal replacement value</name>
+      <type>double</type>
+      <default>0.0</default>
+      <description>If it's determined that a diagonal entry in prolongator smoothing is too small, replace that entry with this value.</description>
+      <visible>false</visible>
+      <name-ML>not supported by ML</name-ML>
+    </parameter>
+
     <parameter>
       <name>interp: interpolation order</name>
       <type>int</type>

packages/muelu/doc/UsersGuide/paramlist_hidden.tex

@@ -241,6 +241,10 @@
 
 \cbb{sa: max eigenvalue}{double}{-1.0}{If not -1, use this value for lambdaMax(Dinv A) in prolongator smoothing instead of calculating it.}
 
+\cbb{sa: rowsumabs diagonal replacement tolerance}{double}{-1.0}{If not -1, use this value for deciding whether a diagonal value in prolongator smoothing is too small and should be replaced.}
+
+\cbb{sa: rowsumabs diagonal replacement value}{double}{0.0}{If it's determined that a diagonal entry in prolongator smoothing is too small, replace that entry with this value.}
+
 \cbb{interp: interpolation order}{int}{1}{Interpolation order used to interpolate values from coarse points to fine points. Possible values are 0 for piece-wise constant interpolation and 1 for piece-wise linear interpolation. This parameter is set to 1 by default.}
 
 \cbb{interp: build coarse coordinates}{bool}{true}{If false, skip the calculation of coarse coordinates.}

packages/muelu/src/Interface/MueLu_ParameterListInterpreter_def.hpp

@@ -664,6 +664,9 @@ namespace MueLu {
     bool useMaxAbsDiagonalScaling = false;
     if (defaultList.isParameter("sa: use rowsumabs diagonal scaling"))
       useMaxAbsDiagonalScaling = defaultList.get<bool>("sa: use rowsumabs diagonal scaling");
+    double chebyReplaceTol = -Teuchos::ScalarTraits<double>::one();
+    if (defaultList.isParameter("sa: rowsumabs diagonal replacement tolerance"))
+      chebyReplaceTol = defaultList.get<double>("sa: rowsumabs diagonal replacement tolerance");
 
     // === Smoothing ===
     // FIXME: should custom smoother check default list too?
@@ -735,6 +738,10 @@ namespace MueLu {
 
         if (preSmootherType == "CHEBYSHEV" && useMaxAbsDiagonalScaling)
           preSmootherParams.set("chebyshev: use rowsumabs diagonal scaling",true);
+        if (preSmootherType == "CHEBYSHEV" && chebyReplaceTol != -Teuchos::ScalarTraits<double>::one())
+          preSmootherParams.set("chebyshev: rowsumabs diagonal replacement tolerance",chebyReplaceTol);
+        if (preSmootherType == "CHEBYSHEV" && defaultList.isParameter("sa: rowsumabs diagonal replacement value"))
+          preSmootherParams.set("chebyshev: rowsumabs diagonal replacement value", defaultList.get<double>("sa: rowsumabs diagonal replacement value"));
 
 #ifdef HAVE_MUELU_INTREPID2
       // Propagate P-coarsening for Topo smoothing
@@ -1694,6 +1701,8 @@ namespace MueLu {
     MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: max eigenvalue", double, Pparams);
     MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: eigenvalue estimate num iterations", int, Pparams);
     MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: use rowsumabs diagonal scaling", bool, Pparams);
+    MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement tolerance", double, Pparams);
+    MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: rowsumabs diagonal replacement value", double, Pparams);
     MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "sa: enforce constraints", bool, Pparams);
     MUELU_TEST_AND_SET_PARAM_2LIST(paramList, defaultList, "tentative: calculate qr", bool, Pparams);
 

packages/muelu/src/MueCentral/MueLu_MasterList.cpp

@@ -136,6 +136,8 @@ namespace MueLu {
     if (name == "sa: use rowsumabs diagonal scaling") { ss << "<Parameter name=\"sa: use rowsumabs diagonal scaling\" type=\"bool\" value=" << value << "/>"; return ss.str(); }      
     if (name == "sa: enforce constraints") { ss << "<Parameter name=\"sa: enforce constraints\" type=\"bool\" value=" << value << "/>"; return ss.str(); }      
     if (name == "sa: max eigenvalue") { ss << "<Parameter name=\"sa: max eigenvalue\" type=\"double\" value=" << value << "/>"; return ss.str(); }      
+    if (name == "sa: rowsumabs diagonal replacement tolerance") { ss << "<Parameter name=\"sa: rowsumabs diagonal replacement tolerance\" type=\"double\" value=" << value << "/>"; return ss.str(); }      
+    if (name == "sa: rowsumabs diagonal replacement value") { ss << "<Parameter name=\"sa: rowsumabs diagonal replacement value\" type=\"double\" value=" << value << "/>"; return ss.str(); }      
     if (name == "pcoarsen: element") { ss << "<Parameter name=\"pcoarsen: element\" type=\"string\" value=" << value << "/>"; return ss.str(); }      
     if (name == "pcoarsen: schedule") { ss << "<Parameter name=\"pcoarsen: schedule\" type=\"string\" value=" << value << "/>"; return ss.str(); }      
     if (name == "pcoarsen: hi basis") { ss << "<Parameter name=\"pcoarsen: hi basis\" type=\"string\" value=" << value << "/>"; return ss.str(); }      
@@ -275,6 +277,8 @@ namespace MueLu {
   "<Parameter name=\"sa: use rowsumabs diagonal scaling\" type=\"bool\" value=\"false\"/>"
   "<Parameter name=\"sa: enforce constraints\" type=\"bool\" value=\"false\"/>"
   "<Parameter name=\"sa: max eigenvalue\" type=\"double\" value=\"-1.0\"/>"
+  "<Parameter name=\"sa: rowsumabs diagonal replacement tolerance\" type=\"double\" value=\"-1.0\"/>"
+  "<Parameter name=\"sa: rowsumabs diagonal replacement value\" type=\"double\" value=\"0.0\"/>"
   "<Parameter name=\"interp: interpolation order\" type=\"int\" value=\"1\"/>"
   "<Parameter name=\"interp: build coarse coordinates\" type=\"bool\" value=\"true\"/>"
   "<ParameterList name=\"transfer: params\"/>"
@@ -746,6 +750,10 @@ namespace MueLu {
 
          ("not supported by ML","sa: max eigenvalue")
 
+         ("not supported by ML","sa: rowsumabs diagonal replacement tolerance")
+
+         ("not supported by ML","sa: rowsumabs diagonal replacement value")
+
          ("interp: interpolation order","interp: interpolation order")
 
          ("interp: build coarse coordinates","interp: build coarse coordinates")

packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp

@@ -643,8 +643,18 @@ namespace MueLu {
         bool doScale = false;
         doScale = paramList.get<bool>("chebyshev: use rowsumabs diagonal scaling");
         paramList.remove("chebyshev: use rowsumabs diagonal scaling");
+        double chebyReplaceTol = Teuchos::ScalarTraits<Scalar>::eps()*100;
+        if (paramList.isParameter("chebyshev: rowsumabs diagonal replacement tolerance")) {
+          paramList.get<double>("chebyshev: rowsumabs diagonal replacement tolerance",chebyReplaceTol);
+          paramList.remove("chebyshev: rowsumabs diagonal replacement tolerance");
+        }
+        double chebyReplaceVal = Teuchos::ScalarTraits<double>::zero();
+        if (paramList.isParameter("chebyshev: rowsumabs diagonal replacement value")) {
+          paramList.get<double>("chebyshev: rowsumabs diagonal replacement value",chebyReplaceVal);
+          paramList.remove("chebyshev: rowsumabs diagonal replacement value");
+        }
         if (doScale) {
-          RCP<Vector> lumpedDiagonal = Utilities::GetLumpedMatrixDiagonal(*(currentLevel.Get<RCP<Matrix> >("A")),true);
+          RCP<Vector> lumpedDiagonal = Utilities::GetLumpedMatrixDiagonal(*(currentLevel.Get<RCP<Matrix> >("A")),true, chebyReplaceTol, chebyReplaceVal);
           const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*lumpedDiagonal);
           paramList.set("chebyshev: operator inv diagonal",tmpVec.getTpetra_Vector());
         }

packages/muelu/src/Transfers/Smoothed-Aggregation/MueLu_SaPFactory_def.hpp

@@ -76,6 +76,8 @@ namespace MueLu {
     SET_VALID_ENTRY("sa: enforce constraints");
     SET_VALID_ENTRY("tentative: calculate qr");
     SET_VALID_ENTRY("sa: max eigenvalue");
+    SET_VALID_ENTRY("sa: rowsumabs diagonal replacement tolerance");
+    SET_VALID_ENTRY("sa: rowsumabs diagonal replacement value");
 #undef  SET_VALID_ENTRY
 
     validParamList->set< RCP<const FactoryBase> >("A",              Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
@@ -159,6 +161,10 @@ namespace MueLu {
     const bool doQRStep         =        pL.get<bool>("tentative: calculate qr");
     const bool enforceConstraints=       pL.get<bool>("sa: enforce constraints");
     const SC   userDefinedMaxEigen =  as<SC>(pL.get<double>("sa: max eigenvalue"));
+    typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
+    double dTol = pL.get<double>("sa: rowsumabs diagonal replacement tolerance");
+    const Magnitude diagonalReplacementTolerance = (dTol == as<double>(-1) ? Teuchos::ScalarTraits<Scalar>::eps()*100 : as<Magnitude>(pL.get<double>("sa: rowsumabs diagonal replacement tolerance")));
+    const SC diagonalReplacementValue =  as<SC>(pL.get<double>("sa: rowsumabs diagonal replacement value"));
 
     // Sanity checking
     TEUCHOS_TEST_FOR_EXCEPTION(doQRStep && enforceConstraints,Exceptions::RuntimeError,
@@ -180,8 +186,8 @@ namespace MueLu {
             const bool returnReciprocal=true;
             const bool replaceSingleEntryRowWithZero=true;
             invDiag = Utilities::GetLumpedMatrixDiagonal(*A,returnReciprocal,
-                                                        Teuchos::ScalarTraits<Scalar>::eps()*100,
-                                                        Teuchos::ScalarTraits<Scalar>::zero(),
+                                                        diagonalReplacementTolerance,
+                                                        diagonalReplacementValue,
                                                         replaceSingleEntryRowWithZero);
             TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(), Exceptions::RuntimeError,
                                        "SaPFactory: eigenvalue estimate: diagonal reciprocal is null.");
@@ -208,8 +214,8 @@ namespace MueLu {
           const bool returnReciprocal=true;
           const bool replaceSingleEntryRowWithZero=true;
           invDiag = Utilities::GetLumpedMatrixDiagonal(*A,returnReciprocal,
-                                                        Teuchos::ScalarTraits<Scalar>::eps()*100,
-                                                        Teuchos::ScalarTraits<Scalar>::zero(),
+                                                        diagonalReplacementTolerance,
+                                                        diagonalReplacementValue,
                                                         replaceSingleEntryRowWithZero);
           TEUCHOS_TEST_FOR_EXCEPTION(invDiag.is_null(), Exceptions::RuntimeError, "SaPFactory: diagonal reciprocal is null.");
         }	

packages/muelu/src/Utils/MueLu_UtilitiesBase_decl.hpp

@@ -163,15 +163,15 @@ namespace MueLu {
 
     NOTE -- it's assumed that A has been fillComplete'd.
     */
-    static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
+    static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
       Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase::GetMatrixDiagonalInverse");
 
       RCP<const Map> rowMap = A.getRowMap();
       RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
 
       A.getLocalDiagCopy(*diag);
 
-      RCP<Vector> inv = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(diag, tol, tolReplacement);
+      RCP<Vector> inv = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(diag, tol, valReplacement);
 
       return inv;
     }
@@ -184,7 +184,7 @@ namespace MueLu {
     */
     static Teuchos::RCP<Vector> GetLumpedMatrixDiagonal(Matrix const & A, const bool doReciprocal = false,
                                                         Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
-                                                        Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
+                                                        Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
                                                         const bool replaceSingleEntryRowWithZero = false) {
 
       typedef Teuchos::ScalarTraits<Scalar> TST;
@@ -194,6 +194,8 @@ namespace MueLu {
       const Scalar one = TST::one();
       const Scalar two = one + one;
 
+tol = 0.;
+
       Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
 
       RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
@@ -231,7 +233,7 @@ namespace MueLu {
               if(TST::magnitude(diagVals[i]) > tol)
                 diagVals[i] = one / diagVals[i];
               else {
-                diagVals[i] = tolReplacement;
+                diagVals[i] = valReplacement;
               }
             }
           }
@@ -262,11 +264,11 @@ namespace MueLu {
     /*! @brief Return vector containing inverse of input vector
      *
      * @param[in] v: input vector
-     * @param[in] tol: tolerance. If entries of input vector are smaller than tolerance they are replaced by tolReplacement (see below). The default value for tol is 100*eps (machine precision)
-     * @param[in] tolReplacement: Value put in for undefined entries in output vector (default: 0.0)
+     * @param[in] tol: tolerance. If entries of input vector are smaller than tolerance they are replaced by valReplacement (see below). The default value for tol is 100*eps (machine precision)
+     * @param[in] valReplacement: Value put in for undefined entries in output vector (default: 0.0)
      * @ret: vector containing inverse values of input vector v
     */
-    static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
+    static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
 
       RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),true);
 
@@ -279,7 +281,7 @@ namespace MueLu {
         for(size_t r = 0; r < bmap->getNumMaps(); ++r) {
           RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
           RCP<const Vector> subvec = submvec->getVector(0);
-          RCP<Vector> subvecinf = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(subvec,tol,tolReplacement);
+          RCP<Vector> subvecinf = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(subvec,tol,valReplacement);
           bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
         }
         return ret;
@@ -292,7 +294,7 @@ namespace MueLu {
         if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
           retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
         else
-          retVals[i] = tolReplacement;
+          retVals[i] = valReplacement;
       }
       return ret;
     }

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

@@ -86,6 +86,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 

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

@@ -83,6 +83,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -188,6 +190,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -293,6 +297,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 

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

@@ -83,6 +83,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -188,6 +190,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -293,6 +297,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -398,6 +404,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]
 
@@ -503,6 +511,8 @@ sa: use rowsumabs diagonal scaling = 0   [default]
 sa: enforce constraints = 0   [default]
 tentative: calculate qr = 1   [default]
 sa: max eigenvalue = -1   [default]
+sa: rowsumabs diagonal replacement tolerance = -1   [default]
+sa: rowsumabs diagonal replacement value = 0   [default]
 matrixmatrix: kernel params -> 
  [empty list]