> restart:
> # 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 system has no unique solution.
> print(`This is Gaussian Elimination with Scaled Partial Pivoting.\n`):
> 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(`Has the input file been created? - enter 1 for yes or 2 for no.`):
> AA := scanf(`%d`)[1]: print(`Your response is`): print(AA):
> 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):
> OK := FALSE:
> 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:
> 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 in position `,I1,J): 
> A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]):
> od:
> od:
> OK := TRUE:
> else print(`The number must be a positive integer.\n`):
> 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:
> M := N+1:
> # Step 1
> for I1 from 1 to N do
> S[I1-1] := abs(A[I1-1,0]):
> # Initialize row pointer
> NROW[I1-1] := I1:
> for J from 1 to N do
> if abs(A[I1-1,J-1]) > S[I1-1] then
> S[I1-1] := abs(A[I1-1,J-1]):
> fi:
> od:
> if S[I1-1] <= 1.0e-20 then
> OK := FALSE:
> fi:
> od:
> NN := N-1:
> ICHG := 0:
> I1 := 1:
> # Step 2
> # Elimination process
> while OK = TRUE  and I1 <= NN do
> # Step 3
> IMAX := NROW[I1-1]:
> AMAX := abs(A[IMAX-1,I1-1])/S[IMAX-1]:
> IMAX := I1:
> JJ := I1+1:
> for IP from JJ to N do
> JP := NROW[IP-1]:
> TEMP := abs(A[JP-1,I1-1]/S[JP-1]):
> if TEMP > AMAX then
> AMAX := TEMP:
> IMAX := IP:
> fi:  
> od:
> # Step 4
> # System has no unique solution
> if AMAX <= 1.0e-20 then
> OK := FALSE:
> else
> # Step 5
> # Simulate row interchange
> if NROW[I1-1] <> NROW[IMAX-1] then
> ICHG := ICHG+1:
> NCOPY := NROW[I1-1]:
> NROW[I1-1] := NROW[IMAX-1]:
> NROW[IMAX-1] := NCOPY:
> fi:
> # Step 6
> I2 := NROW[I1-1]:
> for J from JJ to N do
> J1 := NROW[J-1]:
> # Step 7
> XM := A[J1-1,I1-1]/A[I2-1,I1-1]:
> # Step 8
> for K from JJ to M do
> A[J1-1,K-1] := A[J1-1,K-1]-XM*A[I2-1,K-1]:
> od:
> # Multiplier XM could be saved in A[J1-1,I1-1]
> A[J1-1,I1-1] := 0:
> od:
> fi:
> I1 := I1+1:
> od:
> if OK = TRUE then
> # Step 9
> N1 := NROW[N-1]:
> if abs(A[N1-1,N-1]) <= 1.0e-20 then
> OK := FALSE:
> # System has no unique solution
> else
> # Step 10
> # Start backward substitution
> X[N-1] := A[N1-1,M-1]/A[N1-1,N-1]:
> # Step 11
> for K from 1 to NN do
> I1 := NN-K+1:
> JJ := I1+1:
> N2 := NROW[I1-1]:
> SUM := 0:
> for KK from JJ to N do
> SUM := SUM-A[N2-1,KK-1]*X[KK-1]:
> od:
> X[I1-1] := (A[N2-1,N]+SUM)/A[N2-1,I1-1]:
> od:
> # Step 12
> # Process is complete
> print(`Choice of output method:\n`):
> print(`1. Output to screen\n`):
> print(`2. Output to text file\n`):
> print(`Please enter 1 or 2\n`):
> 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, `GAUSSIAN ELIMINATION WITH SCALED PARTIAL PIVOTING\n\n`):
> fprintf(OUP, `The reduced system - output by rows:\n`):
> for I1 from 1 to N do
> for J from 1 to M do
> fprintf(OUP, `  %11.8f`, A[I1-1,J-1]):
> od:
> fprintf(OUP,`\n`):
> od:
> fprintf(OUP, `\n\nHas solution vector:\n`):
> for I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I1-1]):
> od:
> fprintf(OUP, `\nwith %3d row interchange(s)\n`, ICHG) :
> fprintf(OUP, `\nThe rows have been logically re-ordered to:\n`):
> for I1 from 1 to N do
> fprintf(OUP, ` %2d`, NROW[I1-1]):
> od:
> fprintf(OUP, `\n`):
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi:
> fi:
> if OK = FALSE then
> print(`System has no unique solution\n`):
> fi:
> fi: