Eigendecomposition of a block tridiagonal matrix

Question

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 provide a minimal reproducible example (in python) with the precise structure of the matrix I am interested in:

import scipy.linalg as linalg
import numpy as np

#N is usually in this form and I want to be able to have a code that is still fast, if  I set the 5 to a 50.
N = 2* 5**2

#Creating Matrix in my form...

N_band_max = int(np.sqrt(N/2))*2

#setting bands...
Main = np.random.random(N)
Diag1 = np.random.random(N-1) + 1j*np.random.random(N-1)
Diag2 = np.ones(N-2)
Diag3 = np.ones(N-N_band_max) 

#This is for the block tridiagonal structure

for i in range(len(Diag2)):
    if(i%N_band_max == N_band_max-2 or i%N_band_max == N_band_max-1):
        Diag2[i] = 0.
        
for i in range(len(Diag1)):
    if(i%2 == 1 ):
        Diag1[i] = 0.


H = np.diag(Main) + np.diag(Diag1,1) + np.diag(Diag1.conjugate(),-1) + np.diag(Diag2,2) + np.diag(Diag2.conjugate(),-2) + np.diag(Diag3,N_band_max) + np.diag(Diag3.conjugate(),-N_band_max)


#This is the part I want to be FAAASTER
EigVals,EigVecs = linalg.eigh(H)

I would appreciate any answer to my problem very much! Thx!