restart:
with(LinearAlgebra):

KilicInv := proc(x::Array, y::Array, z::Array, delta, mu)

local n, i, j, jval, ainv, yfact, zfact, p, t, q, ordsign:
local vsign, fact2, fact3, fact4, fact5, temp, sing, det:

n := upperbound(x): # get inverse bound

# Check the rest of the supplied data is consistent
if upperbound(y) <> n-1 or upperbound(z) <> n-1 then
error "x, y and z arrays are of incompatible lengths"
end if:

# Acquire required space for locally used arrays and
# initialize them all to zero
yfact := Array(0..n-1): zfact := Array(0..n-1):
p := Array(-1..n-1): t := Array(1..n): q := Array(1..n):
vsign := Array(1..n+1):
fact2 := Array(1..n): fact3 := Array(1..n):
fact4 := Array(1..n): fact5 := Array(1..n):
ainv := Matrix[1..n,1..n]:

# Assume non-singular matrix
sing := false:

if n mod 2 = 0 then
ordsign := 1:
else
ordsign := -1:
end if:

yfact[0], zfact[0], p[-1], p[0] := 1, 1, 1, x[1]:

for i from 1 to n-1 do
yfact[i] := y[i]*yfact[i-1]:
zfact[i] := z[i]*zfact[i-1]:
p[i] := x[i+1]*p[i-1] - y[i]*z[i]*p[i-2]:
end do:

t[n], q[n] := 0, ordsign:

for i from n-1 by -1 to 1 do
t[i] := t[i+1] + (yfact[i-1]*zfact[i-1])/(p[i-1]*p[i-2]):
q[i] := (delta*p[n-2]*t[i] + ordsign*yfact[n-1])/yfact[i-1]:
end do:

det := p[n-1] - delta*mu*p[n-2]*t[1] -
ordsign*(mu*yfact[n-1] + delta*zfact[n-1]):

if det = 0 then
sing := true;
return:
end if:

fact5[1..n] := det /~ yfact[0..n-1]:
fact2[1..n] := mu *~ q[1..n]:
fact3[1..n] := (ordsign*zfact[n-1]/p[n-2]) *~ q[1..n]:

fact4[1..n-1] := fact5[1..n-1] + fact2[1..n-1]:
fact5[1..n] := (fact5[1..n] *~ t[1..n]) + fact3[1..n]:

for i from 1 by 2 to n+1 do
vsign[i] := 1:
end do:
for i from 2 by 2 to n+1 do
vsign[i] := -1:
end do:


(*
vsign[1..n+1, 2] := 1:
vsign[2..n+1, 2] := -1:
*)

for j from 2 to n do
for i from 1 to j-1 do
ainv[i, j] := t[j]*fact4[i] + fact3[i]:
end do
end do:

for j from 1 to n do
for i from j to n do
ainv[i, j] := t[j]*fact2[i] + fact5[i]:
end do
end do:

jval := 0:
fact5[1..n] := p[-1..n-2]/det:

for j from 1 to n do
temp:= (p[j-2]/zfact[j-1]):
for i from 1 to n do
ainv[i, j] := temp*ainv[i,j] * vsign[i+jval] * fact5[i] :
end do:
jval := 1 - jval:
end do:

return eval(sing), simplify(det), eval(ainv):

end proc:

processTimes:=proc(MapleTimes::Array, OtherTimes::Array,
                   n::integer)
# Produce the desired summary of the execution times obtained
# by running the test matrix given to testInv nRepeats times
local nRep, i: 
nRep := upperbound(MapleTimes): 
print(n, max(MapleTimes), min(MapleTimes),
      add(MapleTimes)/nRep, max(OtherTimes),
      min(OtherTimes), add(OtherTimes)/nRep):
end proc:

testInv:=proc(x::Array, y::Array, z::Array, delta, mu,
nRepeats, invTest::procedure, descript::string)
#
# Invert A using the Maple library procedure and the
# method provided by the user supplied procedure
# invTest. The string descript is used to label the output
# reporting on the procedure invTest.
#
# Check that the max element of the difference in the
# two inverses is zero. This assumes that the difference
# isn't an algebraic quantity which cannot be simplified
# to zero for some reason.
#
# Return 0 if matrices are equal and 1 otherwise.
#
# In the case of non-equality print the two inverses
#
local A, AinvMaple, Ainv, mval, st1, st2, M, n, i, j, 
      Adet, det, sing, MapleTimes, OtherTimes:

n := upperbound(x):
A := Matrix(1..n, 1..n);
Ainv := Matrix(1..n, 1..n);
AinvMaple := Matrix(1..n, 1..n);

MapleTimes := Array(1..nRepeats);
OtherTimes := Array(1..nRepeats);

for i from 1 to n-1 do
A[i,i], A[i,i+1], A[i+1,i] := x[i], y[i], z[i]:
end do:
A[n,n], A[1,n], A[n,1] := x[n], delta, mu:

for i from 1 to nRepeats do
st1:=time():
Adet := Determinant(A):
AinvMaple:=MatrixInverse(A):
st1:=time() - st1;
# print("Time using Maple library routine, MatrixInverse: ",st1);

st2:=time():
sing,det,Ainv:=invTest(x,y,z,delta,mu):
st2:=time() - st2;
# print(cat("Time using procedure ", descript, ":",st2));

MapleTimes[i] := st1:
OtherTimes[i] := st2:
end do:

mval := 0:
# print ("Test for matrix of order ",n); 

