Avoid setting the constant mode as a right hand side

(works also in distributed computation case)
Adjust tolerance and update test

diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index 48c4720e1082e892b38b1bf1e4ac25305973aa01..c21ead9852bbaecf1db69178c84a338e30c2203e 100644 (file)
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -1472,13 +1472,22 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
               ExcMessage ("Need to set at least two iterations to find eigenvalues."));
 
       // set a very strict tolerance to force at least two iterations
-      ReductionControl control (data.eig_cg_n_iterations, 1e-20, 1e-20);
+      ReductionControl control (data.eig_cg_n_iterations, 1e-35, 1e-10);
       GrowingVectorMemory<VECTOR> memory;
       VECTOR *rhs = memory.alloc();
       VECTOR *dummy = memory.alloc();
-      rhs->reinit(data.matrix_diagonal_inverse, true);
+      rhs->reinit(data.matrix_diagonal_inverse);
       dummy->reinit(data.matrix_diagonal_inverse);
+
+      // heuristically, a right hand side close to a constant has been shown
+      // to quickly reveal the largest eigenvalue. however, avoid to use the
+      // exact constant because that might be not in the range space of some
+      // matrices (purely Neumann matrices with constant mode filtered out by
+      // orthogonal projection in the matrix-vector product)
       *rhs = 1./std::sqrt(static_cast<double>(matrix.m()));
+      if (rhs->locally_owned_elements().is_element(0))
+        (*rhs)(0) = 0.;
+      rhs->compress(VectorOperation::insert);
 
       internal::PreconditionChebyshev::EigenvalueTracker eigenvalue_tracker;
       SolverCG<VECTOR> solver (control, memory);
@@ -1494,8 +1503,6 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
       catch (SolverControl::NoConvergence &)
         {
         }
-      Assert(control.last_step() >= 2,
-             ExcMessage("Could not find eigenvalues"));
 
       memory.free(dummy);
       memory.free(rhs);

diff --git a/tests/matrix_free/step-37.with_lapack=true.output b/tests/matrix_free/step-37.with_lapack=true.output
index 2ef2bd9d95db6a844b6cc7129ee251b5e7441228..76f9274e64735ef43805fd1198fffd85c18ae8c4 100644 (file)
--- a/tests/matrix_free/step-37.with_lapack=true.output
+++ b/tests/matrix_free/step-37.with_lapack=true.output
@@ -1,37 +1,37 @@
 
 DEAL:2d::Cycle 0
 DEAL:2d::Number of degrees of freedom: 81
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Convergence step 3 value 0
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Failure step 10 value 0
-DEAL:2d:cg::Starting value 1.00
+DEAL:2d:cg::Starting value 0.943
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.980
+DEAL:2d:cg::Convergence step 8 value 0
+DEAL:2d:cg::Starting value 0.994
 DEAL:2d:cg::Failure step 10 value 0.00110
 DEAL:2d:cg::Starting value 0.132
 DEAL:2d:cg::Convergence step 4 value 0
 DEAL:2d::
 DEAL:2d::Cycle 1
 DEAL:2d::Number of degrees of freedom: 289
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Convergence step 3 value 0
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Failure step 10 value 0
-DEAL:2d:cg::Starting value 1.00
+DEAL:2d:cg::Starting value 0.943
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.980
+DEAL:2d:cg::Convergence step 8 value 0
+DEAL:2d:cg::Starting value 0.994
 DEAL:2d:cg::Failure step 10 value 0.00110
-DEAL:2d:cg::Starting value 1.00
+DEAL:2d:cg::Starting value 0.998
 DEAL:2d:cg::Failure step 10 value 0.153
 DEAL:2d:cg::Starting value 0.0677
 DEAL:2d:cg::Convergence step 5 value 0
 DEAL:2d::
 DEAL:2d::Cycle 2
 DEAL:2d::Number of degrees of freedom: 1089
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Convergence step 3 value 0
-DEAL:2d:cg::Starting value 1.00
-DEAL:2d:cg::Failure step 10 value 0
-DEAL:2d:cg::Starting value 1.00
+DEAL:2d:cg::Starting value 0.943
+DEAL:2d:cg::Convergence step 2 value 0
+DEAL:2d:cg::Starting value 0.980
+DEAL:2d:cg::Convergence step 8 value 0
+DEAL:2d:cg::Starting value 0.994
 DEAL:2d:cg::Failure step 10 value 0.00110
-DEAL:2d:cg::Starting value 1.00
+DEAL:2d:cg::Starting value 0.998
 DEAL:2d:cg::Failure step 10 value 0.153
 DEAL:2d:cg::Starting value 1.00
 DEAL:2d:cg::Failure step 10 value 1.07
@@ -40,31 +40,31 @@ DEAL:2d:cg::Convergence step 5 value 0
 DEAL:2d::
 DEAL:3d::Cycle 0
 DEAL:3d::Number of degrees of freedom: 125
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Convergence step 3 value 0
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Failure step 10 value 0
+DEAL:3d:cg::Starting value 0.981
+DEAL:3d:cg::Convergence step 2 value 0
+DEAL:3d:cg::Starting value 0.996
+DEAL:3d:cg::Convergence step 10 value 0
 DEAL:3d:cg::Starting value 0.125
 DEAL:3d:cg::Convergence step 5 value 0
 DEAL:3d::
 DEAL:3d::Cycle 1
 DEAL:3d::Number of degrees of freedom: 729
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Convergence step 3 value 0
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Failure step 10 value 0
-DEAL:3d:cg::Starting value 1.00
+DEAL:3d:cg::Starting value 0.981
+DEAL:3d:cg::Convergence step 2 value 0
+DEAL:3d:cg::Starting value 0.996
+DEAL:3d:cg::Convergence step 10 value 0
+DEAL:3d:cg::Starting value 0.999
 DEAL:3d:cg::Failure step 10 value 0.000695
 DEAL:3d:cg::Starting value 0.0479
 DEAL:3d:cg::Convergence step 4 value 0
 DEAL:3d::
 DEAL:3d::Cycle 2
 DEAL:3d::Number of degrees of freedom: 4913
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Convergence step 3 value 0
-DEAL:3d:cg::Starting value 1.00
-DEAL:3d:cg::Failure step 10 value 0
-DEAL:3d:cg::Starting value 1.00
+DEAL:3d:cg::Starting value 0.981
+DEAL:3d:cg::Convergence step 2 value 0
+DEAL:3d:cg::Starting value 0.996
+DEAL:3d:cg::Convergence step 10 value 0
+DEAL:3d:cg::Starting value 0.999
 DEAL:3d:cg::Failure step 10 value 0.000695
 DEAL:3d:cg::Starting value 1.00
 DEAL:3d:cg::Failure step 10 value 0.0952