Merge pull request trilinos#9095 from cgcgcg/ifpack2ChebyReuse

Ifpack2 cheby reuse

packages/ifpack2/src/Ifpack2_AdditiveSchwarz_decl.hpp

@@ -645,6 +645,10 @@ class AdditiveSchwarz :
   void
   setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist);
 
+  bool supportsZeroStartingSolution() { return true; }
+
+  void setZeroStartingSolution (bool zeroStartingSolution) { ZeroStartingSolution_ = zeroStartingSolution; };
+
   /// \brief Get a list of the preconditioner's default parameters.
   ///
   /// See the documentation of setParameters() for a list of the

packages/ifpack2/src/Ifpack2_BlockRelaxation_decl.hpp

@@ -209,6 +209,10 @@ class BlockRelaxation :
   */
   void setParameters(const Teuchos::ParameterList& params);
 
+  bool supportsZeroStartingSolution() { return true; }
+
+  void setZeroStartingSolution (bool zeroStartingSolution) { ZeroStartingSolution_ = zeroStartingSolution; };
+
   //! Return a list of all the parameters that this class accepts.
   Teuchos::RCP<const Teuchos::ParameterList>
   getValidParameters () const;

packages/ifpack2/src/Ifpack2_Chebyshev_decl.hpp

@@ -425,6 +425,10 @@ class Chebyshev :
   /// always provide a row Map or range Map vector for this parameter.
   void setParameters (const Teuchos::ParameterList& params);
 
+  bool supportsZeroStartingSolution() { return true; }
+
+  void setZeroStartingSolution (bool zeroStartingSolution);
+
   /// \brief Initialize the preconditioner.
   ///
   /// The compute() method will call initialize() automatically if it

packages/ifpack2/src/Ifpack2_Chebyshev_def.hpp

@@ -97,6 +97,13 @@ Chebyshev<MatrixType>::setParameters (const Teuchos::ParameterList& List)
 }
 
 
+template<class MatrixType>
+void
+Chebyshev<MatrixType>::setZeroStartingSolution (bool zeroStartingSolution)
+{
+  impl_.setZeroStartingSolution(zeroStartingSolution);
+}
+
 template<class MatrixType>
 Teuchos::RCP<const Teuchos::Comm<int> >
 Chebyshev<MatrixType>::getComm () const

packages/ifpack2/src/Ifpack2_Details_Chebyshev_decl.hpp

