> 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)\n`):
>    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.\n`):
>    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:\\name.ext`):
>       print(`for example:   A:\\DATA.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 
              print(`The number must be a positive integer.\n`):
>          fi:
>       od:
>    else 
>       print(`The program will end so the input file can be created.\n`):
>    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 
           print(`The number must be a positive integer.`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OUP := default:
>    fprintf(OUP, `The original system - output by rows:\n`):
>    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, `\n`):
>    od:
>    fprintf(OUP, `The initial approximation:\n`):
>    for I1 from 1 to N do
>    fprintf(OUP, ` %11.8f`, X1[I1-1]):
>    od:
>    fprintf(OUP, `\n`):
>    fprintf(OUP, `The inverse - output by rows:\n`):
>    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, `\n`):
>    od:
> fi:
> if OK = TRUE then
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input the tolerance.\n`):
>       TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):
>       if TOL > 0 then
>          OK := TRUE:
>       else
>          print(`Tolerance must be a positive number.\n`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input maximum number of iterations.\n`):
>       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.\n`):
>       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:\\name.ext\n`):
>          print(`for example:   A:\\OUTPUT.DTA\n`):
>          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\n\n`):
>       fprintf(OUP, `The solution vector is :\n`):
>       for I1 from 1 to N do
>       fprintf(OUP, ` %11.8f`, X1[I1-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 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: