[eigen] A different stop criterion for the Conjugate Gradient solver

[Marcel Weiler <weiler@xxxxxxxxxxxxxxxxxxx>]

Hi,

I am currently implementing an algorithm for the particle-based 
simulation of elastic solids [1]. It requires the solution of a SPD 
linear system, for which the authors use a matrix-free Conjugate 
Gradient (CG) solver. I want to use the Eigen's CG solver since I 
successfully used it in other parts of my project before.
However, the authors of the paper use a different stopping criterion 
than the one Eigen uses. Eigen stops the CG iterations when the relative 
residual error falls below a given threshold, meaning |Ax-b|/|b| < 
tolerance. The authors of the paper, on the other hand, exploit the fact 
that in their linear system three consecutive entries of the solution 
vector correspond to the x-, y- and z-velocity of a particle. They stop 
the iteration when the "average absolute per particle [velocity] error 
of the system" reaches 1e-3. This error is computed as
1/n \sum_j^n |r_j|,
where r_j is the j'th Vector3d segmentof the solution vector. This 
results in a different overall behavior of the simulation.

How can I replace Eigen's stop criterion for the CG solver by a custom 
one without modifying ConjugateGradient.h? Should I write my own 
template specialization of void conjugate_gradient(...)? If so, can you 
provide me with an example (my previous attempts failed)?

I think a custom stop criterion might come in handy for other people as 
well, so might be a nice feature to include in Eigen itself.

I am looking forward to your suggestions!

Best,
Marcel

[1] "An Implicit SPH Formulation for Incompressible Linearly Elastic 
Solids", 
https://cg.informatik.uni-freiburg.de/publications/2017_CGF_elasticSolidsSPH.pdf

--
Marcel Weiler
Technische Universität Darmstadt
Graduate School of Computational Engineering
Dolivostraße 15, S4|10-204
64293 Darmstadt
Germany

Office: +49 6151 16 - 24387
Fax: +49 6151 16 - 24404