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.
