I am currently solving the harmonic equation using a P1 FEM discretisation. The resulting matrix $A$ is SPD and fairly sparse so I use a preconditioned conjugate gradients (CG) solver to find a solution. I have tried to accelerate the convergence both through a Jacobi preconditioner and an incomplete Cholesky preconditioner (using the sparsity pattern of $A$: IC(0)). Both work great for decreasing the number of CG iterations substantially, however even though the Choleksy preconditioned iterations are 3 times fewer (30) compared to the Jacobi preconditioned iterations (90) it is nevertheless a little slower. I have narrowed this down to the solving of the triangular systems $L y = x, DL^T z = y$, which seem to account for the 3 times increase of runtime per iteration (I am assuming this is because they are inherently unparallelizable by the compiler). My question is whether there are ways that this can be sped up. I have read about the Eisenstat trick in the context of SSOR, and in Saad's book the following is mentioned:
Also, for certain types of matrices, the IC(0) preconditioner can be expressed in this manner, where D can be obtained by a recurrence formula.
though I haven't found more details on this. From my profiling the computation of the incomplete Cholesky factorisation is negligible in term of runtime. Are there any relevant optimisations I can do? Is it a good idea to look into incomplete Cholesky factorisation with fill-in instead? Or maybe it would be better to focus on an algebraic multigrid preconditioner?
Edit: I have added an image of a zoom-in on the diagonal of a similar matrix to the one I have: The original system is 16 million by 16 million elements so I cannot store the full image of it, but I have tested several smaller variants and all the non-zero entries are clustered around the diagonal, which is to be expected since vertices in the FEM mesh are enumerated left to right and top to bottom (spatially), so there is at least some ordering, although I doubt it's optimal (also Dirichlet nodes are taken out since they only affect the right-hand side).
Edit 2: It seems like this is a known bottleneck, but the solutions to it are by no means simple: https://arxiv.org/pdf/1908.00741.pdf I will probably look into multigrid instead, after trying how an RCM reordering will affect the runtime.
Edit 3: I implemented the suggestions from Neil Lindquist, however those didn't produce a notable improvement in runtime, but they did reduce memory requirements and memory access substantially. I am assuming that for my problem I am limited by the computations and not the memory access, and mainly the lack of parallelization for the triangular solves. For this I guess I should refer to the linked paper, or use a domain decomposition or multigrid preconditioner instead (since both can be parallelized).