Merge pull request #9645 from tjhei/precondition_chebyshev_report_eig

PreconditionChebyshev: report stats with output_details

doc/news/changes/minor/20200323TimoHeister

@@ -0,0 +1,3 @@
+New: PreconditionChebychev can now report information about eigenvalue and
+degree estimation when calling PreconditionChebychev::estimate_eigenvalues().
+<br> (Timo Heister, 2020/03/23)

include/deal.II/lac/precondition.h

@@ -914,7 +914,7 @@ class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
  * which are computed during the first invocation of vmult(). The algorithm
  * invokes a conjugate gradient solver (i.e., Lanczos iteration) so symmetry
  * and positive definiteness of the (preconditioned) matrix system are
- * requirements. The eigenvalue algorithm can be controlled by
+ * required. The eigenvalue algorithm can be controlled by
  * PreconditionChebyshev::AdditionalData::eig_cg_n_iterations specifying how
  * many iterations should be performed. The iterations are started from an
  * initial vector that depends on the vector type. For the classes
@@ -929,16 +929,20 @@ class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
  * length of the vector, except for the very first entry that is zero,
  * triggering high-frequency content again.
  *
- * The computation of eigenvalues happens the first time one of the
- * vmult(), Tvmult(), step() or Tstep() functions is called. This is because
- * temporary vectors of the same layout as the source and destination vectors
- * are necessary for these computations and this information gets only
- * available through vmult().
+ * The computation of eigenvalues happens the first time one of the vmult(),
+ * Tvmult(), step() or Tstep() functions is called or when
+ * estimate_eigenvalues() is called directly. In the latter case, it is
+ * necessary to provide a temporary vector of the same layout as the source
+ * and destination vectors used during application of the preconditioner.
  *
- * Due to the cost of the eigenvalue estimate in the first vmult(), this class
- * is most appropriate if it is applied repeatedly, e.g. in a smoother for a
- * geometric multigrid solver, that can in turn be used to solve several
- * linear systems.
+ * The estimates for minimum and maximum eigenvalue are taken from SolverCG
+ * (even if the solver did not converge in the requested number of
+ * iterations). Finally, the maximum eigenvalue is multiplied by a safety
+ * factor of 1.2.
+ *
+ * Due to the cost of the eigenvalue estimate, this class is most appropriate
+ * if it is applied repeatedly, e.g. in a smoother for a geometric multigrid
+ * solver, that can in turn be used to solve several linear systems.
  *
  * <h4>Bypassing the eigenvalue computation</h4>
  *
@@ -1160,6 +1164,31 @@ class PreconditionChebyshev : public Subscriptor
   size_type
   n() const;
 
+  /**
+   * A struct that contains information about the eigenvalue estimation
+   * performed by the PreconditionChebychev class.
+   */
+  struct EigenvalueInformation
+  {
+    /**
+     * Estimate for the smallest eigenvalue.
+     */
+    double min_eigenvalue_estimate;
+    /**
+     * Estimate for the largest eigenvalue.
+     */
+    double max_eigenvalue_estimate;
+    /**
+     * Number of CG iterations performed or 0.
+     */
+    unsigned int cg_iterations;
+    /**
+     * The degree of the Chebyshev polynomial (either as set using
+     * AdditionalData::degree or estimated as described there).
+     */
+    unsigned int degree;
+  };
+
   /**
    * Compute eigenvalue estimates required for the preconditioner.
    *
@@ -1172,7 +1201,7 @@ class PreconditionChebyshev : public Subscriptor
    * in AdditionalData, no computation is performed and the information given
    * by the user is used.
    */
-  void
+  EigenvalueInformation
   estimate_eigenvalues(const VectorType &src) const;
 
 private:
@@ -2242,20 +2271,24 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
 
 
 template <typename MatrixType, typename VectorType, typename PreconditionerType>
-inline void
+inline typename PreconditionChebyshev<MatrixType,
+                                      VectorType,
+                                      PreconditionerType>::EigenvalueInformation
 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
   estimate_eigenvalues(const VectorType &src) const
 {
   Assert(eigenvalues_are_initialized == false, ExcInternalError());
   Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
 
+  PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+    EigenvalueInformation info{};
+
   solution_old.reinit(src);
   temp_vector1.reinit(src, true);
 
   // calculate largest eigenvalue using a hand-tuned CG iteration on the
   // matrix weighted by its diagonal. we start with a vector that consists of
   // ones only, weighted by the length.
-  double max_eigenvalue, min_eigenvalue;
   if (data.eig_cg_n_iterations > 0)
     {
       Assert(data.eig_cg_n_iterations > 2,
@@ -2297,25 +2330,28 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
 
       // read the eigenvalues from the attached eigenvalue tracker
       if (eigenvalue_tracker.values.empty())
-        min_eigenvalue = max_eigenvalue = 1;
+        info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
       else
         {
-          min_eigenvalue = eigenvalue_tracker.values.front();
+          info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
 
           // include a safety factor since the CG method will in general not
           // be converged
-          max_eigenvalue = 1.2 * eigenvalue_tracker.values.back();
+          info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
         }
+
+      info.cg_iterations = control.last_step();
     }
   else
     {
-      max_eigenvalue = data.max_eigenvalue;
-      min_eigenvalue = data.max_eigenvalue / data.smoothing_range;
+      info.max_eigenvalue_estimate = data.max_eigenvalue;
+      info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
     }
 
   const double alpha = (data.smoothing_range > 1. ?
-                          max_eigenvalue / data.smoothing_range :
-                          std::min(0.9 * max_eigenvalue, min_eigenvalue));
+                          info.max_eigenvalue_estimate / data.smoothing_range :
+                          std::min(0.9 * info.max_eigenvalue_estimate,
+                                   info.min_eigenvalue_estimate));
 
   // in case the user set the degree to invalid unsigned int, we have to
   // determine the number of necessary iterations from the Chebyshev error
@@ -2324,7 +2360,7 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
   // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
   if (data.degree == numbers::invalid_unsigned_int)
     {
-      const double actual_range = max_eigenvalue / alpha;
+      const double actual_range = info.max_eigenvalue_estimate / alpha;
       const double sigma        = (1. - std::sqrt(1. / actual_range)) /
                            (1. + std::sqrt(1. / actual_range));
       const double eps = data.smoothing_range;
@@ -2333,16 +2369,18 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
         this)
         ->data.degree =
         1 + static_cast<unsigned int>(
-              std::log(1. / eps + std::sqrt(1. / eps / eps - 1)) /
+              std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
               std::log(1. / sigma));
+
+      info.degree = data.degree;
     }
 
   const_cast<
     PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
-    ->delta = (max_eigenvalue - alpha) * 0.5;
+    ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
   const_cast<
     PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
-    ->theta = (max_eigenvalue + alpha) * 0.5;
+    ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
 
   // We do not need the second temporary vector in case we have a
   // DiagonalMatrix as preconditioner and use deal.II's own vectors
@@ -2368,6 +2406,8 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
   const_cast<
     PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
     ->eigenvalues_are_initialized = true;
+
+  return info;
 }