ICCG negative residual products $r^TM^{-1}r$

Question

I have a linear system $Ax=b$ resulting from a finite element discretization of the Poisson equation. I am applying an IC0 (incomplete Cholesky ($LDL^T$) with the same sparsity as the original matrix) preconditioned conjugate gradients solver to it. However, when this system grows large enough I occasionally get negative preconditioned residual products $r_k^TM^{-1}r_k$ (where $r_k = b-Ax_k$, and $M=LDL^T$ is the preconditioner) in some iterations.

I am assuming that this is caused by floating point error (and I use doubles at this point so using higher precision is not currently an option) either during the construction of the incomplete $LDL^T$ decomposition or during the forward and backsubstitution. I use the Cholesky–Banachiewicz algorithm for the $LDL^T$ decomposition (the $LL^T$ decomposition fails due to negative square roots due to floating point precision). I have ensured manually that $A$ is exactly symmetric, and I have also tried using the update $r_k = b-Ax_k$ instead of $r_k = r_{k-1} - \alpha Aq_k$ in each PCG iteration, but this didn't rectify the issue. Are there any ways that this can be made more robust? Or any other suggestions regarding preconditioners that would not involve such large errors? I am currently considering domain decomposition and block Jacobi as it would allow me to perform the decomposition on multiple smaller systems.

Answer 1

For reference on what I did: I tried more precise summation algorithms, however this helped only to some extent and having a slightly larger system or a slightly more ill-conditioned problem still resulted in CG diverging. I have since switched to geometric multigrid and it works extremely well even in single precision. It also turned out to be much faster than Choleksy preconditioned CG, and there are extensions that make it work for boundaries with more complex shape: A parallel multigrid Poisson solver for fluids simulation on large grids