program ALG075; { 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, the entries XO(I), 1<=I<=n, of x(0); tolerance TOL; maximum number of iterations N. OUTPUT: the approximate solution X(1),...,X(N) or a message that the number of iterations was exceeded. } var INP,OUP : text; A : array [ 1..10, 1..11 ] of real; CI, CT : array [1..10,1..10] of real; R,V,W,X1,U,Z : array [ 1..10 ] of real; ALPHA,BETA,SS,S,ERR,TOL,ERR1,T,QERR : real; FLAG,N,I,J,NN,K : integer; OK : boolean; AA : char; NAME : string [ 30 ]; procedure INPUT; begin writeln('This is the Conjugate Gradient Method for Linear Systems.'); OK := false; writeln ('The array will be input from a text file in the order: '); writeln('A(1,1), A(1,2), ..., A(1,n+1), A(2,1), A(2,2), ..., A(2,n+1),'); write ('..., A(n,1), '); writeln ('A(n,2), ..., A(n,n+1) '); writeln; write ('Place as many entries as desired on each line, but separate '); writeln ('entries with '); writeln ('at least one blank.'); write ('The preconditioner, C inverse, should follow in '); writeln(' the same way.'); writeln ('The initial approximation should also follow in same format.' ); writeln; writeln; writeln ('Has the input file been created? - enter Y or N. '); readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y' ) then begin writeln ('Input the file name in the form - drive:name.ext, '); writeln ('for example: A:DATA.DTA '); readln ( NAME ); assign ( INP, NAME ); reset ( INP ); OK := false; while ( not OK ) do begin writeln ('Input the number of equations - an integer. '); readln ( N ); if ( N > 0 ) then begin for I := 1 to N do for J := 1 to N + 1 do read ( INP, A[I,J] ); for I := 1 to N do for J := 1 to N do read ( INP, CI[I,J] ); for I := 1 to N do read ( INP, X1[I]); OK := true; close ( INP ) end else writeln ('The number must be a positive integer. ') end; OK := false; while ( not OK) do begin writeln ('Input the tolerance.'); readln ( TOL ); if (TOL > 0) then OK := true else writeln('Tolerance must be a positive number.') end; OK := false; while ( not OK) do begin writeln('Input maximum number of iterations.'); readln ( NN ); if (NN > 0) then OK := true else writeln('Number must be a positive integer.') end end else writeln ('The program will end so the input file can be created. ') end; procedure OUTPUT; begin writeln ('Choice of output method: '); writeln ('1. Output to screen '); writeln ('2. Output to text file '); writeln ('Please enter 1 or 2. '); readln ( FLAG ); if ( FLAG = 2 ) then begin writeln ('Input the file name in the form - drive:name.ext, '); writeln('for example: A:OUTPUT.DTA'); readln ( NAME ); assign ( OUP, NAME ) end else assign ( OUP, 'CON' ); rewrite ( OUP ); writeln(OUP,'CONJUGATE GRADIENT METHOD FOR LINEAR SYSTEMS'); writeln ( OUP ); writeln ( OUP, 'The solution vector is : '); for I := 1 to N do write ( OUP, X1[I]:12:8 ); writeln ( OUP ); writeln ( OUP, 'using ',K,' iterations with '); writeln ( OUP, 'Tolerance', TOL); writeln ( OUP, 'The Residual vector is : '); for I := 1 to N do write ( OUP, R[I]:12:8 ); close ( OUP ) end; begin INPUT; if ( OK ) then begin { STEP 1 } for I := 1 to N do for J := 1 to N do begin CT[I,J] := CI[J,I] end; for I := 1 to N do begin R[I] := A[I,N+1]; for J := 1 to N do R[I] := R[I] - A[I,J] * X1[J]; end; for I := 1 to N do begin W[I] := 0; for J := 1 to N do W[I] := W[I] + CI[I,J] * R[J]; end; for I := 1 to N do begin V[I] := 0; for J := 1 to N do V[I] := V[I]+CT[I,J]*W[J]; end; ALPHA := 0; for I := 1 to N do ALPHA := ALPHA + W[I]*W[I]; { STEP 2 } K := 1; OK := false; while ( not OK ) and ( K <= NN ) do begin ERR := 0.0; for I := 1 to N do ERR := ERR + V[I]*V[I]; { STEP 4 } if ( sqrt(ERR) < TOL) then begin OK := true; K := K-1 end else begin { STEP 5 } for I := 1 to N do begin U[I] := 0; for J := 1 to N do U[I] := U[I] + A[I,J]*V[J] end; S := 0; for I := 1 to N do S := S + V[I]*U[I]; T := ALPHA/S; for I := 1 to N do begin X1[I] := X1[I] + T*V[I]; R[I] := R[I] - T*U[I] end; for I := 1 to N do begin W[I] := 0; for J := 1 to N do W[I] := W[I] + CI[I,J]*R[J] end; BETA := 0; for I := 1 to N do BETA := BETA + W[I]*W[I]; { STEP 6 } if ( sqrt(BETA) <= TOL ) then begin ERR := 0; for I := 1 to N do ERR := ERR + R[I]*R[I]; ERR := sqrt(ERR); if (ERR < TOL) then OK := true end; if (not OK) then begin { STEP 7 } K := K + 1; S := BETA/ALPHA; for I := 1 to N do begin Z[I] := 0; for J := 1 to N do Z[I] := Z[I] + CT[I,J]*W[J] end; for I := 1 to N do V[I] := Z[I] + S*V[I]; ALPHA := BETA ; end end end; { STEP 8 } if ( not OK ) then writeln ('Maximum Number of Iterations Exceeded. ') else OUTPUT end end.