Solving a symmetric banded linear system with Cholesky factorisation

Question

I'm trying to solve a linear system Ax=b, where A has been factorised by Cholesky decomposition, and A is stored using symmetric banded storage.

Given a symmetric, positive definite matrix, A, I've managed to implement a function that first transforms the matrix A into a symmetric banded format via the transformation A[i,j] = Ab[i-j,j] with j <= i < min(n,j+1+p) and p is the lower bandwidth. As an example, inputting A, we get output Ab

  [[1, 1, 0],   Ab=[[1, 2, 3],
A=[1, 2, 4],      [1, 4, 0]]
  [0, 4, 3]]

I have also implemented a function that applies Cholesky decomposition to Ab. I now need to implement forward and backward substitution to the banded Cholesky matrix. Given L (banded Cholesky factorisation of Ab) and b, I need to solve Lx=b, and I have based my algorithm off of a non-banded forward-backward substitution algorithm that is tested and works.

My current banded algorithm is as follows, but I think my indices for the array slicing are wrong.

(n,m) = L.shape

# forward sub
for j in range(m):
    m = min(j+1+p,n)
    x[j] = x[j]/L[0,j]
    x[(j+1):m] = x[(j+1):m] - L[(j+1):m,j]*x[j]

# back sub
for j in range(n-1, -1, -1):
    m = max(0,j-p)
    x[j] = x[j]/L[0,j]
    diagonal = np.fliplr(L).diagonal(n-j-1)[1:]
    x[m:j] = x[m:j] - np.flip(diagonal*b[j])
return x

---

EDIT

I've updated the forward substitution algorithm to use a for loop instead of slicing.

for j in range(n):
    c = min(j+1+p,n)
     x[j] = x[j]/L[0,j]
     temp = x[j]
     for i in range(j+1,c):
          x[i] = x[i] - L[(j+i),j]*x[j]

Testing this code using A (thus p=1), and b (row sum, so x should be a vector of ones)

  [[12.0, 4.0, 0.0],   
A=[4.0, 2.0, -1.0],   b=[16.0, 5.0, 6.0]
  [0.0, -1.0, 7.0]]

Using non working symmetrically-stored algorithms we get the following outputs L (Cholesky factorisation) and x

  [[ 3.46410162,  0,          0        ]
L=[ 1.15470054,  0.81649658  0         ]   x=[1,1,1]
  [ 0,         -1.22474487,  2.34520788]]

and the output of the forward substitution is x_forward=[4.61880215,-0.40824829,2.34520788]. Testing this now with the symmetrically stored algorithm we get the Cholesky factorisation we get

L_sym=[[ 3.46410162,0.81649658,2.34520788]
       [ 1.15470054,-1.22474487,0        ]]

as expected, but the output of the forward substitution is x_sym_forward=[4.61880215,-0.40824829,6] so my updated algorithm seems to compute the first two entries correctly but then fails to compute the final entry. I've yet to attempt fixing the back substitution algorithm.