alg074.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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. 
> alg074 := proc() local CHIP, A1, OK, NAME, INP, N, I, J, A, NN, RND, D, TOL, FLAG, OUP, M, NROW, B, KK, IS, C, L, X, S, K, XX, LL, R, I1, J1, XXMAX, YMAX, ERR1, TEMP, COND;
> 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;
> printf(`This is the Iterative Refinement Method.\n`);
> printf(`The array will be input from a text file in the order\n`);
> printf(`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`);
> printf(`Place as many entries as desired on each line, but separate\n`);
> printf(`entries with `);
> printf(`at least one blank.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> A1 := scanf(`%c`)[1];
> OK := FALSE;
> if A1 = "Y" or A1 = "y" then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number of equations - an integer.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> for I from 1 to N do
> for J from 1 to N+1 do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> od;
> od;
> OK := TRUE;
> fclose(INP);
> else
> printf(`The number must be a positive integer\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations.\n`);
> # NN is used for the maximum number of iterations
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Number must be a positive integer.\n`);
> fi;
> od;
> OK := FALSE;
> printf(`Choice of rounding or chopping:\n`);
> printf(`1. Rounding\n`);
> printf(`2. Chopping\n`);
> printf(`Enter 1 or 2.\n`);
> RND := scanf(`%d`)[1];
> while OK = FALSE do
> printf(`Input number of digits D <= 8 of rounding\n`);
> D := scanf(`%d`)[1];
> if D > 0 then
> OK := TRUE;
> else
> printf(`D must be a positive integer.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input tolerance, which is usually 10^(-D).\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be a positive.\n`);
> fi;
> od;
> else
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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 I from 1 to N do
> for J from 1 to M do
> fprintf(OUP,` %.10e`,A[I-1,J-1]);
> od;
> fprintf(OUP,`\n`);
> od;
> if RND = 1 then
> fprintf(OUP,`Rounding to %d Digits.\n`,D);
> else fprintf(OUP,`Chopping to %d Digits.\n`,D);
> fi;
> fprintf(OUP,`\n Modified System \n`);
> for I from 1 to N do
> NROW[I-1] := I;
> for J from 1 to M do
> A[I-1,J-1] := CHIP(RND,D,A[I-1,J-1]);
> B[I-1,J-1] := A[I-1,J-1];
> fprintf(OUP,`  %.10e`, A[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> # NROW and B have been initialized, Gauss elimination will begin.
> # Step 0
> I := 1;
> while I <= N-1 and OK = TRUE do
> KK := I;
> while abs(A[KK-1,I-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 <> I then
> # Row interchange is necessary
> IS := NROW[I-1];
> NROW[I-1] := NROW[KK-1];
> NROW[KK-1] := IS;
> for J from 1 to M do
> C := A[I-1,J-1];
> A[I-1,J-1] := A[KK-1,J-1];
> A[KK-1,J-1] := C;
> od;
> fi;
> for J from I+1 to N do
> A[J-1,I-1] := CHIP(RND,D,A[J-1,I-1]/A[I-1,I-1]);
> for L from I+1 to M do
> A[J-1,L-1] := CHIP(RND,D,A[J-1,L-1]-CHIP(RND,D,A[J-1,I-1]*A[I-1,L-1]));
> od;
> od;
> fi;
> I := I+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 I from 1 to N do
> for J from 1 to M do
> fprintf(OUP, `  %.10e`, A[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> X[N-1] := CHIP(RND,D,A[N-1,M-1]/A[N-1,N-1]);
> for I from 1 to N-1 do
> J := N-I;
> S := 0.0;
> for L from J+1 to N do
> S := CHIP(RND,D,S-CHIP(RND,D,A[J-1,L-1]*X[L-1]));
> od;
> S := CHIP(RND,D,A[J-1,M-1]+S);
> X[J-1] := CHIP(RND,D,S/A[J-1,J-1]);
> od;
> fi;
> fprintf(OUP, `Initial solution\n`);
> for I from 1 to N do
> fprintf(OUP,`  %.10e`, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> # Refinement begins
> # Step 1
> if OK = TRUE then
> K := 1;
> for I from 1 to N do
> XX[I-1] := X[I-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 I from 1 to N do
> R[I-1] := 0;
> for J from 1 to N do
> R[I-1] := CHIP(RND,2*D,R[I-1]-CHIP(RND,2*D,B[I-1,J-1]*XX[J-1]));
> od;
> R[I-1] := CHIP(RND,2*D,B[I-1,M-1]+R[I-1]);
> od;
> fprintf(OUP, `Residual number %d\n`, K);
> for I from 1 to N do
> R[I-1] := CHIP(RND,D,R[I-1]);
> fprintf(OUP, `%18.10e `, R[I-1]);
> od;
> fprintf(OUP, `\n`);
> # Step 4
> # Solve the linear system in the same order as in Step 0.
> for I from 1 to N-1 do
> I1 := NROW[I-1];
> for J from I+1 to N do
> J1 := NROW[J-1];
> R[J1-1] := CHIP(RND,D,R[J1-1]-CHIP(RND,D,A[J-1,I-1]*R[I1-1]));
> od;
> od;
> X[N-1] := CHIP(RND,D,R[NROW[N-1]-1]/A[N-1,N-1]);
> for I from 1 to N-1 do
> J := N-I;
> S := 0;
> for L from J+1 to N do
> S := CHIP(RND,D,S-CHIP(RND,D,A[J-1,L-1]*X[L-1]));
> od;
> S := CHIP(RND,D,S+R[NROW[J-1]-1]);
> X[J-1] := CHIP(RND,D,S/A[J-1,J-1]);
> od;
> fprintf(OUP, `Vector Y\n`);
> for I from 1 to N do
> fprintf(OUP,`%18.10e `, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> # Steps 5 and 6
> XXMAX := 0;
> YMAX := 0;
> ERR1 := 0;
> for I from 1 to N do
> # If not sufficiently accurate, set LL to 1.
> if abs(X[I-1]) > TOL then
> LL := 1;
> fi;
> if K = 1 then
> if abs(X[I-1]) > YMAX then
> YMAX := abs(X[I-1]);
> fi;
> if abs(XX[I-1]) > XXMAX then
> XXMAX := abs(XX[I-1]);
> fi;
> fi;
> TEMP := XX[I-1];
> XX[I-1] := CHIP(RND,D,XX[I-1]+X[I-1]);
> TEMP := abs(TEMP-XX[I-1]);
> if TEMP > ERR1 then
> ERR1 := TEMP;
> fi;
> od;
> if ERR1 <= TOL then
> LL := 2;
> fi;
> if K = 1 then
> COND := YMAX/XXMAX*10^D;
> fi;
> fprintf(OUP, `New approximation\n`);
> for I from 1 to N do
> fprintf(OUP, `%18.10e `, XX[I-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
> fprintf(OUP, `Maximum Number of Iterations Exceeded.\n`);
> fi;
> fprintf(OUP, `Condition number is %.10e\n`, COND);
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg074();