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.\134n`):

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\134n`):

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)\134n`):

print(`Place as many entries as desired on each line, but separate\134n`):

print(`entries with `):

print(`at least one blank.\134n\134n\134n`):

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:\134\134name.ext\134n`):

print(`for example: A:\134\134DATA.DTA\134n`):

NAME := scanf(`%s`)[1]:

INP := fopen(NAME,READ,TEXT):

OK := FALSE:

while OK = FALSE do

print(`Input the number of equations - an integer.\134n`):

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:

OK := TRUE:

fclose(INP):

else

print(`The number must be a positive integer\134n`):

fi:

od:

else

print(`The program will end so the input file can be created.\134n`):

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:

OK := TRUE:

else print(`The number must be a positive integer.\134n`):

fi:

od:

fi:

if OK = TRUE then

OK := FALSE:

while OK = FALSE do

print(`Input maximum number of iterations.\134n`):

# 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.\134n`):

fi:

od:

OK := FALSE:

print(`Choice of rounding or chopping:\134n`):

print(`1. Rounding\134n`):

print(`2. Chopping\134n`):

print(`Enter 1 or 2.\134n`):

RND := scanf(`%d`)[1]:print(`Your input is `):print(RND):

while OK = FALSE do

print(`Input number of digits D <= 8 of rounding\134n`):

DD := scanf(`%d`)[1]:print(`Number of digits = `):print(DD):

if DD > 0 then

OK := TRUE:

else

print(`D must be a positive integer.\134n`):

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input tolerance, which is usually 10^(-D).\134n`):

TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):

if TOL > 0 then

OK := TRUE:

else

print(`Tolerance must be a positive.\134n`):

fi:

od:

if OK = TRUE then

print(`Choice of output method:\134n`):

print(`1. Output to screen\134n`):

print(`2. Output to text file\134n`):

print(`Please enter 1 or 2.\134n`):

FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext\134n`):

print(`for example: A:\134\134OUTPUT.DTA\134n`):

NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):

OUP := fopen(NAME,WRITE,TEXT):

else

OUP := default:

fi:

fprintf(OUP, `ITERATIVE REFINEMENT METHOD\134n\134n`):

M := N+1:

fprintf(OUP, `Original system\134n`):

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,`\134n`):

od:

if RND = 1 then

fprintf(OUP,`Rounding to %d Digits.\134n`,DD):

else fprintf(OUP,`Chopping to %d Digits.\134n`,DD):

fi:

fprintf(OUP,`\134n Modified System \134n`):

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, `\134n`):

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.\134n`):

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:

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\134n`):

fi:

if OK = TRUE then

fprintf(OUP, `Reduced system\134n`):

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, `\134n`):

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\134n`):

for I1 from 1 to N do

fprintf(OUP,` %.10e`, X[I1-1]):

od:

fprintf(OUP, `\134n`):

# 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\134n`, 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, `\134n`):

# 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:

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\134n`):

for I1 from 1 to N do

fprintf(OUP,`%18.10e `, X[I1-1]):

od:

fprintf(OUP, `\134n`):

# 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:

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\134n`):

for I1 from 1 to N do

fprintf(OUP, `%18.10e `, XX[I1-1]):

od:

fprintf(OUP, `\134n`):

# Step 7

if LL = 0 then

fprintf(OUP, `The above vector is the solution.\134n`):

OK := FALSE:

else

if LL = 2 then

fprintf(OUP,`The above vector is the best possible\134n`):

fprintf(OUP,`with TOL := %18.10e \134n`,TOL):

OK := FALSE:

else

K := K+1:

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\134n`, COND):

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created successfully`):

fi:

SUpUaGlzfmlzfnRoZX5JdGVyYXRpdmV+UmVmaW5lbWVudH5NZXRob2QufCtHNiI=

STdDaG9pY2V+b2Z+aW5wdXR+bWV0aG9kRzYi

SWVuMS5+aW5wdXR+ZnJvbX5rZXlib2FyZH4tfm5vdH5yZWNvbW1lbmRlZH5mb3J+bGFyZ2V+c3lzdGVtc0c2Ig==

SToyLn5pbnB1dH5mcm9tfmF+dGV4dH5maWxlRzYi

STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig==

SS5Zb3VyfmlucHV0fmlzRzYi

IiIj

SVdUaGV+YXJyYXl+d2lsbH5iZX5pbnB1dH5mcm9tfmF+dGV4dH5maWxlfmlufnRoZX5vcmRlcnwrRzYi

