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

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