(* ITERATIVE REFINEMENT ALGORITHM 7.4
*
*  To approximate the solution to the linear system Ax = b 
*    when A is suspected to be ill-conditioned: 
*
*   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), i<=i<=n of the inhomogeneous term b;
*	     tolerance TOL; maximum number of iterations NN.
*
*  Output:  the approximate solution XX(1), ..., XX(n) or a message
*	     that the maximum number of iterations was exceeded.
*)
CHIP[RND_,R1_,XXX_] := Module[{Y,OK,Z1,Z2,TEMP,Z,TEMP1,
     TN,SL,TT,ZZ,NN,MM,i,J,AT,RR},
   If[Abs[N[XXX]] < zero,
      Z = 0.0,
      i = 0;
      Y = N[XXX];
      RR = Floor[R1];
      AT = N[Log[10.0]];
      OK = 1;
      Z1 = Exp[N[R1]*AT];
      Z2 = Exp[(N[R1]-1)*AT];
      If[Abs[N[XXX]] >= Z1,
         While[OK == 1,
            Y = Y/10.0;
            i++;
            If[Abs[Y] < Z1,
               OK = 0;
            ];
         ],
         If[Abs[N[XXX]] < Z2,
            While[OK == 1,
               Y = Y*10.0;
               i--;
               If[Abs[N[Y]] >= Z2,
                  OK = 0;
               ];
            ];
         ];
      ];
      SL = Exp[-R1*AT];
      If[RND == 1,
         If[Y >= 0.0,
            TEMP = Y+0.5+0.1*SL,
            TEMP = Y-0.5-0.1*SL;
           ],
         If[Y >= 0.0,
            TEMP = Y+0.1*SL,
            TEMP = Y-0.1*SL;
           ];
        ];
      NN = Floor[RR/4.0];
      MM = Floor[4*NN];
      TEMP1[0] = TEMP;
      ZZ[0] = 0;
      TN[0] = 0;
      TT := TEMP;
      For[J=1, J <= NN, J++,
         TN[J] = Exp[MM*AT];
         TEMP1[J] = TT/TN[J];
         ZZ[J] = Floor[TEMP1[J]];
         If[ZZ[J] < 0,
            ZZ[J] = ZZ[J]+1
           ];
         TT = TT - ZZ[J]*TN[J];
         MM = MM - 4;
         ];
      ZZ[NN+1] = Floor[TT];
      If[ZZ[NN+1] < 0,
        ZZ[NN+1] = ZZ[NN+1] + 1
        ];
      TN[NN+1] = 1;
      Z = ZZ[1]*TN[1];
      For[J=2, J <= NN+1, J++,
         Z = Z + ZZ[J]*TN[J];
          ];
      Z = Z*Exp[i*AT];
    ];
  Z
];
zero = 0.0000000001;
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["Place as many entries as desired on each line, but \n"];
Print["separate entries with at least one blank\n"];
Print["\n"];
Print["\n"];
A1 = InputString["This is the Iterative Refinement Method\n
   The array 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[A1 == "yes" || A1 == "y" || A1 == "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];
   	    ];
	 ];
         OK = 1;
	 Close[INP],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   (* NN is used for the maximum number of iterations *)
   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
      ];
   ];
   OK = 0;
   RND = Input["Choice of rounding or chopping:\n
      1. Rounding\n
      2. Chopping\n
    Enter 1 or 2\n"];
   While[OK == 0,
      d = Input["Input the number of digits d <= 8 of rounding\n"];
      If[d > 0,
	 OK = 1,
	 Input["d must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      TOL = Input["Input tolerance, which is usually 10^(-d).\n"];
      If[TOL > 0,
	 OK = 1,
	 Input["Tolerance must be positive\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   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,"ITERATIVE REFINEMENT METHOD\n"];
   Write[OUP,"\n"];
   M = n+1; 
   Write[OUP,"Original system\n"];
   For[i = 1, i <= n, i++,
      For[J = 1, J <= M, J++,
	Write[OUP,SetPrecision[A[i-1,J-1],9]];
      ];
      Write[OUP,"\n"];
   ];
   If[RND == 1,
      Write[OUP,"Rounding to ",d," digits\n"],
      Write[OUP,"Chopping to ",d," digits\n"];
   ];
   Write[OUP,"\n"];
   Write[OUP,"Modified System\n"];
   For[i = 1, i <= n, i++,
      NROW[i-1] = i;
      For[J = 1, J <= M, J++,
         B[i-1,J-1] = CHIP[RND,d,A[i-1,J-1]];
         A[i-1,J-1] = B[i-1,J-1];
         Write[OUP,SetPrecision[A[i-1,J-1],9]];
      ];
      Write[OUP,"\n"];
   ];
   (* NROW and B have been initialized, Gauss elimination will begin *)
   (* step 0 *)
   i = 1;
   While[i <= n-1 && OK == 1,
      KK = i;
      While[Abs[A[KK-1,i-1]] < zero && KK <= n,
         KK = KK+1;
      ];
      If[KK > n,
         OK = 0;
         Write[OUP,"System does not have a unique solution.\n"],
         If[KK != i,
            (* rows interchange necessary *)
            IS = NROW[i-1];
            NROW[i-1] = NROW[KK-1];
            NROW[KK-1] = IS;
            For[J = 1,J <= M,J++,
               c = A[i-1,J-1];
               A[i-1,J-1] = A[KK-1,J-1];
               A[KK-1,J-1] = c;
            ];
         ];
         For[J = i+1,J <= n,J++,
            A[J-1,i-1] =
                CHIP[RND,d,A[J-1,i-1]/A[i-1,i-1]];
            For[L = i+1,L <= M,L++,
               Q1 = CHIP[RND,d,A[J-1,i-1]*
                         A[i-1,L-1]];
               A[J-1,L-1] = 
                    CHIP[RND,d,A[J-1,L-1]-Q1];
               ];
         ];
      ]; 
      i = i+1;
   ];
   If[Abs[A[n-1,n-1]] < zero && OK == 1,
      OK = 0;
      Write[OUP,"System has a singular matrix\n"];
   ];
   If[OK == 1,
      Write[OUP,"Reduced System\n"];
      For[i = 1,i <= n,i++,
         For[J = 1,J <= M,J++,
            Write[OUP,SetPrecision[A[i-1,J-1],9]];
         ];
         Write[OUP,"\n"];
      ];
      X[n-1] = CHIP[RND,d,A[n-1,M-1]/A[n-1,n-1]];
      For[i = 1,i <= n-1,i++,
         J = n-i;
         S = 0;
         For[L = J+1,L <= n,L++,
            S = CHIP[RND,d,S-
                CHIP[RND,d,A[J-1,L-1]*X[L-1]]];
         ];
         S = CHIP[RND,d,A[J-1,M-1]+S];
         X[J-1] = CHIP[RND,d,S/A[J-1,J-1]];
      ];
   ];
   Write[OUP,"Initial solution\n"];
   For[i = 1,i <= n,i++,
      Write[OUP,SetPrecision[X[i-1],9]];
   ];
   Write[OUP,"\n"];
   Write[OUP,"Row ordering is:\n"];
   For[i = 1, i <= n,i++,
      Write[OUP,NROW[i-1]];
   ];
   Write[OUP,"\n"];
   (* Refinement begins *)
   (* step 1 *)
   If[OK == 1,
      K = 1;
      For[i = 1,i <= n,i++,
         XX[i-1] = X[i-1];
      ];
   ];
   (* step 2 *)
   While[OK == 1 && K <= NN,
   (* LL is set to 1 if the desired accuracy in any component is not
      achieved.  Thus LL is initially 0 for each iteration *)
      LL = 0;
   (* step 3 *)
      For[i = 1,i <= n,i++,
         R[i-1] = 0.0;
         For[J = 1,J <= n,J++,
            Q1 = CHIP[RND,2*d,B[i-1,J-1]*XX[J-1]];
            R[i-1] = CHIP[RND,2*d,R[i-1]-Q1];
            ];
         R[i-1] = CHIP[RND,2*d,B[i-1,M-1]+R[i-1]];
      ];
      Write[OUP,"Residual Number ",K];
      For[i = 1,i <= n,i++,
         R[i-1] = CHIP[RND,d,R[i-1]];
         Write[OUP,SetPrecision[R[i-1],9]];
      ];
      Write[OUP,"\n"];
      (* step 4 *)
      (* Solve the linear system in the same order as in step 0. *)
      (*  The solution will be placed in X instead of Y. *)
      For[i = 1,i <= n-1,i++,
         I1 = NROW[i-1]; 
         For[J = i+1,J <= n,J++,
            J1 = NROW[J-1]; 
            Q1 = CHIP[RND,d,A[J-1,i-1]*R[I1-1]];
            R[J1-1] = CHIP[RND,d,R[J1-1]-Q1];
            ];
      ];
      I1 = NROW[n-1];
      X[n-1] = CHIP[RND,d,R[I1-1]/A[n-1,n-1]];
      For[i = 1,i <= n-1,i++,
         J = n-i;
         S = 0.0;
         For[L = J+1,L <= n,L++,
            Q1 = CHIP[RND,d,A[J-1,L-1]*X[L-1]];
            S = CHIP[RND,d,S-Q1];
         ];
         J1 = NROW[J-1]; 
         S = CHIP[RND,d,S+R[J1-1]];
         X[J-1] = CHIP[RND,d,S/A[J-1,J-1]];
      ];
      Write[OUP,"Vector Y\n"];
      For[i = 1,i <= n,i++,
         Write[OUP,SetPrecision[X[i-1],9]];
      ];
      Write[OUP,"\n"];
      (* steps 5 and 6 *)
      XXMAX = 0;
      YMAX = 0;
      ERR1 = 0;
      For[i = 1,i <= n,i++,
      (* if not accurate set LL to 1 *)
         If[Abs[X[i-1]] > TOL,
            LL = 1;
         ];
         If[K == 1,
            If[Abs[X[i-1]] > YMAX,
               YMAX = Abs[X[i-1]];
            ];
            If[Abs[XX[i-1]] > XXMAX,
               XXMAX = Abs[XX[i-1]];
            ];
         ];
         TEMP = XX[i-1];
         XX[i-1] = CHIP[RND,d,XX[i-1]+X[i-1]];
         TEMP = Abs[TEMP-XX[i-1]];
         If[TEMP > ERR1,
            ERR1 = TEMP;
         ];
      ];
      If[ERR1 <= TOL,
         LL = 2;
         ];
      If[K == 1,
         COND = YMAX/XXMAX*10.0^d;
      ];
      Write[OUP,"New approximation\n"];
      For[i = 1,i <= n,i++,
         Write[OUP,SetPrecision[XX[i-1],9]];
      ];
      Write[OUP,"\n"];
      (* step 7 *)
      If[LL == 0,
         Write[OUP,"The above vector is the solution.\n"];
         OK = 0,
         If[LL == 2,
            Write[OUP,"The above vector is the best possible with TOL = ",TOL];
            OK = 0,
            K = K+1;
         ];
      ];
      (* step 8 is not used in this implementation *)
   ];
   If[K > NN,
      Write[OUP,"Method fails\n"];
   ];
   Write[OUP,"Condition number is ",SetPrecision[COND,9]];
];
Quit[];