SVxxQSgxLDEpLH5BKDEsMiksfi4uLix+QSgxLG4rMSksfkEoMiwxKSx+QSgyLDIpLH4uLi4sfnwrQSgyLG4rMSksLi4uLH5BKG4sMSksfkEobiwyKSx+Li4uLH5BKG4sbisxKXwrRzYi

SWduUGxhY2V+YXN+bWFueX5lbnRyaWVzfmFzfmRlc2lyZWR+b25+ZWFjaH5saW5lLH5idXR+c2VwYXJhdGV8K0c2Ig==

SS5lbnRyaWVzfndpdGh+RzYi

STdhdH5sZWFzdH5vbmV+YmxhbmsufCt8K3wrRzYi

SWpuSGFzfnRoZX5pbnB1dH5maWxlfmJlZW5+Y3JlYXRlZD9+LX5lbnRlcn4xfmZvcn55ZXN+b3J+Mn5mb3J+bm8uRzYi

STFZb3VyfnJlc3BvbnNlfmlzRzYi

IiIi

SVNJbnB1dH50aGV+ZmlsZX5uYW1lfmlufnRoZX5mb3Jtfi1+ZHJpdmU6XG5hbWUuZXh0fCtHNiI=

STxmb3J+ZXhhbXBsZTp+fn5BOlxEQVRBLkRUQXwrRzYi

SU1JbnB1dH50aGV+bnVtYmVyfm9mfmVxdWF0aW9uc34tfmFufmludGVnZXIufCtHNiI=

SSVOfj1+RzYi

IiIk

SUVJbnB1dH5tYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zLnwrRzYi

SUBNYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zfj1+RzYi

IiM6

SUFDaG9pY2V+b2Z+cm91bmRpbmd+b3J+Y2hvcHBpbmc6fCtHNiI=

SS0xLn5Sb3VuZGluZ3wrRzYi

SS0yLn5DaG9wcGluZ3wrRzYi

SS9FbnRlcn4xfm9yfjIufCtHNiI=

SS9Zb3VyfmlucHV0fmlzfkc2Ig==

IiIi

SUtJbnB1dH5udW1iZXJ+b2Z+ZGlnaXRzfkR+PD1+OH5vZn5yb3VuZGluZ3wrRzYi

STROdW1iZXJ+b2Z+ZGlnaXRzfj1+RzYi

IiIm

SUxJbnB1dH50b2xlcmFuY2UsfndoaWNofmlzfnVzdWFsbHl+MTBeKC1EKS58K0c2Ig==

SS1Ub2xlcmFuY2V+PX5HNiI=

JCIiIiEiJg==

STpDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDp8K0c2Ig==

STUxLn5PdXRwdXR+dG9+c2NyZWVufCtHNiI=

STgyLn5PdXRwdXR+dG9+dGV4dH5maWxlfCtHNiI=

STZQbGVhc2V+ZW50ZXJ+MX5vcn4yLnwrRzYi

SS9Zb3VyfmlucHV0fmlzfkc2Ig==

IiIi

ITERATIVE REFINEMENT METHOD Original system 3.3330000000e+00 1.5920000000e+04 -1.0333000000e+01 1.5913000000e+04 2.2220000000e+00 1.6710000000e+01 9.6120000000e+00 2.8544000000e+01 1.5611000000e+00 5.1791000000e+00 1.6852000000e+00 8.4254000000e+00 Rounding to 5 Digits. Modified System 3.3330000000e+00 1.5920000000e+04 -1.0333000000e+01 1.5913000000e+04 2.2220000000e+00 1.6710000000e+01 9.6120000000e+00 2.8544000000e+01 1.5611000000e+00 5.1791000000e+00 1.6852000000e+00 8.4254000000e+00 Reduced system 3.3330000000e+00 1.5920000000e+04 -1.0333000000e+01 1.5913000000e+04 6.6667000000e-01 -1.0596000000e+04 1.6501000000e+01 -1.0580000000e+04 4.6838000000e-01 7.0323000000e-01 -5.0790000000e+00 -4.7000000000e+00 Initial solution 1.2001000000e+00 9.9991000000e-01 9.2538000000e-01 Residual number 1 -5.1800000000e-03 2.7413000000e-01 -1.8616000000e-01 Vector Y -2.0008000000e-01 8.9989000000e-05 7.4607000000e-02 New approximation 1.0000000000e+00 1.0000000000e+00 9.9999000000e-01 Residual number 2 -1.0000000000e-04 9.6120000000e-05 1.6852000000e-05 Vector Y -1.5002000000e-09 2.0951000000e-10 1.0000000000e-05 New approximation 1.0000000000e+00 1.0000000000e+00 1.0000000000e+00 The above vector is the best possible with TOL := 1.0000000000e-05 Condition number is 1.6671944000e+04

JSFH