rem system of linear equations solver
rem compact crout

print "Linear Equations - Compact Crout"
print "to solve multiple symmetric equations of the form"
print "A1*x1+A2*x2+...+An*xn=B"
print "data is input by columns using the upper triangle of matrix A"
print "enter 'q' at any prompt to quit"

dim a(120)
dim b(15)
e5=.00000001

gosub [inputA]
gosub [factorA]
print
print "Determinant of A"
print
print "Det = ";d8
print
[again]
rem read vector B
gosub [inputB]
rem solution of A*X=B
gosub [getX]
rem output vector X
gosub [outputX]
print
print "Do you want to input another vector B (y/n)? ";
input " ";w$
if w$="n" or w$="q" then end
goto [again]

[getX]
rem forward solution
j0=1
rem loop over j columns
for j=2 to n
 j2=j-1
 rem reduce component b(j)
 for k=1 to j2
  j4=j0+k
  b(j)=b(j)-a(j4)*b(k)
 next k
 j0=j0+k
next j
rem back substitution
j0=0
rem perform b(j)/a(j,j)
for j=1 to n
 j0=j0+j
 b(j)=b(j)/a(j0)
next j
rem backward loop over j columns
j=n
j0=n*(n+1)/2
[djump]
if (j<=1) goto [cjump]
j2=j-1
rem subtract columns
j0=j0-j
for i=1 to j2
 i1=j0+i
 b(i)=b(i)-a(i1)*b(j)
next i
j=j2
goto [djump]
[cjump]
return

[factorA]
j0=1
rem check for zero pivot
d8=a(j0)
if abs(d8)<e5 goto [fatal]
rem loop over j columns
for j=2 to n
 j2=j-1
 if j=2 goto [bjump]
 i0=1
 rem loop over i rows
 for i=2 to j2
  i2=i-1
  i1=j0+i
  rem reduce element a(i,j)
  for k=1 to i2
   k1=i0+k
   k2=j0+k
   a(i1)=a(i1)-1*a(k1)*a(k2)
  next k
  i0=i0+i
 next i
 rem reduce column j
 rem reduce element (aj,j)
 [bjump]
 j3=j0+j
 i3=0
 d1=a(j3)
 for i=1 to j2
  i3=i3+i
  i1=j0+i
  u=a(i1)/a(i3)
  d1=d1-1*u*a(i1)
  a(i1)=u
 next i
 a(j3)=d1
 rem check for zero pivot
 if abs(d1)<e5 goto [fatal]
 j0=j0+j
 d8=d8*d1
next j
return

[inputA]
print
input "Enter A matrix dimension: ";n$
if n$="q" then end
n=val(n$)
j0=0
print
print "Enter A matrix upper triangle"
for j=1 to n
 [ajump]
 print
 print "upper part of column ";j
 print
 for i=1 to j
  i1=i+j0
  print "A(";i;",";j;")= ";
  input " ";a$
  if a$="q" then end
  a(i1)=val(a$)
 next i
 print
 input "data correct (y/n)? ";w$
 if w$="q" then end
 if w$="n" goto [ajump]
 j0=j0+j
next j
print
return

[inputB]
print
print "input vector B (constants) "
print
for i=1 to n
 print "B(";i;")= ";
 input " ";b$
 if b$="q" then end
 b(i)=val(b$)
next i
print
return

[outputX]
print
print "Solution"
print
for i=1 to n
 print "X(";i;")= ";b(i)
next i
return

[fatal]
print
print "Sorry, can't proceed..."
print "matrix A is not positive-definite"
print
end