> restart:
> # CROUT FACTORIZATION FOR TRIDIAGONAL LINEAR SYSTEMS ALGORITHM 6.7
> #
> # To solve the n x n linear system
> #
> # E1:  A[1,1] X[1] + A[1,2] X[2]                  = A[1,n+1]
> # E2:  A[2,1] X[1] + A[2,2] X[2] + A[2,3] X[3]    = A[2,n+1]
> # :
> # .
> # E(n):          A[n,n-1] X[n-1] + A[n,n] X[n]    = A[n,n+1]
> #
> # INPUT:   the dimension n: the entries of A.
> #
> # OUTPUT:  the solution X(1), ..., X(N).
> print(`This is Crout Method for tridiagonal linear systems.\n`):
> print(`Choice of input method`):
> print(`1. input from keyboard - not recommended for large systems`):
> print(`2. input from a text file`):
> print(`Please enter 1 or 2.`):
> FLAG := scanf(`%d`)[1]: print(`Your input is`): print(FLAG):
> if FLAG = 2 then
> print(`The array will be input from a text file in the order:\n`):
> print(`diagonal entries, lower sub-diagonal entries,\n`):
> print(`upper sub-diagonal entries, inhomogeneous term.\n\n`):
> print(`Place as many entries as desired on each line, but separate `):
> print(`entries \n`):
> print(`with at least one blank.\n\n\n`):
> print(`Has the input file been created? - enter 1 for yes or 2 for no.`):
> AA := scanf(`%d`)[1]: print(`Your response is`): print(AA):
> if AA = 1 then
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example:   A:\\DATA.DTA`):
> NAME := scanf(`%s`)[1]: print(`The file name is`): print(NAME):
> INP := fopen(NAME,READ,TEXT):
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number of equations - an integer.`):
> N := scanf(`%d`)[1]: print(`N is`): print(N):
> if N > 0 then
> # A(I,I) is stored in A(I), 1 <= I <= n 
> for I1 from 1 to N do
> A[I1-1] := fscanf(INP, `%f`)[1]:
> od:
> # the lower sub-diagonal A(I,I-1) is stored
> #in B(I), 2 <= I <= n 
> for I1 from 2 to N do
> B[I1-1] := fscanf(INP, `%f`)[1]:
> od:
> # the upper sub-diagonal A(I,I+1) is stored
> #in C(I), 1 <= I <= n-1 
> NN := N-1:
> for I1 from 1 to NN do
> C[I1-1] := fscanf(INP, `%f`)[1]:
> od:
> # A(I,N+1) is stored in BB(I), 1 <= I <= n 
> for I1 from 1 to N do
> BB[I1-1] := fscanf(INP, `%f`)[1]:
> od:
> OK := TRUE:
> fclose(INP):
> else print(`The number must be a positive integer.\n`):
> fi:
> od:
> else 
> print(`The program will end so the input file can be created.\n`):
> fi:
> else
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number of equations - an integer.`):
> N := scanf(`%d`)[1]: print(`N= `): print(N):
> if N > 0 then
> # A(I,I) is stored in A(I), 1 <= I <= n 
> for I1 from 1 to N do
> print(`Input the diagonal entry in row = `):print(I1):
> A[I1-1] := scanf( `%f`)[1]:print(`Entry = `,A[I1-1]):
> od:
> # the lower sub-diagonal A(I,I-1) is stored
> #in B(I), 2 <= I <= n 
> for I1 from 2 to N do
> print(`Input the lower sub-diagonal diagonal entry in row = `):print(I1):print(` and column = `):print(I1-1):
> B[I1-1] := scanf(`%f`)[1]:print(`Entry = `,B[I1-1]):
> od:
> # the upper sub-diagonal A(I,I+1) is stored
> #in C(I), 1 <= I <= n-1 
> NN := N-1:
> for I1 from 1 to NN do
> print(`Input the upper sub-diagonal diagonal entry in row = `):print(I1):print(` and column = `):print(I1+1):
> C[I1-1] := scanf(`%f`)[1]:print(`Entry = `,C[I1-1]):
> od:
> # A(I,N+1) is stored in BB(I), 1 <= I <= n 
> for I1 from 1 to N do
> print(`Input the entry on the right-hand side of equation = `):print(I1):
> BB[I1-1] := scanf(`%f`)[1]:print(`Entry = `,BB[I1-1]):
> od:
> OK:=TRUE:
> else print(`The number must be a positive integer.\n`):
> fi:
> od:
> fi:
> if OK = TRUE then
> OUP := default:
> fprintf(OUP, `The original system \n`):
> fprintf(OUP,`Diagonal Entries`):
> for I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, A[I1-1]):
> od:
> fprintf(OUP, `\n`):
> fprintf(OUP,`Lower Sub-Diagonal Entries`):
> for I1 from 2 to N do
> fprintf(OUP, ` %11.8f`, B[I1-1]):
> od:
> fprintf(OUP, `\n`):
> fprintf(OUP,`Upper Sub-Diagonal Entries`):
> for I1 from 1 to NN do
> fprintf(OUP, ` %11.8f`, C[I1-1]):
> od:
> fprintf(OUP, `\n`):
> fprintf(OUP,`Inhomogeneous Entries`):
> for I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, BB[I1-1]):
> od:
> fprintf(OUP, `\n`):
> if OK = TRUE then
> # Steps 1 - 3 set up and solve Lz = 0
> # Step 1
> # The entries of U overwrite C and the entries of L overwrite A
> C[0] := C[0]/A[0]:
> Z[0] := BB[0]/A[0]:
> # Step 2
> for I1 from 2 to NN do
> A[I1-1] := A[I1-1]-B[I1-1]*C[I1-2]:
> C[I1-1] := C[I1-1]/A[I1-1]:
> Z[I1-1] := (BB[I1-1]-B[I1-1]*Z[I1-2])/A[I1-1]:
> od:
> # Step 3
> A[N-1] := A[N-1]-B[N-1]*C[N-2]:
> Z[N-1] := (BB[N-1]-B[N-1]*Z[N-2])/A[N-1]:
> # Step 4
> # Steps 4 and 5 solve Ux = z
> X[N-1] := Z[N-1]:
> # Step 5
> for I2 from 1 to NN do
> I1 := NN-I2+1:
> X[I1-1] := Z[I1-1]-C[I1-1]*X[I1]:
> od:
> # Step 6
> print(`Choice of output method:\n`):
> print(`1. Output to screen\n`):
> print(`2. Output to text file\n`):
> print(`Please enter 1 or 2.\n`):
> FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive:\\name.ext\n`):
> print(`for example:   A:\\OUTPUT.DTA\n`):
> NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `CROUT METHOD FOR TRIDIAGONAL LINEAR SYSTEMS\n\n`):
> fprintf(OUP, `The solution is\n`):
> for I1 from 1 to N do
> fprintf(OUP, `  %12.8f`, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi:
> fi: