(* 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 CI(i,j), 1<=i, j<=n of the matrix C inverse;
*	     The entries B(i), i<=i<=n of the inhomogeneous term b;
*	     the entries XO(i), 1<=i<=n, of x(0);
*	     tolerance TOL; maximum number of iterations NN;
*
*  Output:  the approximate solution X(1), ..., X(n) or a message
*	     that the maximum number of iterations was exceeded.
*)
Print["\n"];
Print["A(1,1), A(2,1), ..., A(1,n+1), A(2,1), A(2,2), ...\n"];
Print["A(2,n+1), ..., A(n,1), A(n,2), ..., A(n,n+1)\n"];
Print["\n"];
Print["The preconditioner, C inverse, should follow in\n"];
Print["the same way.  The initial approximation should also\n"];
Print["follow in the same format.\n"];
Print["\n"];
Print["Place as many entries as desired on each line, but \n"];
Print["separate entries with at least one blank\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is the Conjugate Gradient Method for Linear Systems\n
     The data will be input from a text file in the order:(see screen)\n
     Has the input file been created?\n
     Enter 'yes' or 'no'\n"];
OK = 0;
If[AA == "yes" || AA == "y" || AA == "Y",
   NAME = InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA.DTA\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the number of equations n - an integer\n"];
      If[n > 0,
	 For[i = 1,i <= n,i++,
	    For[J = 1,J <= n+1,J++,
	       A[i-1,J-1] = Read[INP,Number];
   	    ];
	 ];
	 For[i = 1,i <= n,i++,
            For[J = 1,J <= n, J++,
               CI[i-1,J-1] = Read[INP,Number];
   	    ];
	 ];
	 For[i = 1, i <= n, i++,
	    X1[i-1] = Read[INP,Number];
	 ];
         (* Use X1 for X0 *)
         OK = 1;
	 Close[INP],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      TOL = Input["Input the tolerance\n"];
      If[TOL <= 0,
         Input["Tolerance must be positive.\n
         Enter 1 [enter] to continue\n"],
         OK = 1;
      ];
   ];
   OK = 0;
   While[OK == 0,
      NN = Input["Input maximum number of iterations\n
        no decimal points\n"];
      If[NN <= 0,
         Input["Must be a positive integer.\n
         Enter 1 [enter] to continue\n"],
         OK = 1
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   (* Step 1 *)
   OUP = "stdout";
   For[i = 1,i <= n,i++,
      For[J = 1,J <= n, J++,
         CT[i-1,J-1] = CI[J-1,i-1];
      ];
   ];
   For[i = 1,i <= n,i++,
      R[i-1] = A[i-1,n];
      For[J = 1,J <= n, J++,
         R[i-1] = R[i-1] - A[i-1,J-1]*X1[J-1];
      ];
   ];
   For[i = 1,i <= n,i++,
      W[i-1] = 0.0;
      For[J = 1,J <= n, J++,
         W[i-1] = W[i-1] + CI[i-1,J-1]*R[J-1];
      ];
   ];
   For[i = 1,i <= n,i++,
      V[i-1] = 0.0;
      For[J = 1,J <= n, J++,
         V[i-1] = V[i-1] + CT[i-1,J-1]*W[J-1];
      ];
   ];
   ALPHA = 0.0;
   For [i = 1, i <= n, i++,
       ALPHA = ALPHA + W[i-1]*W[i-1];
   ];
   (* Step 2 *)
   K = 1;
   OK = 0;
   (* Step 3 *)
   While[OK == 0 && K <= NN,
      (* ERR is used to test accuracy - it measures the 2-norm *)
      ERR = 0;
      (* Step 4 *)
      For[i = 1, i <= n, i++,
         ERR = ERR + V[i-1]*V[i-1];
      ];
      If[Sqrt[ERR] < TOL,
         OK = 1;
         K = K -1,
      (* Step 5 *)
      For[i = 1,i <= n,i++,
         U[i-1] = 0.0;
         For[J = 1,J <= n, J++,
            U[i-1] = U[i-1] + A[i-1,J-1]*V[J-1];
         ];
      ];
      S = 0.0;
      For[i = 1, i <= n, i++,
         S = S + V[i-1]*U[i-1];
      ];
      T = ALPHA/S;
      For[i = 1, i <= n, i++,
         X1[i-1] = X1[i-1] + T*V[i-1];
         R[i-1]  = R[i-1]  - T*U[i-1];
      ];
      Write[OUP,"The approximate solution vector is:\n"];
      For[i = 1, i <= n, i++,
         Write[OUP,SetPrecision[X1[i-1],9]];
      ];
      Write[OUP,"\n The residual vector is:\n"];
      For[i = 1, i <= n, i++,
         Write[OUP,SetPrecision[R[i-1],9]];
      ];
      Write[OUP,"\n"];
      For[i = 1,i <= n,i++,
         W[i-1] = 0.0;
         For[J = 1,J <= n, J++,
            W[i-1] = W[i-1] + CI[i-1,J-1]*R[J-1];
         ];
      ];
   BETA = 0.0;
   For [i = 1, i <= n, i++,
       BETA = BETA + W[i-1]*W[i-1];
   ];
   (* Step 6 *)
   If[Sqrt[BETA] < TOL,
     ERR = 0.0;
     For[i=1, i <=n, i++,
       ERR = ERR + R[i-1]*R[i-1];
     ];
     If[Sqrt[ERR] < TOL,
       OK = 1;
     ];
   ];
   (* Step 7 *)
   If[OK == 0,
      K = K + 1;
      S = BETA/ALPHA;
      For[i = 1,i <= n,i++,
         Z[i-1] = 0.0;
         For[J = 1,J <= n, J++,
            Z[i-1] = Z[i-1] + CT[i-1,J-1]*W[J-1];
         ];
      ];
     For[i=1, i <=n, i++,
       V[i-1] = Z[i-1] + S*V[i-1];
     ];
     ALPHA = BETA;
   ];
 ];
 ];
   If[OK == 0,
      Print["Maximum number of iterations exceeded\n"],
      (* Step 8 - Procedure completed unsuccessfully *)
      FLAG = Input["Select output destination\n
        1. Screen\n
        2. Text file\n
        Enter 1 or 2\n"];
      If[FLAG == 2,
         NAME = InputString["Input the file name\n
            For example:   output.dta\n"];
         OUP = OpenWrite[NAME,FormatType->OutputForm],
         OUP = "stdout";
      ];
      Write[OUP,"CONJUGATE GRADIENT METHOD FOR LINEAR SYSTEMS\n"];
      Write[OUP,"\n"];
      Write[OUP,"The solution vector is:\n"];
      For[i = 1, i <= n, i++,
         Write[OUP,SetPrecision[X1[i-1],9]];
      ];
      Write[OUP,"\n The residual vector is:\n"];
      For[i = 1, i <= n, i++,
         Write[OUP,SetPrecision[R[i-1],9]];
      ];
      Write[OUP,"\n"];
      Write[OUP,"Using ",K," iterations\n"];
      Write[OUP,"with tolerance ",SetPrecision[TOL,9]," in the 2-norm\n"];
      If[OUP == "OutputStream[",NAME," 3]",
         Print["Output file: ",NAME," created successfully\n"];
         Close[OUP];
      ];
   ];
];	
Quit[];