# Fast solution of a heptadiagonal linear system

I have a linear system of the form $\mathbf{A}\mathbf{x} = \mathbf{f}$. If the length of the vector $\mathbf{x}$ is $N$, meaning that there are $N$ unknowns, then the matrix $\mathbf{A}$ has seven bands and the matrix is sparse. The matrix is not symmetric. It looks like this:

$\pmatrix{a_1&a_2&a_3&a_4&&&&&&\\ b_1&b_2&b_3&b_4&b_5&&&&&\\ c_1&c_2&c_3&c_4&c_5&c_6&&&&\\ d_1&d_2&d_3&d_4&d_5&d_6&d_7&&&\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&&\\ &&&&&&\\ &&&k_1&k_2&k_3&k_4&k_5&k_6&k_7&&&\\ &&&&l_1&l_2&l_3&l_4&l_5&l_6&l_7&&&\\ &&&&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&&\\ &&&&&&\\ &&&&&&&w_1&w_2&w_3&w_4&w_5&w_6&w_7\\ &&&&&&&&x_1&x_2&x_3&x_4&x_5&x_6\\ &&&&&&&&&y_1&y_2&y_3&y_4&y_5\\ &&&&&&&&&&z_1&z_2&z_3&z_4\\ }$

and the vector $\mathbf{f}$ looks like this: $\pmatrix{f_1\\f_2\\f_3\\f_4\\\vdots\\\\f_i\\f_{i+1}\\\vdots\\\\f_{N-3}\\f_{N-2}\\f_{N-1}\\f_{N}}$

I know that the general method is to use a LU decomposition and as the matrix is constant for me, I can pre-compute its decomposition and store it and use it to solve for the unknowns. But I would like to know if there is any method to solve such systems faster than the LU decomposition. Methods may not be restricted to direct methods only.

• The standard approach (implemented in LAPACK's gbsv and hence in SciPy and Matlab)) is to compute a LU decomposition using the knowledge that the matrix is banded and therefore most elements stay zero and can be skipped. In SciPy, this is exposed as scipy.linalg.solve_banded. (Note that this works for dense matrices; your matrix is sparse, so a general sparse LU factorization might be faster.) – Christian Clason Jan 11 '17 at 15:49
• @ChristianClason -- no, the banded solver is the fastest unless a significant number of entries within the heptadiagonal band are zero. If they are not, then the LU factors will be nonzero only in the band, the band of the LU factors will consist of nonzeros, and consequently computing as a banded matrix (rather than a sparse matrix) is definitely fastest. Note that any banded matrix by definition is sparse if the bandwidth is constant as the size of the matrix increases. – Wolfgang Bangerth Jan 12 '17 at 4:12
• @WolfgangBangerth Computationally, you are of course right; I was thinking more of the storage overhead if the (dense) matrix is large. But if the routine stores both the matrix and the factors as banded matrices (and not as arbitrary dense matrices), that will indeed always be faster. (I'm sure SciPy does; not sure about Matlab.) So we're back at the two general mantras: 1) Storage matters, and 2) If at all possible, use a high-quality library and not your own implementation. – Christian Clason Jan 12 '17 at 8:21
• @ChristianClason -- we're in violent agreement! :-) – Wolfgang Bangerth Jan 12 '17 at 18:28