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.print(`This is the Conjugate Gradient Method for Linear Systems.`):OK := FALSE: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`): print(`A(1,1), A(1,2), ..., A(1,N+1), A(2,1), A(2,2), ..., A(2,N+1)`): print(`..., A(N,1), A(N,2), ..., A(N,N+1)\134n`): print(`Place as many entries as desired on each line, but separate `): print(`entries with`): print(`at least one blank.`): print(`Do the same for the input of the inverse of C.\134n`): print(`The initial approximation should follow in same format`): 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`): print(`for example: A:\134\134DATA.DTA`): NAME := scanf(`%s`)[1]: print(`The file name is`): print(NAME): INP := fopen(NAME,READ,TEXT): while OK = FALSE do print(`Input the number of equations - an integer.`): N := scanf(`%d`)[1]: print(`N is`): 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: for I1 from 1 to N do for J from 1 to N do CI[I1-1,J-1] := fscanf(INP, `%f`)[1]: CT[J-1,I1-1] := CI[I1-1,J-1]: od: od: for I1 from 1 to N do X1[I1-1] := fscanf(INP, `%f`)[1]: od: OK := TRUE: fclose(INP): else 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 of A in position `,I1,J): A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]): od: od: for I1 from 1 to N do for J from 1 to N do print(`input entry of the inverse matrix in position `,I1,J): CI[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(CI[I1-1,J-1]): CT[J-1,I1-1] := CI[I1-1,J-1]: od: od: for I1 from 1 to N do print(`input initial approximation entry in position `,I1): X1[I1-1] := scanf(`%f`)[1]:print(`Data is `):print(X1[I1-1]): od: OK := TRUE: else fi: od:fi:if OK = TRUE then OUP := default: fprintf(OUP, `The original system - output by rows:\134n`): for I1 from 1 to N do for J from 1 to N+1 do fprintf(OUP, ` %11.8f`, A[I1-1,J-1]): od: fprintf(OUP, `\134n`): od: fprintf(OUP, `The initial approximation:\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, X1[I1-1]): od: fprintf(OUP, `\134n`): fprintf(OUP, `The inverse - output by rows:\134n`): for I1 from 1 to N do for J from 1 to N do fprintf(OUP, ` %11.8f`, CI[I1-1,J-1]): od: fprintf(OUP, `\134n`): od:fi:if OK = TRUE then OK := FALSE: while OK = FALSE do print(`Input the tolerance.\134n`): TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL): if TOL > 0 then OK := TRUE: else print(`Tolerance must be a positive number.\134n`): fi: od:fi:if OK = TRUE then OK := FALSE: while OK = FALSE do print(`Input maximum number of iterations.\134n`): 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:fi:if OK = TRUE then# Step 1 for I1 from 1 to N do R[I1-1] := A[I1-1,N]: for J from 1 to N do R[I1-1] := R[I1-1]-A[I1-1,J-1]*X1[J-1]: od: od: for I1 from 1 to N do W[I1-1] := 0: for J from 1 to N do W[I1-1] := W[I1-1]+CI[I1-1,J-1]*R[J-1]: od: od: for I1 from 1 to N do V[I1-1] := 0: for J from 1 to N do V[I1-1] := V[I1-1]+CT[I1-1,J-1]*W[J-1]: od: od: ALPHA := 0.0: for I1 from 1 to N do ALPHA := ALPHA + W[I1-1]*W[I1-1]: od:# Step 2 K := 1: OK := FALSE:# Step 3 while (OK = FALSE) and (K <= NN) do ERR := 0: for I1 from 1 to N do ERR := ERR + V[I1-1]*V[I1-1]: od:# Step 4 if sqrt(ERR) < TOL then K := K -1: OK := TRUE: else# Step 5 for I1 from 1 to N do U[I1-1] := 0.0: for J from 1 to N do U[I1-1] := U[I1-1]+A[I1-1,J-1]*V[J-1]: od: od: SS := 0.0: for I1 from 1 to N do SS := SS + V[I1-1]*U[I1-1]: od: T := ALPHA/SS: for I1 from 1 to N do X1[I1-1] := X1[I1-1]+T*V[I1-1]: R[I1-1] := R[I1-1] - T*U[I1-1]: od: for I1 from 1 to N do W[I1-1] := 0.0: for J from 1 to N do W[I1-1] := W[I1-1]+CI[I1-1,J-1]*R[J-1]: od: od: BETA := 0.0: for I1 from 1 to N do BETA := BETA + W[I1-1]*W[I1-1]: od: ERR1 := sqrt(BETA):# Step 6 if ERR1 <= TOL then ERR := 0.0: for I1 from 1 to N do ERR := ERR + R[I1-1]*R[I1-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 I1 from 1 to N do Z[I1-1] := 0: for J from 1 to N do Z[I1-1] := Z[I1-1]+CT[I1-1,J-1]*W[J-1]: od: od: for I1 from 1 to N do V[I1-1] := Z[I1-1]+S*V[I1-1]: od: ALPHA := BETA: fi: fi: od:# Step 8 if OK = FALSE then print(`Maximum Number of Iterations Exceeded.`): else print(`Choice of output method:`): print(`1. Output to screen`): print(`2. Output to text file`): print(`Please enter 1 or 2.`): 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(`The output file is `):print(NAME): OUP := fopen(NAME, WRITE,TEXT): else OUP := default: fi: fprintf(OUP, `PRECONDITIONED CONJUGATE GRADIENT METHOD\134n\134n`): fprintf(OUP, `The solution vector is :\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, X1[I1-1]): od: fprintf(OUP, `\134nusing %d iterations with\134n`, K): fprintf(OUP, `Tolerance %.10e in infinity-norm\134n`, TOL): fprintf(OUP, `The residual vector is :\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, R[I1-1]): od: if OUP <> default then fclose(OUP): print(`Output file `,NAME,` created successfully`): fi: fi:fi:JSFH