3
$\begingroup$

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).

Improve this question
$\endgroup$
8
  • $\begingroup$ How do you determine the sparsity pattern for your incomplete decomposition? Is it substantially larger than for your original matrix? Because if it isn't, then applying one decomposition is equivalent in work to one matrix-vector product. $\endgroup$
    – Wolfgang Bangerth
    Oct 22, 2021 at 22:27
  • 1
    $\begingroup$ You should experiment with reorderings of your matrix. It is entirely possibly that a perfect elimination order exists or that an iincomplete Cholesky factorization will produce a better preconditioner for the reordered system. If this does not produce a satisfactory solution, then combine reorderings with different nonzero drop tolerances for the incomplete factorization. If this does not work out, then please add a link to the matrix. $\endgroup$
    – Carl Christian
    Oct 23, 2021 at 0:28
  • $\begingroup$ Perhaps it is implied, but are you using a fill in reducing permutation? That can make a huge difference for the matrix vector products and the triangular solves (if you increase the sparsity of the preconditioner) $\endgroup$
    – Charlie S
    Oct 23, 2021 at 0:37
  • $\begingroup$ Carl beat me to it! $\endgroup$
    – Charlie S
    Oct 23, 2021 at 0:41
  • $\begingroup$ @WolfgangBangerth I am using IC(0) which means that I use the exactly same sparsity pattern as for the system matrix $A$. That is I re-use the CSR sparsity arrays from $A$. Applying the Cholesky factorization is cheap, applying the forward and backsubstitution is what is slowing down things in each iteration. $\endgroup$
    – lightxbulb
    Oct 23, 2021 at 7:05

1 Answer 1

Reset to default
3
$\begingroup$

Iterative solvers are (almost always) memory-bandwidth bound. So, how you access memory is very important.

Modern computers work in cache lines, where a set of contiguous bytes will be fetched together, regardless of how many of those bytes will be accessed. On Intel machines, for example, cache lines are 64 bytes. That's eight 64-bit floats, alternatively, sixteen ints or 32-bit floats.

From your comment

That is I re-use the CSR sparsity arrays from A.

you're going to have cachelines shared between $L$ and $L^T$ (all of them if you have less than 8 nonzeros per row). So, each of your triangular solves will need as much data from the main memory as the original matrix-vector product. Thus, the $3\times$ increase in per iteration cost.

If you only store $L$ in a lower-triangular sparsity pattern, you'll halve data movement for each of the triangular solves. This would put you at about the performance of 60 iterations of the Jacobi-preconditioned solver.

Improve this answer
$\endgroup$
9
  • $\begingroup$ So you would suggest keeping a purely lower and upper triangular (sparse) arrays and using that for the forward and backward substitutions instead? Did I understand this correctly? $\endgroup$
    – lightxbulb
    Oct 23, 2021 at 16:35
  • 1
    $\begingroup$ Yes, that's the gist of it. $\endgroup$
    – Neil Lindquist
    Oct 23, 2021 at 18:48
  • 1
    $\begingroup$ The back-solve is: for $i = n, n-1, \dots, 1$ do $b[i] /= L^T[i, i]$; $b[0:i-1] -= L^T[0:i-1, i]b[i] = L[i, 0:i-1]b[i]$. So, at the $i$th step, you're only accessing the $i$th row of $L$ (or column of $L^T$) which is contiguous in memory. $\endgroup$
    – Neil Lindquist
    Oct 23, 2021 at 21:21
  • 1
    $\begingroup$ Thank you. I hadn't realized I could rewrite this to solve a tridiagonal system with an optimized access pattern (in hindsight I use a similar rewriting to evaluate $A^T x$ products when I only have $A$). I'll implement these, test them, and see what performance I get. $\endgroup$
    – lightxbulb
    Oct 23, 2021 at 21:47
  • 1
    $\begingroup$ I implemented the Cholesky only while storing $L$ and with the suggested backsubstitution. Unfortunately it doesn't seem to have been a memory access issue since the performance is almost the same, i.e. I assume it's the FLOPs that is the issue, and likely the fact that the backsubstitution is not parallelized (I am on a laptop CPU, so multiple weak cores). Nevertheless your suggestion were very useful for reducing memory allocation (and access), and I plan to also make my symmetric mat-vec multiplies use a triangular matrix. As far as parallelization goes, the paper I linked may be an option. $\endgroup$
    – lightxbulb
    Oct 24, 2021 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.