This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

No (reasonable) default stopping criteria for CG in Eigen?

0 votes

Solving magnetostatic problems, I found some strange behaviour concerning
stopping criteria for Conjugate Gradients (and may be also for other Krylov methods)
using the Eigen linear algebra backend in version 2016.1.
Starting the solver with

problem = LinearVariationalProblem(...)
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = 'cg'
solver.parameters["preconditioner"] = 'default'
solver.solve()

and solver.parameters['krylov_solver']

('absolute_tolerance', None)
('divergence_limit', None)
('error_on_nonconvergence', None)
('maximum_iterations', None)
('monitor_convergence', None)
('nonzero_initial_guess', None)
('relative_tolerance', None)
('report', None)

always ends up with

*** Warning: Krylov solver did not converge in ... iterations

The same code works fine with PETSc in version 2016.1 (and also
older versions).
Explicitely setting

solver.parameters["krylov_solver"]["relative_tolerance"]

to some reasonable value leads to a proper behaviour of the solver.
It seems, that in the case of Eigen there is no reasonable default
stopping criterium defined whereas for PETSc it is.
This is pretty strange to me.

Is there any reason for this difference between Eigen and PETSc?

This is the first time I am using Eigen, so may be this difference
is something well known to the 'experienced user'.

asked Jul 1, 2016 by mre FEniCS Novice (240 points)

1 Answer

+1 vote

If you don't specify the solver parameters you will get backend default settings. For PETSc, the default relative tolerance is $10^{-5}$. For Eigen, the default tolerances are apparently much more strict.

Also, you've left the preconditioner to backed-dependent defaults. For PETSc, this will be incomplete LU, for Eigen it seems to be one Jacobi iteration.

answered Jul 1, 2016 by Magne Nordaas FEniCS Expert (13,820 points)

Thank you for your explanation.
Does there exist some documentation concerning the default values? I checked the online-documentation of Eigen (version 3.2.8), which states "The defaults are the size of the problem for the maximal number of iterations and NumTraits::epsilon() for the tolerance."
In my "non-converging" examples, the iteration is stopped after 2*n steps, where n is the
dimension of the system. So obviously, the defaults used there are not identical to those
listed in the Eigen-documentation.

...