The Setup

Using finite differences to discretize the 2d diffusion equation $$\partial_tu=\partial_x\left(A\partial_xu+B\partial_yu\right)+\partial_y\left(B\partial_xu+C\partial_yu\right)$$ we get a block-tridiagonal matrix ($n$ by $n$ blocks), in which each block is a tridiagonal ($m$ by $m$) matrix, with ($x$,$y$ dependent) diffusion coefficients $A,B,C$. These diffusion coefficients, and hence the whole matrix, only update every ~50 time steps. We have solved the resulting system with Lapack's banded LU matrix solver, as well as Intel's MKL sparse QR solver.

The Problem

Performence is bottlenecked by the matrix factorization. Because the matrix updates so rarely, those matrix decompositions could be skipped most of the time, increasing performance quite dramatically. Because we solve about 1000 of these systems per time step, each with about $n=200$ diagonal-blocks, with each $m=201$ diagonal elements, storing the decompositions in memory isn't an option. I think the LU decomposition is $O(n^2m)$ in space, and experiments with storing the sparse MKL handles after the factorization step also aren't very promising.

I am looking for a way to reuse those computations done for the same matrix, where only the right hand side changes. I considered the block version of the Thomas algorithm, but here we would have to solve tridiagonal systems with tridiagonal matrix right hand sides. Those solutions are dense $m$ by $m$ matrices. Storing these would bring me back to $O(nm^2)$.

Maybe there is a way of tackling this that I'm missing here. Any suggestion would be very welcome.


1 Answer 1


From my quick experiments in python, I find that a LU decomposition with the permutation strategy MMD_AT_PLUS_A¹ yields to an $\mathcal{O}(12n^{2.268})$ number of nonzeros. With m=n=200, this is reduces the weight of 1000 LU's to about 15GB of memory, which I think should be manageable on a modern workstation? Would that work for you?

enter image description here

¹ https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html

  • $\begingroup$ Thanks! This made me take a second look at the MKL Pardiso solver which uses a sparse LU decomposition, just as scipy, instead of QR. When parallelizing Pardiso with OpenMP I still get pretty hefty memory costs (would be around 40GB instead of your proposed 15), but it seems Pardiso offers an option to store the decompositions on disk. I shall give that a go. $\endgroup$
    – Michael
    May 25, 2023 at 16:06

Your Answer

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

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