# 
# Numerical Mathematics and Computing, Fifth Edition
# Ward Cheney & David Kincaid
# Brooks/Cole Publ. Co.
# (c) 2003
#
# file:  penta
#
# Pentadiagonal Solver
# To test, take size n, either even or odd, and
# create a random pentadiagonal system.
# 
# Contents of vector  a,b,c,d,e,f will change so
# Original data kept in aa,bb,cc,dd,ee,ff
#
with(linalg):
# Select size of pentadiagonal matrix
 n:=8;
 randvector(n):a:=map(evalf,%):aa:=evalm(a):
 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):e:=map(evalf,%):ee:=evalm(e):
 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]:=ee[j] end if;
       if j=i-1 then M[i,j]:=aa[j] end if;
       if j=i then M[i,j]:=dd[i] end if;
       if j=i+1 then M[i,j]:=cc[i] end if;
       if j=i+2 then M[i,j]:=ff[i] end if;
     end do:
     M[i,n+1]:=bb[i]:
 end do:
 map(evalf,M);
# This is the procedure that solves the system.
# 
# 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]*a[2*i-3]):
    Xml[1,1]:=ood*(e[2*i-3]*d[2*i-2]-a[2*i-2]*a[2*i-3]):
    Xml[1,2]:=ood*(-e[2*i-3]*c[2*i-3]+a[2*i-2]*d[2*i-3]):
    Xml[2,1]:=ood*(-e[2*i-2]*a[2*i-3]):
    Xml[2,2]:=ood*(e[2*i-2]*d[2*i-3]):
# Temporary 2 x 2 storage matrix, since 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]:=a[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.`)
      end if:         
    # Done working, store results 
    d[2*i-1]:=t1[1,1]:
    c[2*i-1]:=t1[1,2]:
    a[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, record
    b[2*i-1]:=t2[1,1]:
    b[2*i]:=t2[2,1]:
 end do:
 if n/2<floor(n/2) then 
    ood:=1/(d[n-2]*d[n-1]-c[n-2]*a[n-2]):
    Xml[1,1]:=ood*(e[n-2]*d[n-1]-a[n-2]*a[n-1]):
    Xml[1,2]:=ood*(-e[n-2]*c[n-2]+a[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.`) end if:
    b[n]:=b[n]-Xml[1,1]*b[n-2]-Xml[1,2]*b[n-1]:
 end if:
# 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:
 end if:
# If n is even, then
 ood:=1/(d[nn-1]*d[nn]-c[nn-1]*a[nn-1]):
 X[nn-1]:=ood*(d[nn]*b[nn-1]-c[nn-1]*b[nn]):
 X[nn]:=ood*(-a[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]*a[i-1]):
    X[i-1]:=ood*(d[i]*t2[1,1]-c[i-1]*t2[2,1]):
    X[i]:=ood*(-a[i-1]*t2[1,1]+d[i-1]*t2[2,1]):
 end do:
# 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)));