(*GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING ALGORITHM 6.3
*
*   To solve the n by n linear system
*
* E1: A[1,1] X[1] + A[1,2] X[2] + ... + A[1,n] X[n] = A[1,n+1] 
* E2: A[2,1] X[1] + A[2,2] X[2] + ... + A[2,n] X[n] = A[2,n+1] 
*  :
*  .
* EN: A[n,1] X[1] + A[n,2] X[2] + ... + A[n,n] X[n] = A[n,n+1] 
*
*   Input:  number of unknowns and equations n; augmented
*  	     matrix A = (A(i,J)) where 1<=i<=n and 1<=J<=n+1.
*
*  Output:  solution x(1), x(2), ..., x(n) or a message that the
* 	     linear equation has no unique solution.
*)
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"];
AA = InputString["This is Gaussian Elimination with Scaled Partial Pivoting\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"];
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 - 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"];
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
   OK = 0;
];
If[OK == 1,
   M = n+1;
   (* Step 1 *)
   For[i = 1,i <= n,i++,
      S[i-1] = Abs[A[i-1,0]];
      (* Initialize row pointer *)
      NROW[i-1] = i;
      For[J = 1,J <= n,J++,
	 If[Abs[A[i-1,J-1]] > S[i-1],
	    S[i-1] = Abs[A[i-1,J-1]];
	 ];
      ];
      If[S[i-1] <= 10^-20,
         OK = 0;
      ];
   ];
   NN = n-1;
   ICHG = 0;
   i = 1;
   (* Step 2 - Elimination process *)
   While[OK == 1 && i <= NN,
      (* Step 3 *)
      IMAX = NROW[i-1];
      AMAX = Abs[A[IMAX-1,i-1]]/S[IMAX-1];
      IMAX = i;
      JJ = i+1;
      For[IP = JJ,IP <= n,IP++,
	 JP = NROW[IP-1];
	 TEMP = Abs[A[JP-1,i-1]/S[JP-1]];
	 If[TEMP > AMAX,
	    AMAX = TEMP;
	    IMAX = IP;
	 ];
      ];
      (* Step 4 - System has no unique solution *)
      If[AMAX <= 10^-20,
         OK = 0,
         (* Step 5 - Simulate row interchange *)
	 If[NROW[i-1] != NROW[IMAX-1],
	    ICHG = ICHG+1;
	    NCOPY = NROW[i-1];
	    NROW[i-1] = NROW[IMAX-1];
	    NROW[IMAX-1] = NCOPY;
	 ];
         (* Step 6 *)
	 I1 = NROW[i-1];
	 For[J = JJ,J <= n,J++,
	    J1 = NROW[J-1];
            (* Step 7 *)
	    XM = A[J1-1,i-1]/A[I1-1,i-1];
            (* Step 8 *)
	    For[K = JJ,K<=M,K++,
	       A[J1-1,K-1] = A[J1-1,K-1]-XM*A[I1-1,K-1];
	    ];
            (* Multiplier XM could be saved in A[J1-I,I-1] *)
	    A[J1-1,i-1] = 0;
	 ];
      ];
      i = i+1;
   ];
   If[OK == 1,
      (* Step 9 *)
      N1 = NROW[n-1];
      If[Abs[A[N1-1,n-1]] <= 10^-20,
	 OK = 0,
         (* System has no unique solution *)
         (* Step 10 - Start backward substitution *)
	 X[n-1] = A[N1-1,M-1]/A[N1-1,n-1];
         (* Step 11 *)
	 For[K = 1,K <= NN,K++,
	    i = NN-K+1;
	    JJ = i+1;
	    N2 = NROW[i-1];
	    SUM = 0;
	    For[KK = JJ,KK <= n,KK++,
	       SUM = SUM-A[N2-1,KK-1]*X[KK-1];
	    ];
     	    X[i-1] = (A[N2-1,n]+SUM)/A[N2-1,i-1];
   	 ];
         (* Step 12 - Procedure completed successfully *)
         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,"GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING\n"]; 
         Write[OUP,"\n"];
         Write[OUP,"The reduced system - output by rows:\n"];
         For[i = 1,i <= n,i++,
	    For[J = 1,J <= n,J++,
               Write[OUP,SetPrecision[A[i-1,J-1],9]];
            ];
            Write[OUP,"\n"];
         ];
	 Write[OUP,"\n"];
	 Write[OUP,"Has solution vector:\n"];
	 For[i = 1,i <= n,i++,
	    Write[OUP,SetPrecision[X[i-1],9]];
         ];
	 Write[OUP,"With ",ICHG," row interchange(s)\n"];
	 Write[OUP,"\n"];
	 Write[OUP,"The rows have been logically re-ordered to:\n"];	
	 For[i = 1,i <= n,i++,
	    Write[OUP,NROW[i-1]];
	 ];
	 Write[OUP,"\n"];
         If[OUP == "OutputStream[",NAME," 3]",
            Print["Output file: ",NAME," created successfully\n"];
            Close[OUP];
         ];
      ];
   ];
   If[OK == 0,
      Print["System has no unique solution\n"];
   ];
];
Quit[];