#----------------------
print("MAPLE:"):
print(AinvMaple[1,1]):
print("Kilic:"):
print(Ainv[1,1]):
#----------------------

for i from 1 to n do
for j from 1 to n do
if evalb(simplify(AinvMaple[i,j] - Ainv[i,j]) <> 0) then
mval := mval + 1:
print ("[",i,",",j,"] entries differ:");
print (AinvMaple[i,j]); print(Ainv[i,j]);
end if:
end do:
end do: 

if mval <> 0 then
print ("Test results differ");
print ("Inverse by Maple procedure, MatrixInverse:");
print (AinvMaple);
print (cat("Inverse by ", descript, " procedure:"));
print (Ainv);
else
 print ("Inverses agree");
end if:

if evalb(simplify(det-Adet) = 0) then
# print("Determinants agree");
else
print("Determinants differ");
print("Using Maple procedure, Determinant: ", Adet);
print(cat("Determinant returned by procedure ",descript, " :",det)):
end if:

processTimes(MapleTimes, OtherTimes, n):

end proc:

# Test using the tridiagonal, inverse Lehmer matrix with
# rational elements
Test1 := proc(invTest::procedure, descript::string,
nstart::integer, nend::integer, nstep::integer,
nRepeats::integer)
local n, x, y, z, delta, mu, i:

# output test routine header 
print(): 
print("Test 1: tridiagonal, inverse Lehmer matrix; rational elements"): 
print():
for n from nstart by nstep to nend do
x := Array(1..n): y := Array(1..n-1): z := Array(1..n-1):
delta, mu = 0, 0:

for i from 1 to n-1 do
x[i] := i/(1-1/(4*i*i)):
y[i] := -i*(i+1)/(2*i+1):
z[i] := y[i]:
end do:
x[n] := n*n/(2*n - 1):
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):
end do:

end proc:

# As Test1 but with non-zero delta and mu
Test2 := proc(invTest::procedure, descript::string,
nstart::integer, nend::integer, nstep::integer,
nRepeats::integer)
local n, x, y, z, delta, mu, i:

# output test routine header 
print(): 
print("Test 2: tridiagonal, inverse Lehmer matrix with", " non-zero mu and delta; rational elements"): 
print():

for n from nstart by nstep to nend do
x := Array(1..n): y := Array(1..n-1): z := Array(1..n-1):
delta, mu = 3*n/4, 2*n/3:

for i from 1 to n-1 do
x[i] := i/(1-1/(4*i*i)):
y[i] := -i*(i+1)/(2*i+1):
z[i] := y[i]:
end do:
x[n] := n*n/(2*n - 1):
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):
end do:

end proc:

# As Test1 but now make the matrix unsymmetric
Test3 := proc(invTest::procedure, descript::string,
          nstart::integer, nend::integer, nstep::integer,
          nRepeats::integer)
local n, x, y, z, delta, mu, i:

# output test routine header 
print(): 
print("Test 3: periodic non-symmetric tridiagonal with", " non-zero mu and delta; rational elements"): 
print():

for n from nstart by nstep to nend do
x := Array(1..n): y := Array(1..n-1): z := Array(1..n-1):
delta, mu = 0, 0:

for i from 1 to n-1 do
x[i] := i/(1-1/(4*i*i)):
y[i] := -i*(i+1)/(2*i+1):
z[i] := -y[i]/2:
end do:
x[n] := n*n/(2*n - 1):
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):
end do:

end proc:

# Algebraic test matrix -- generalization of the example used
# by El-Shehawey
Test4 := proc(invTest::procedure, descript::string,
     nstart::integer, nend::integer, nstep::integer,
     nRepeats::integer)
local n, x, y, z, delta, mu, i:

# output test routine header 
print(): 
print("Test 4: Generalized version of test example used", " by El-Shehawey; algebraic elements"): 
print():

for n from nstart by nstep to nend do
x := Array(1..n, fill=a):
y := Array(1..n-1, fill=1):
z := Array(1..n-1, fill=-1):
delta, mu = -1, 1:

testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):
end do:

end proc:
(*
x := Array(1..4,[1,2,3,4]):
y := Array(1..3,[1/5, 1/6, 1/7]):
z := Array(1..3,[1/2, 1/3, 1/4]):
delta:= 1/8:
mu := 1/9:
nRepeats := 10:
MTimes := Array(1..nRepeats):
OTimes := Array(1..nRepeats):
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):

x := Array(1..4,[a,a,a,a]):
y := Array(1..3,[1,1,1]):
z := Array(1..3,[-1,-1,-1]):
delta := -1:
mu := 1:
nRepeats := 10:
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):

for n from 10 by 10 to 50 do
x := Array(1..n, fill=a):
y := Array(1..n-1, fill=1):
z := Array(1..n-1, fill=-1):
delta := -1:
mu := 1:
nRepeats := 10:
testInv(x,y,z,delta,mu,nRepeats,KilicInv,"KilicInv"):
end do:

nRepeats := 10:
Test1(KilicInv, "KilicInv", 10, 30, 10, nRepeats);
Test2(KilicInv, "KilicInv", 10, 30, 10, nRepeats);
Test3(KilicInv, "KilicInv", 10, 30, 10, nRepeats); 
*)
#nRepeats := 10:
#Test2(KilicInv, "KilicInv", 10, 50, 10, nRepeats);
nRepeats := 1:
Test4(KilicInv, "KilicInv", 20, 20, 10, nRepeats);