# Symmetric Pentadiagonal Solver
# To test, take size n, either even or odd, and 
# create a random symmetric pentadiagonal system.
# 
# Contents of vector  a,b,c,d,e,f will change,  
# so original data saved in aa,bb,cc,dd,ee,ff
 with(linalg):
# select size of symmetric pentadiagonal matrix
 n:=8;
 randvector(n):b:=map(evalf,%):bb:=evalm(b):
 randvector(n):c:=map(evalf,%):cc:=evalm(c):
 randvector(n):d:=map(evalf,%):dd:=evalm(d):
 randvector(n):f:=map(evalf,%):ff:=evalm(f):
# Fix original augmented matrix as M.
 M:=matrix(n,n+1,0):
 for i from 1 to n do
    for j from 1 to n+1 do
       if j=i-2 then M[i,j]:=ff[j] fi;
       if j=i-1 then M[i,j]:=cc[j] fi;
       if j=i then M[i,j]:=dd[i] fi;
       if j=i+1 then M[i,j]:=cc[i] fi;
       if j=i+2 then M[i,j]:=ff[i] fi;
     od:
     M[i,n+1]:=bb[i]:
 od:
 map(evalf,M);
# This is the procedure that solves the system.
# 
# Remark: The symmetric solver was created from the nonsymmetric version
# by changing each a to c, each e to  f.
# 
# forward solve
 for i from 2 to floor(n/2) do
    # ood = one over determinant; then use adjoint matrix
    ood:=1/(d[2*i-3]*d[2*i-2]-c[2*i-3]*c[2*i-3]):
    Xml[1,1]:=ood*(f[2*i-3]*d[2*i-2]-c[2*i-2]*c[2*i-3]):
    Xml[1,2]:=ood*(-f[2*i-3]*c[2*i-3]+c[2*i-2]*d[2*i-3]):
    Xml[2,1]:=ood*(-f[2*i-2]*c[2*i-3]):
    Xml[2,2]:=ood*(f[2*i-2]*d[2*i-3]):
    # temporary 2 x 2 storage matrix, still working with elements in block
    t1[1,1]:=d[2*i-1]-Xml[1,1]*f[2*i-3]-Xml[1,2]*c[2*i-2]:
    t1[1,2]:=c[2*i-1]-Xml[1,2]*f[2*i-2]:
    t1[2,1]:=c[2*i-1]-Xml[2,1]*f[2*i-3]-Xml[2,2]*c[2*i-2]:
    t1[2,2]:=d[2*i]-Xml[2,2]*f[2*i-2]:
    if t1[1,1]*t1[2,2]-t1[1,2]*t1[2,1]=0 then ERROR(`Singular matrix.`)
 fi:         
    # done working, store results 
    d[2*i-1]:=t1[1,1]:
    c[2*i-1]:=t1[1,2]:
    c[2*i-1]:=t1[2,1]:
    d[2*i]:=t1[2,2]:
    # temporary 2 x 1 storage matrix
    t2[1,1]:=b[2*i-1]-Xml[1,1]*b[2*i-3]-Xml[1,2]*b[2*i-2]:
    t2[2,1]:=b[2*i]-Xml[2,1]*b[2*i-3]-Xml[2,2]*b[2*i-2]:
    # done working, store results
    b[2*i-1]:=t2[1,1]:
    b[2*i]:=t2[2,1]:
 od:
 if n/2<floor(n/2) then 
    ood:=1/(d[n-2]*d[n-1]-c[n-2]*c[n-2]):
    Xml[1,1]:=ood*(f[n-2]*d[n-1]-c[n-2]*c[n-1]):
    Xml[1,2]:=ood*(-f[n-2]*c[n-2]+c[n-1]*d[n-2]): 
    d[n]:=d[n]-Xml[1,1]*f[n-2]-Xml[1,2]*c[n-1]:
    if d[n]=0 then ERROR(`Singular matrix.`) fi:
    b[n]:=b[n]-Xml[1,1]*b[n-2]-Xml[1,2]*b[n-1]:
 fi:
# start backsolve
 nn:=n:
# if n is odd, then
 if n/2<floor(n/2) then 
    X[n]:=b[n]/d[n]:
    b[n-2]:=b[n-2]-X[n]*f[n-2]:
    b[n-1]:=b[n-1]-X[n]*c[n-1]:
    nn:=n-1:
 fi:
# if n is even, then
 ood:=1/(d[nn-1]*d[nn]-c[nn-1]*c[nn-1]):
 X[nn-1]:=ood*(d[nn]*b[nn-1]-c[nn-1]*b[nn]):
 X[nn]:=ood*(-c[nn-1]*b[nn-1]+d[nn-1]*b[nn]):
 for i from (nn-2) by -2 to 2 do
    t2[1,1]:=b[i-1]-f[i-1]*X[i+1]:
    t2[2,1]:=b[i]-c[i]*X[i+1]-f[i]*X[i+2]:
    ood:=1/(d[i-1]*d[i]-c[i-1]*c[i-1]):
    X[i-1]:=ood*(d[i]*t2[1,1]-c[i-1]*t2[2,1]):
    X[i]:=ood*(-c[i-1]*t2[1,1]+d[i-1]*t2[2,1]):
 od:
# Show the solution given by the procedure.
 print(X);
# Check the procedure by Maple's linsolve to see if it gives the same
# answer.
 linsolve(submatrix(M,1..n,1..n),submatrix(M,1..n,(n+1)..(n+1)));