Is there an efficient way to diagonalize a block tridiagonal $N\times N$ matrix of the following form:
\begin{matrix} A_0 & B & 0 & 0 & \ldots \\ B & A_1 & B & 0 & \ldots \\ 0 & B & A_2 & B & \ldots \\ \vdots & \vdots & \vdots & \vdots & \ddots&B \\ &&&&B&A_{K-1} \end{matrix}
where $A_i$ and $B$ are $L \times L$ matrices with $N = K\cdot L$ ? I need an algorithm that scales way better than the standard $\mathcal{O}(N^3)$ scaling, since I am usually interested in $N\approx 10.000$ and $K \approx L \approx 100$ and need to diagonalize a lot of them. I need both the eigenvalues AND the eigenvectors, if possible to full precision.
Further information, if useful: The matrix I am considering is hermitian, the matrices $B$ are just the unit matrix and the matrices $A_i$ are pentagonal. Furthermore I would in principle be interested in the same matrix with an additional $B$ matrix in the upper right and lower left corner, but I'd be grateful for any help. I've also read this paper:
https://arxiv.org/abs/1306.0217
but while it in principle provides an algorithm for my problem, I was not able to implement it efficiently in python and the authors did not publish their code. I would appreciate any answer to my problem very much! Thx!!