> restart:
> # ITERATIVE REFINEMENT ALGORITHM 7.4
> #
> # To approximate the solution to the linear system Ax=b when A is
> # suspected to be ill-conditioned:
> #
> # INPUT:  The number of equations and unknowns n: the entries
> #         A(i,j), 1<=i, j<=n, of the matrix A: the entries b(i),
> #         1<=i<=n, of the inhomogeneous term b: the maximum number
> #         of iterations N.
> #
> # OUTPUT: The approximation XX(1),...,XX(n) or a message that the
> #         number of iterations was exceeded. 
> CHIP := proc(RND,R,X) local x,e:
> if RND = 1 then
> Digits := R:
> x := evalf(X):
> RETURN(x):
> else
> if X = 0 then 
> x := 0:
> else
> e := trunc(evalf(log10(abs(X)))):
> if abs(X) > 1 then
> e := e+1:
> fi:
> x := evalf(trunc(X*10^(R-e))*10^(e-R)):
> fi:
> RETURN(x):
> fi:
> end:
> print(`This is the Iterative Refinement Method.\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(`A(1,1), A(1,2), ..., A(1,n+1), A(2,1), A(2,2), ..., 
> A(2,n+1),..., A(n,1), A(n,2), ..., A(n,n+1)\n`):
> print(`Place as many entries as desired on each line, but separate\n`):
> print(`entries with `):
> print(`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):
> OK := FALSE:
> if AA = 1 then
> print(`Input the file name in the form - drive:\\name.ext\n`):
> print(`for example:   A:\\DATA.DTA\n`):
> NAME := scanf(`%s`)[1]:
> INP := fopen(NAME,READ,TEXT):
> OK := FALSE:
> while OK = FALSE do
> print(`Input the number of equations - an integer.\n`):
> N := scanf(`%d`)[1]:print(`N = `):print(N):
> if N > 0 then
> for I1 from 1 to N do
> for J from 1 to N+1 do
> A[I1-1,J-1] := fscanf(INP, `%f`)[1]:
> od:
> 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
> for I1 from 1 to N do
> for J from 1 to N+1 do
> print(`input entry in position `,I1,J): 
> A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]):
> od:
> od:
> OK := TRUE:
> else print(`The number must be a positive integer.\n`):
> fi:
> od:
> fi:
> if OK = TRUE then
> OK := FALSE:
> while OK = FALSE do
> print(`Input maximum number of iterations.\n`):
> # NN is used for the maximum number of iterations
> NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):
> if NN > 0 then
> OK := TRUE:
> else
> print(`Number must be a positive integer.\n`):
> fi:
> od:
> OK := FALSE:
> print(`Choice of rounding or chopping:\n`):
> print(`1. Rounding\n`):
> print(`2. Chopping\n`):
> print(`Enter 1 or 2.\n`):
> RND := scanf(`%d`)[1]:print(`Your input is `):print(RND):
> while OK = FALSE do
> print(`Input number of digits D <= 8 of rounding\n`):
> DD := scanf(`%d`)[1]:print(`Number of digits = `):print(DD):
> if DD > 0 then
> OK := TRUE:
> else
> print(`D must be a positive integer.\n`):
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input tolerance, which is usually 10^(-D).\n`):
> TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):
> if TOL > 0 then
> OK := TRUE:
> else
> print(`Tolerance must be a positive.\n`):
> fi:
> od:
> if OK = TRUE then
> 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, `ITERATIVE REFINEMENT METHOD\n\n`):
> M := N+1:
> fprintf(OUP, `Original system\n`):
> for I1 from 1 to N do
> for J from 1 to M do
> fprintf(OUP,` %.10e`,A[I1-1,J-1]):
> od:
> fprintf(OUP,`\n`):
> od:
> if RND = 1 then
> fprintf(OUP,`Rounding to %d Digits.\n`,DD):
> else fprintf(OUP,`Chopping to %d Digits.\n`,DD):
> fi:
> fprintf(OUP,`\n Modified System \n`):
> for I1 from 1 to N do
> NROW[I1-1] := I1:
> for J from 1 to M do
> A[I1-1,J-1] := CHIP(RND,DD,A[I1-1,J-1]):
> B[I1-1,J-1] := A[I1-1,J-1]:
> fprintf(OUP,`  %.10e`, A[I1-1,J-1]):
> od:
> fprintf(OUP, `\n`):
> od:
> # NROW and B have been initialized, Gauss elimination will begin.
> # Step 0
> I1 := 1:
> while I1 <= N-1 and OK = TRUE do
> KK := I1:
> while abs(A[KK-1,I1-1]) < 1.0e-20 and KK <= N do
> KK := KK+1:
> od:
> if KK > N then
> OK := false:
> fprintf(OUP, `System does not have a unique solution.\n`):
> else 
> if KK <> I1 then
> # Row interchange is necessary
> IS := NROW[I1-1]:
> NROW[I1-1] := NROW[KK-1]:
> NROW[KK-1] := IS:
> for J from 1 to M do
> C := A[I1-1,J-1]:
> A[I1-1,J-1] := A[KK-1,J-1]:
> A[KK-1,J-1] := C:
> od:
> fi:
> for J from I1+1 to N do
> A[J-1,I1-1] := CHIP(RND,DD,A[J-1,I1-1]/A[I1-1,I1-1]):
> for L from I1+1 to M do
> A[J-1,L-1] := CHIP(RND,DD,A[J-1,L-1]-CHIP(RND,DD,A[J-1,I1-1]*A[I1-1,L-1])):
> od:
> od:
> fi:
> I1 := I1+1:
> od:
> if abs(A[N-1,N-1]) < 1.0e-20 and OK = TRUE then
> OK := FALSE:
> fprintf(OUP, `System has singular matrix\n`):
> fi:
> if OK = TRUE then
> fprintf(OUP, `Reduced system\n`):
> for I1 from 1 to N do
> for J from 1 to M do
> fprintf(OUP, `  %.10e`, A[I1-1,J-1]):
> od:
> fprintf(OUP, `\n`):
> od:
> X[N-1] := CHIP(RND,DD,A[N-1,M-1]/A[N-1,N-1]):
> for I1 from 1 to N-1 do
> J := N-I1:
> S := 0.0:
> for L from J+1 to N do
> S := CHIP(RND,DD,S-CHIP(RND,DD,A[J-1,L-1]*X[L-1])):
> od:
> S := CHIP(RND,DD,A[J-1,M-1]+S):
> X[J-1] := CHIP(RND,DD,S/A[J-1,J-1]):
> od:
> fi:
> fprintf(OUP, `Initial solution\n`):
> for I1 from 1 to N do
> fprintf(OUP,`  %.10e`, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> # Refinement begins
> # Step 1
> if OK = TRUE then
> K := 1:
> for I1 from 1 to N do
> XX[I1-1] := X[I1-1]:
> od:
> fi:
> # Step 2
> while OK = TRUE and K <= NN do
> # LL is set to 1 if the desired accuracy in any component is not
> # achieved.
> LL := 0:
> # Step 3
> for I1 from 1 to N do
> R[I1-1] := 0:
> for J from 1 to N do
> R[I1-1] := CHIP(RND,2*DD,R[I1-1]-CHIP(RND,2*DD,B[I1-1,J-1]*XX[J-1])):
> od:
> R[I1-1] := CHIP(RND,2*DD,B[I1-1,M-1]+R[I1-1]):
> od:
> fprintf(OUP, `Residual number %d\n`, K):
> for I1 from 1 to N do
> R[I1-1] := CHIP(RND,DD,R[I1-1]):
> fprintf(OUP, `%18.10e `, R[I1-1]):
> od:
> fprintf(OUP, `\n`):
> # Step 4
> # Solve the linear system in the same order as in Step 0.
> for I1 from 1 to N-1 do
> I2 := NROW[I1-1]:
> for J from I1+1 to N do
> J1 := NROW[J-1]:
> R[J1-1] := CHIP(RND,DD,R[J1-1]-CHIP(RND,DD,A[J-1,I1-1]*R[I2-1])):
> od:
> od:
> X[N-1] := CHIP(RND,DD,R[NROW[N-1]-1]/A[N-1,N-1]):
> for I1 from 1 to N-1 do
> J := N-I1:
> S := 0:
> for L from J+1 to N do
> S := CHIP(RND,DD,S-CHIP(RND,DD,A[J-1,L-1]*X[L-1])):
> od:
> S := CHIP(RND,DD,S+R[NROW[J-1]-1]):
> X[J-1] := CHIP(RND,DD,S/A[J-1,J-1]):
> od:
> fprintf(OUP, `Vector Y\n`):
> for I1 from 1 to N do
> fprintf(OUP,`%18.10e `, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> # Steps 5 and 6
> XXMAX := 0:
> YMAX := 0:
> ERR1 := 0:
> for I1 from 1 to N do
> # If not sufficiently accurate, set LL to 1.
> if abs(X[I1-1]) > TOL then
> LL := 1:
> fi:
> if K = 1 then
> if abs(X[I1-1]) > YMAX then
> YMAX := abs(X[I1-1]):
> fi:
> if abs(XX[I1-1]) > XXMAX then
> XXMAX := abs(XX[I1-1]):
> fi:
> fi:
> TEMP := XX[I1-1]:
> XX[I1-1] := CHIP(RND,DD,XX[I1-1]+X[I1-1]):
> TEMP := abs(TEMP-XX[I1-1]):
> if TEMP > ERR1 then
> ERR1 := TEMP:
> fi:
> od:
> if ERR1 <= TOL then
> LL := 2:
> fi:
> if K = 1 then
> COND := YMAX/XXMAX*10^DD:
> fi:
> fprintf(OUP, `New approximation\n`):
> for I1 from 1 to N do
> fprintf(OUP, `%18.10e `, XX[I1-1]):
> od:
> fprintf(OUP, `\n`):
> # Step 7
> if LL = 0 then
> fprintf(OUP, `The above vector is the solution.\n`):
> OK := FALSE:
> else
> if LL = 2 then
> fprintf(OUP,`The above vector is the best possible\n`):
> fprintf(OUP,`with TOL := %18.10e \n`,TOL):
> OK := FALSE:
> else
> K := K+1:
> fi
> fi:
> # Step 8 is not used in this implementation.
> od:
> if K > NN then
> print( `Maximum Number of Iterations Exceeded.`):
> fi:
> fprintf(OUP, `Condition number is %.10e\n`, COND):
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi:
> fi: