alg075.txt

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

> restart;
> # CONJUGATE GRADIENT ALGORITHM 7.5
> #
> # To solve Ax = b given the preconditioning matrix C inverse
> # and an initial approximation
> # x(0):
> #
> # 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
> #          entries C(I,J), 1<=I, J<=n, of the preconditioning
> #          matrix C inverse, entries XO(I), 1<=I<=n, of x(0);
> #
> # OUTPUT:  the approximate solution X(1),...,X(n) and its
> #          residual vector R(1),...,R(N) or a message
> #          that the number of iterations was exceeded.
> alg075 := proc() local OK,AA,NAME,INP,N,I,J,A,X1,TOL,NN,W,K,ERR,S,FLAG,OUP,R,T,ALPHA,BETA,U,V,CI,QERR,ERR1,CT,SS,Z;
> printf(`This is the Conjugate Gradient Method for Linear Systems.\n`);
> OK := FALSE;
> 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),\n`);
> printf(`..., A(n,1), A(n,2), ..., A(n,n+1)\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with at least one blank.\n`);
> printf(`Do the same for the input of the inverse of C.\n`);
> printf(`The initial approximation should follow in same format.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> AA := scanf(`%c`)[1];
> if AA = "Y" or AA = "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;
> for I from 1 to N do
> for J from 1 to N do
> CI[I-1,J-1] := fscanf(INP, `%f`)[1];
> CT[J-1,I-1] := CI[I-1,J-1];
> od;
> od;
> for I from 1 to N do
> X1[I-1] := fscanf(INP, `%f`)[1];
> 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 the tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be a positive number.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Number must be a positive integer.\n`);
> fi;
> od;
> else
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> # Step 1
> for I from 1 to N do
> R[I-1] := A[I-1,N];
> for J from 1 to N do
> R[I-1] := R[I-1]-A[I-1,J-1]*X1[J-1];
> od;
> od;
> for I from 1 to N do
> W[I-1] := 0;
> for J from 1 to N do
> W[I-1] := W[I-1]+CI[I-1,J-1]*R[J-1];
> od;
> od;
> for I from 1 to N do
> V[I-1] := 0;
> for J from 1 to N do
> V[I-1] := V[I-1]+CT[I-1,J-1]*W[J-1];
> od;
> od;
> ALPHA := 0.0;
> for I from 1 to N do
> ALPHA := ALPHA + W[I-1]*W[I-1];
> od;
> # Step 2
> K := 1;
> OK := FALSE;
> # Step 3
> while (OK = FALSE) and (K <= NN) do
> ERR := 0;
> for I from 1 to N do
> ERR := ERR + V[I-1]*V[I-1];
> od;
> # Step 4
> if sqrt(ERR) < TOL then
> K := K -1;
> OK := TRUE:
> else
> # Step 5
> for I from 1 to N do
> U[I-1] := 0.0;
> for J from 1 to N do
> U[I-1] := U[I-1]+A[I-1,J-1]*V[J-1];
> od;
> od;
> SS := 0.0;
> for I from 1 to N do
> SS := SS + V[I-1]*U[I-1];
> od;
> T := ALPHA/SS:
> for I from 1 to N do
> X1[I-1] := X1[I-1]+T*V[I-1];
> R[I-1] := R[I-1] - T*U[I-1];
> od;
> for I from 1 to N do
> W[I-1] := 0.0;
> for J from 1 to N do
> W[I-1] := W[I-1]+CI[I-1,J-1]*R[J-1];
> od;
> od;
> BETA := 0.0;
> for I from 1 to N do
> BETA := BETA + W[I-1]*W[I-1];
> od;
> ERR1 := sqrt(BETA);
> # Step 6
> if ERR1 <= TOL then
> ERR := 0.0;
> for I from 1 to N do
> ERR := ERR + R[I-1]*R[I-1];
> od;
> ERR := sqrt(ERR);
> if ERR < TOL then
> OK := TRUE;
> fi;
> fi;
> if OK = FALSE then
> # Step 7
> K := K + 1;
> S := BETA/ALPHA;
> for I from 1 to N do
> Z[I-1] := 0;
> for J from 1 to N do
> Z[I-1] := Z[I-1]+CT[I-1,J-1]*W[J-1];
> od;
> od;
> for I from 1 to N do
> V[I-1] := Z[I-1]+S*V[I-1];
> od;
> ALPHA := BETA;
> fi;
> fi;
> od;
> # Step 8
> if OK = FALSE then
> printf(`Maximum Number of Iterations Exceeded.\n`);
> else
> 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, `PRECONDITIONED CONJUGATE GRADIENT METHOD\n\n`);
> fprintf(OUP, `The solution vector is :\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X1[I-1]);
> od;
> fprintf(OUP, `\nusing %d iterations with\n`, K);
> fprintf(OUP, `Tolerance  %.10e in infinity-norm\n`, TOL);
> fprintf(OUP, `The residual vector is :\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, R[I-1]);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> fi;
> RETURN(0);
> end;
> alg075();