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

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.

• Kahan(-Babuska-Neumaier) summation algorithm may help: en.wikipedia.org/wiki/Kahan_summation_algorithm . Probably, only the dot product computation is wrong. (If this doesn't work, you can try ORO (Ogita-Rump-Oishi) dot product as well.) Nov 20, 2021 at 22:34