From 9d149330831bc23b2fd9c03f2c439fe8237149f0 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 15 Aug 2015 10:46:57 +0200 Subject: [PATCH] Avoid setting the constant mode as a right hand side (works also in distributed computation case) Adjust tolerance and update test --- include/deal.II/lac/precondition.h | 15 +++-- .../step-37.with_lapack=true.output | 62 +++++++++---------- 2 files changed, 42 insertions(+), 35 deletions(-) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 48c4720e10..c21ead9852 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -1472,13 +1472,22 @@ PreconditionChebyshev::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 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(matrix.m())); + if (rhs->locally_owned_elements().is_element(0)) + (*rhs)(0) = 0.; + rhs->compress(VectorOperation::insert); internal::PreconditionChebyshev::EigenvalueTracker eigenvalue_tracker; SolverCG solver (control, memory); @@ -1494,8 +1503,6 @@ PreconditionChebyshev::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 2ef2bd9d95..76f9274e64 100644 --- 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 -- 2.39.5