# Block-Tridiagonal Matrices with tridiagonal blocks

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

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?