(*GAUSSIAN ELIMINATION WITH BACKWARD SUBSTITUTION ALGORITHM 6.1
*
*   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 to solve a linear system\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"];
];
If[OK == 1,
   (* Step 1 *)
   NN = n-1;
   M = n+1;
   ICHG = 0;
   i = 1;
   While[OK == 1 && i <= NN,
      (* Step 2 - Use IP instead of P *)
      IP = i;
      While[Abs[A[IP-1,i-1]] <= 10^-20 && IP <= n,
	 IP = IP+1;
      ];
      If[IP == M,
	 OK = 0,
         (* Step 3 *)
	 If[IP != i,
	    For[JJ = 1,JJ <= M,JJ++,
               c = A[i-1,JJ-1];
	       A[i-1,JJ-1] = A[IP-1,JJ-1];
	       A[IP-1,JJ-1] = c;
            ];
	    ICHG = ICHG+1;
	 ];
         (* Step 4 *)
	 JJ = i+1;
	 For[J = JJ,J <= n,J++,
            (* Step 5 - Use XM in place of m(J,I) *)
	    XM = A[J-1,i-1]/A[i-1,i-1];
            (* Step 6 *)
	    For[K = JJ,K <= M,K++,
	       A[J-1,K-1] = A[J-1,K-1]-XM*A[i-1,K-1];
	    ];
            (* Multiplier XM could be saved in A[J,I] *)
	    A[J-1,i-1] = 0;
	 ];
      ];
      i = i+1;
   ];
   If[OK == 1,
      (* Step 7 *)
      If[Abs[A[n-1,n-1]] <= 10^-20,
	 OK = 0,
         (* Step 8 - Start backward substitution *)
	 X[n-1] = A[n-1,M-1]/A[n-1,n-1];
         (* Step 9 *)
	 For[K = 1,K <= NN,K++,
     	    i = NN-K+1;
	    JJ = i+1;
	    SUM = 0;
	    For[KK = JJ,KK <= n,KK++,
	       SUM = SUM-A[i-1,KK-1]*X[KK-1];
	    ];
	    X[i-1] = (A[i-1,M-1]+SUM)/A[i-1,i-1];
	 ];
         (* Step 10 - 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\n"]; 
	 Write[OUP,"\n"];
	 Write[OUP,"The reduced system - output by rows:\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"];
	 ];
	 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"];
	 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[];