alg063.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg063 := proc() local AA, NAME, INP, OK, N, I, J, A, M, S, NROW, NN, ICHG, IMAX, AMAX, JJ, IP, JP, TEMP, NCOPY, I1, J1, XM, K, N1, X, N2, SUM, KK, FLAG, OUP;
> printf(`This is Gaussian Elimination with Scaled Partial Pivoting.\n`);
> printf(`The array will be input from a text file in the order:\n`);
> printf(`A(1,1), A(1,2), ..., A(1,N+1), A(2,1), A(2,2), ..., 
> A(2,N+1),\n`);
> printf(`..., A(N,1), A(N,2), ..., A(N,N+1)\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> AA := scanf(`%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number of equations - an integer.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> for I from 1 to N do
> for J from 1 to N+1 do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> od;
> od;
> OK := TRUE;
> fclose(INP);
> else printf(`The number must be a positive integer.\n`);
> fi;
> od;
> else 
> printf(`The program will end so the input file can be created.\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> M := N+1;
> # Step 1
> for I from 1 to N do
> S[I-1] := abs(A[I-1,0]);
> # Initialize row pointer
> NROW[I-1] := I;
> for J from 1 to N do
> if abs(A[I-1,J-1]) > S[I-1] then
> S[I-1] := abs(A[I-1,J-1]);
> fi;
> od;
> if S[I-1] <= 1.0e-20 then
> OK := FALSE;
> fi;
> od;
> NN := N-1;
> ICHG := 0;
> I := 1;
> # Step 2
> # Elimination process
> while OK = TRUE  and I <= NN do
> # Step 3
> IMAX := NROW[I-1];
> AMAX := abs(A[IMAX-1,I-1])/S[IMAX-1];
> IMAX := I;
> JJ := I+1;
> for IP from JJ to N do
> JP := NROW[IP-1];
> TEMP := abs(A[JP-1,I-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[I-1] <> NROW[IMAX-1] then
> ICHG := ICHG+1;
> NCOPY := NROW[I-1];
> NROW[I-1] := NROW[IMAX-1];
> NROW[IMAX-1] := NCOPY;
> fi;
> # Step 6
> I1 := NROW[I-1];
> for J from JJ to N do
> J1 := NROW[J-1];
> # Step 7
> XM := A[J1-1,I-1]/A[I1-1,I-1];
> # Step 8
> for K from JJ to M do
> A[J1-1,K-1] := A[J1-1,K-1]-XM*A[I1-1,K-1];
> od;
> # Multiplier XM could be saved in A[J1-1,I-1]
> A[J1-1,I-1] := 0;
> od;
> fi;
> I := I+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
> I := NN-K+1;
> JJ := I+1;
> N2 := NROW[I-1];
> SUM := 0;
> for KK from JJ to N do
> SUM := SUM-A[N2-1,KK-1]*X[KK-1];
> od;
> X[I-1] := (A[N2-1,N]+SUM)/A[N2-1,I-1];
> od;
> # Step 12
> # Process is complete
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`For example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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 I from 1 to N do
> for J from 1 to M do
> fprintf(OUP, `  %11.8f`, A[I-1,J-1]);
> od;
> fprintf(OUP,`\n`);
> od;
> fprintf(OUP, `\n\nHas solution vector:\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\nwith %3d row interchange(s)\n`, ICHG) ;
> fprintf(OUP, `\nThe rows have been logically re-ordered to:\n`);
> for I from 1 to N do
> fprintf(OUP, ` %2d`, NROW[I-1]);
> od;
> fprintf(OUP, `\n`);
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> fi;
> if OK = FALSE then
> printf(`System has no unique solution\n`);
> fi;
> fi;
> RETURN(0);
> end;
> alg063();