Block Levinson solver

Efficiently solve symmetric block Toeplitz matrix equations using Block Levinson recursion

BLOCK_LEVINSON(Y, L) solves the matrix equation T * x = y, where T is a symmetric matrix with block Toeplitz structure, and returns the solution vector x. The matrix T is never stored in full (because it is large and mostly redundant), so the input parameter L is actually the leftmost "block column" of T (the leftmost d columns where d is the block dimension).

By Keenan Pepper; uploaded with permission.