alg062.txt

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

> restart;
> # GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING ALGORITHM 6.2
> #
> # 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.
> alg062 := proc() local AA, NAME, INP, OK, N, I, J, A, M, NROW, NN, ICHG, IMAX, AMAX, JJ, IP, JP, NCOPY, I1, J1, XM, K, N1, X, N2, SUM, KK, FLAG, OUP;
> printf(`This is Gaussian Elimination with 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`);
> fi;
> if OK = TRUE then
> M := N+1;
> # Step 1
> for I from 1 to N do
> NROW[I-1] := I;
> od;
> # Initialize row pointer
> NN := N-1;
> ICHG := 0;
> I := 1;
> # Step 2
> while OK = TRUE and I <= NN do
> # Step 3
> IMAX := NROW[I-1];
> AMAX := abs(A[IMAX-1,I-1]);
> IMAX := I;
> JJ := I+1;
> for IP from JJ to N do
> JP := NROW[IP-1];
> if abs(A[JP-1,I-1]) > AMAX then
> AMAX := abs(A[JP-1,I-1]);
> IMAX := IP;
> fi;
> od;
> # Step 4
> 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;
> I1 := NROW[I-1];
> # Step 6
> 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 - 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, `  %12.8f`, X[I-1]);
> od;
> fprintf (OUP, `\n\nwith %d 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;
> alg062();