@@ -262,6 +262,8 @@ class Chebyshev : public Teuchos::Describable {
   /// from Ifpack.
   void setParameters (Teuchos::ParameterList& plist);
 
+  void setZeroStartingSolution (bool zeroStartingSolution) { zeroStartingSolution_ = zeroStartingSolution; }
+
   /// \brief (Re)compute the left scaling D_inv, and estimate min and
   ///   max eigenvalues of D_inv * A.
   ///
@@ -474,6 +476,16 @@ class Chebyshev : public Teuchos::Describable {
   //! Number of power method iterations for estimating the max eigenvalue.
   int eigMaxIters_;
 
+  //! Relative tolerance for power method iterations for estimating the max eigenvalue.
+  MT eigRelTolerance_;
+
+  //! Whether the iteration vectors of the power method should be saved.
+  bool eigKeepVectors_;
+
+  //! Iteration vectors of the power method. Can be saved for the purpose of multiple setups.
+  Teuchos::RCP<V> eigVector_;
+  Teuchos::RCP<V> eigVector2_;
+
   //! Frequency of normalization in the power method.
   int eigNormalizationFreq_;
 

packages/ifpack2/src/Ifpack2_Details_Chebyshev_def.hpp

@@ -286,6 +286,8 @@ Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
   minDiagVal_ (STS::eps ()),
   numIters_ (1),
   eigMaxIters_ (10),
+  eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
+  eigKeepVectors_(false),
   eigNormalizationFreq_(1),
   zeroStartingSolution_ (true),
   assumeMatrixUnchanged_ (false),
@@ -314,6 +316,8 @@ Chebyshev (Teuchos::RCP<const row_matrix_type> A,
   minDiagVal_ (STS::eps ()),
   numIters_ (1),
   eigMaxIters_ (10),
+  eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
+  eigKeepVectors_(false),
   eigNormalizationFreq_(1),
   zeroStartingSolution_ (true),
   assumeMatrixUnchanged_ (false),
@@ -357,6 +361,8 @@ setParameters (Teuchos::ParameterList& plist)
   const ST defaultMinDiagVal = STS::eps ();
   const int defaultNumIters = 1;
   const int defaultEigMaxIters = 10;
+  const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero ();
+  const bool defaultEigKeepVectors = false;
   const int defaultEigNormalizationFreq = 1;
   const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
   const bool defaultAssumeMatrixUnchanged = false;
@@ -376,6 +382,8 @@ setParameters (Teuchos::ParameterList& plist)
   ST minDiagVal = defaultMinDiagVal;
   int numIters = defaultNumIters;
   int eigMaxIters = defaultEigMaxIters;
+  MT eigRelTolerance = defaultEigRelTolerance;
+  bool eigKeepVectors = defaultEigKeepVectors;
   int eigNormalizationFreq = defaultEigNormalizationFreq;
   bool zeroStartingSolution = defaultZeroStartingSolution;
   bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
@@ -581,6 +589,15 @@ setParameters (Teuchos::ParameterList& plist)
     "\" parameter (also called \"eigen-analysis: iterations\") must be a "
     "nonnegative integer.  You gave a value of " << eigMaxIters << ".");
 
+  if (plist.isType<double>("chebyshev: eigenvalue relative tolerance"))
+    eigRelTolerance = Teuchos::as<MT>(plist.get<double> ("chebyshev: eigenvalue relative tolerance"));
+  else if (plist.isType<MT>("chebyshev: eigenvalue relative tolerance"))
+    eigRelTolerance = plist.get<MT> ("chebyshev: eigenvalue relative tolerance");
+  else if (plist.isType<ST>("chebyshev: eigenvalue relative tolerance"))
+    eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<ST> ("chebyshev: eigenvalue relative tolerance"));
+
+  eigKeepVectors = plist.get ("chebyshev: eigenvalue keep vectors", eigKeepVectors);
+
   eigNormalizationFreq = plist.get ("chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
   TEUCHOS_TEST_FOR_EXCEPTION(
     eigNormalizationFreq < 0, std::invalid_argument,
@@ -650,6 +667,8 @@ setParameters (Teuchos::ParameterList& plist)
   minDiagVal_ = minDiagVal;
   numIters_ = numIters;
   eigMaxIters_ = eigMaxIters;
+  eigRelTolerance_ = eigRelTolerance;
+  eigKeepVectors_ = eigKeepVectors;
   eigNormalizationFreq_ = eigNormalizationFreq;
   zeroStartingSolution_ = zeroStartingSolution;
   assumeMatrixUnchanged_ = assumeMatrixUnchanged;
@@ -696,6 +715,8 @@ Chebyshev<ScalarType, MV>::reset ()
   W_ = Teuchos::null;
   computedLambdaMax_ = STS::nan ();
   computedLambdaMin_ = STS::nan ();
+  eigVector_ = Teuchos::null;
+  eigVector2_ = Teuchos::null;
 }
 
 
@@ -1359,9 +1380,16 @@ powerMethodWithInitGuess (const op_type& A,
   const ST zero = static_cast<ST> (0.0);
   const ST one = static_cast<ST> (1.0);
   ST lambdaMax = zero;
-  ST RQ_top, RQ_bottom = one, norm;
-
-  V y (A.getRangeMap ());
+  ST lambdaMaxOld = one;
+  ST norm;
+
+  Teuchos::RCP<V> y;
+  if (eigVector2_.is_null()) {
+    y = rcp(new V(A.getRangeMap ()));
+    if (eigKeepVectors_)
+      eigVector2_ = y;
+  } else
+    y = eigVector2_;
   norm = x.norm2 ();
   TEUCHOS_TEST_FOR_EXCEPTION
     (norm == zero, std::runtime_error,
@@ -1386,16 +1414,31 @@ powerMethodWithInitGuess (const op_type& A,
     if (debug_) {
       *out_ << "  Iteration " << iter << endl;
     }
-    A.apply (x, y);
-    solve (x, D_inv, y);
+    A.apply (x, *y);
+    solve (x, D_inv, *y);
 
     if (((iter+1) % eigNormalizationFreq_ == 0) && (iter < numIters-2)) {
       norm = x.norm2 ();
       if (norm == zero) { // Return something reasonable.
         if (debug_) {
           *out_ << "   norm is zero; returning zero" << endl;
+          *out_ << "   Power method terminated after "<< iter << " iterations." << endl;
         }
         return zero;
+      } else {
+        lambdaMaxOld = lambdaMax;
+        lambdaMax = pow(norm, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq_);
+        if (Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld) < eigRelTolerance_ * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
+          if (debug_) {
+            *out_ << "  lambdaMax: " << lambdaMax << endl;
+            *out_ << "  Power method terminated after "<< iter << " iterations." << endl;
+          }
+          return lambdaMax;
+        } else if (debug_) {
+          *out_ << "  lambdaMaxOld: " << lambdaMaxOld << endl;
+          *out_ << "  lambdaMax: " << lambdaMax << endl;
+          *out_ << "  |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld)/Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
+        }
       }
       x.scale (one / norm);
     }
@@ -1408,19 +1451,19 @@ powerMethodWithInitGuess (const op_type& A,
   if (norm == zero) { // Return something reasonable.
     if (debug_) {
       *out_ << "   norm is zero; returning zero" << endl;
+      *out_ << "   Power method terminated after "<< numIters << " iterations." << endl;
     }
     return zero;
   }
   x.scale (one / norm);
-  A.apply (x, y);
-  solve (y, D_inv, y);
-  RQ_top = y.dot (x);
-  RQ_bottom = one;
+  A.apply (x, *y);
+  solve (*y, D_inv, *y);
+  lambdaMax = y->dot (x);
   if (debug_) {
-    *out_ << "   RQ_top: " << RQ_top
-          << ", RQ_bottom: " << RQ_bottom << endl;
+    *out_ << "  lambdaMax: " << lambdaMax << endl;
+    *out_ << "  Power method terminated after "<< numIters << " iterations." << endl;
   }
-  lambdaMax = RQ_top / RQ_bottom;
+
   return lambdaMax;
 }
 
@@ -1456,13 +1499,19 @@ powerMethod (const op_type& A, const V& D_inv, const int numIters)
   }
 
   const ST zero = static_cast<ST> (0.0);
-  V x (A.getDomainMap ());
-  // For the first pass, just let the pseudorandom number generator
-  // fill x with whatever values it wants; don't try to make its
-  // entries nonnegative.
-  computeInitialGuessForPowerMethod (x, false);
-
-  ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
+  Teuchos::RCP<V> x;
+  if (eigVector_.is_null()) {
+    x = rcp(new V(A.getDomainMap ()));
+    if (eigKeepVectors_)
+      eigVector_ = x;
+    // For the first pass, just let the pseudorandom number generator
+    // fill x with whatever values it wants; don't try to make its
+    // entries nonnegative.
+    computeInitialGuessForPowerMethod (*x, false);
+  } else
+    x = eigVector_;
+
+  ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
 
   // mfh 07 Jan 2015: Taking the real part here is only a concession
   // to the compiler, so that this class can build with ScalarType =
@@ -1483,8 +1532,8 @@ powerMethod (const op_type& A, const V& D_inv, const int numIters)
 
     // For the second pass, make all the entries of the initial guess
     // vector have nonnegative real parts.
-    computeInitialGuessForPowerMethod (x, true);
-    lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
+    computeInitialGuessForPowerMethod (*x, true);
+    lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
   }
   return lambdaMax;
 }
@@ -1672,6 +1721,7 @@ describe (Teuchos::FancyOStream& out,
           << "userEigRatio_: " << userEigRatio_ << endl
           << "numIters_: " << numIters_ << endl
           << "eigMaxIters_: " << eigMaxIters_ << endl
+          << "eigRelTolerance_: " << eigRelTolerance_ << endl
           << "eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
           << "zeroStartingSolution_: " << zeroStartingSolution_ << endl
           << "assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl

packages/ifpack2/src/Ifpack2_Hiptmair_decl.hpp

@@ -134,6 +134,10 @@ namespace Ifpack2 {
     ///   - "hiptmair: zero starting solution" (\c bool)
     void setParameters (const Teuchos::ParameterList& params);
 
+    bool supportsZeroStartingSolution() { return true; }
+
+    void setZeroStartingSolution (bool zeroStartingSolution) { ZeroStartingSolution_ = zeroStartingSolution; };
+
     //! Do any initialization that depends on the input matrix's structure.
     void initialize ();
 

packages/ifpack2/src/Ifpack2_Parameters.cpp

@@ -73,6 +73,8 @@ void getValidParameters(Teuchos::ParameterList& params)
   params.set("chebyshev: min eigenvalue", 30.0);
   params.set("chebyshev: degree", 1);
   params.set("chebyshev: eigenvalue max iterations", 10);
+  params.set("chebyshev: eigenvalue relative tolerance", 0.0);
+  params.set("chebyshev: eigenvalue keep vector", false);
   params.set("chebyshev: assume matrix does not change", false);
   // params.set("chebyshev: operator inv diagonal",Teuchos::null);
   params.set("chebyshev: min diagonal value", STS::eps());

packages/ifpack2/src/Ifpack2_Preconditioner.hpp

@@ -146,6 +146,11 @@ class Preconditioner :
   //! Set this preconditioner's parameters.
   virtual void setParameters (const Teuchos::ParameterList& List) = 0;
 
+  virtual bool supportsZeroStartingSolution () { return false; };
+
+  //! Set this preconditioner's parameters.
+  virtual void setZeroStartingSolution (bool zeroStartingSolution) { TEUCHOS_ASSERT(false); };
+
   /// @brief Set up the graph structure of this preconditioner.
   ///
   /// If the graph structure of the constructor's input matrix has

packages/ifpack2/src/Ifpack2_Relaxation_decl.hpp

@@ -403,6 +403,10 @@ class Relaxation :
   /// then magnitude_type is \c T.
   void setParameters (const Teuchos::ParameterList& params);
 
+  bool supportsZeroStartingSolution() { return true; }
+
+  void setZeroStartingSolution (bool zeroStartingSolution) { ZeroStartingSolution_ = zeroStartingSolution; };
+
   //! Return a list of all the parameters that this class accepts.
   Teuchos::RCP<const Teuchos::ParameterList>
   getValidParameters () const;

packages/muelu/src/Smoothers/MueLu_Ifpack2Smoother_def.hpp

@@ -769,53 +769,13 @@ namespace MueLu {
     //        initial value at the end but there is no way right now to get
     //        the current value of the "zero starting solution" in ifpack2.
     //        It's not really an issue, as prec_  can only be used by this method.
-    // TODO: When https://software.sandia.gov/bugzilla/show_bug.cgi?id=5283#c2 is done
-    // we should remove the if/else/elseif and just test if this
-    // option is supported by current ifpack2 preconditioner
     Teuchos::ParameterList paramList;
     bool supportInitialGuess = false;
     const Teuchos::ParameterList params = this->GetParameterList();
-    if (type_ == "CHEBYSHEV") {
-      const std::string paramName = "chebyshev: zero starting solution";
-      if (!params.isType<bool>(paramName) ||
-          (params.get<bool>(paramName) != InitialGuessIsZero)) {
-        paramList.set(paramName, InitialGuessIsZero);
-        SetPrecParameters(paramList);
-      }
-      supportInitialGuess = true;
 
-    } else if (type_ == "RELAXATION"       ||
-               type_ == "BLOCK_RELAXATION" ||
-               type_ == "BLOCK RELAXATION" ||
-               type_ == "BLOCKRELAXATION"  ||
-               // Banded
-               type_ == "BANDED_RELAXATION" ||
-               type_ == "BANDED RELAXATION" ||
-               type_ == "BANDEDRELAXATION"  ||
-               // Tridiagonal
-               type_ == "TRIDI_RELAXATION"       ||
-               type_ == "TRIDI RELAXATION"       ||
-               type_ == "TRIDIRELAXATION"        ||
-               type_ == "TRIDIAGONAL_RELAXATION" ||
-               type_ == "TRIDIAGONAL RELAXATION" ||
-               type_ == "TRIDIAGONALRELAXATION") {
-      const std::string paramName = "relaxation: zero starting solution";
-      if (!params.isType<bool>(paramName) ||
-          (params.get<bool>(paramName) != InitialGuessIsZero)) {
-        paramList.set(paramName, InitialGuessIsZero);
-        SetPrecParameters(paramList);
-      }
+    if (prec_->supportsZeroStartingSolution()) {
+      prec_->setZeroStartingSolution(InitialGuessIsZero);
       supportInitialGuess = true;
-
-    } else if (type_ == "KRYLOV") {
-      const std::string paramName = "krylov: zero starting solution";
-      if (!params.isType<bool>(paramName) ||
-          (params.get<bool>(paramName) != InitialGuessIsZero)) {
-        paramList.set(paramName, InitialGuessIsZero);
-        SetPrecParameters(paramList);
-      }
-      supportInitialGuess = true;
-
     } else if (type_ == "SCHWARZ") {
       paramList.set("schwarz: zero starting solution", InitialGuessIsZero);
       //Because additive Schwarz has "delta" semantics, it's sufficient to