alg094.txt

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

> restart;
> # WIELANDT'S DEFLATION ALGORITHM 9.4
> #
> # To approximate the second most dominant eigenvalue and an
> # associated eigenvector of the n by n matrix A given an
> # approximation LAMBDA to the dominant eigenvalue, an
> # approximation V to a corresponding eigenvector and a vector X
> # belonging to R^(n-1), tolerance TOL, maximum number of
> # iterations N.
> #
> # INPUT:   Dimension n; matrix A; approximate eigenvalue LAMBDA;
> #          approximate eigenvector V belonging to R^n; vector X
> #          belonging to R^(n-1).
> #
> # OUTPUT:  Approximate eigenvalue MU; approximate eigenvector U or
> #          a message that the method fails.
> alg094 := proc() local POWER, OK, AA, NAME, INP, N, TOL, NN, I, J, A, V, XMU, YMU, M, X, FLAG, OUP, AMAX, K, B, W, S, VV, L1, L2, Y;
> POWER := proc(X,M,OK,Y,B,YMU,TOL,NN,OUP) local K, LP, AMAX, I, DONE, J, ERR, T;
> K := 1;
> LP := 1;
> AMAX := abs(X[0]);
> for I from 2 to M do
> if abs(X[I-1]) > AMAX then
> AMAX := abs(X[I-1]);
> LP := I;
> fi;
> od;
> DONE := FALSE;
> for I from 1 to M do
> X[I-1] := X[I-1] / AMAX;
> od;
> while K <= NN and OK = TRUE and DONE = FALSE do
> for I from 1 to M do
> Y[I-1] := 0;
> for J from 1 to M do
> Y[I-1] := Y[I-1] + B[I-1,J-1] * X[J-1];
> od;
> od;
> YMU := Y[LP-1];
> LP := 1;
> AMAX := abs(Y[0]);
> for I from 2 to M do
> if abs(Y[I-1]) > AMAX then
> AMAX := abs(Y[I-1]);
> LP := I;
> fi;
> od;
> if AMAX <= 1.0e-20 then
> printf(`Zero eigenvalue - B is singular\n`);
> OK := FALSE;
> else
> ERR := 0;
> for I from 1 to M do
> T := Y[I-1]/Y[LP-1];
> if abs(X[I-1]-T) > ERR then
> ERR := abs(X[I-1]-T);
> fi;
> X[I-1] := T;
> od;
> if ERR < TOL then
> for I from 1 to M do
> Y[I-1] := X[I-1];
> od;
> DONE := TRUE;
> else
> K := K+1;
> fi;
> fi;
> od;
> if K > NN and OK = TRUE then
> printf(`Power Method did not converge in %d iterations.\n`,NN);
> OK := FALSE;
> else
> fprintf(OUP, `Number Iterations for Power Method = %d\n \n`, K);
> fi;
> end;
> printf(`This is Wielandt Deflation.\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), A(2,1), A(2,2), ..., A(2,n),\n`);
> printf(`..., A(n,1), A(n,2), ..., A(n,n)\n\n`);
> printf(`Next place the approximate eigenvector V(1), ..., `);
> printf(`V(n) and follow it\n`);
> printf(`by the approximate eigenvalue. Finally, an `);
> printf(`initial approximate\n`);
> printf(`eigenvector of dimension n-1: X(1), ..., X(n-1) `);
> printf(`should follow.\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\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 dimension n.\n`);
> N := scanf(`%d`)[1];
> if N > 1 then
> OK := TRUE
> else
> printf(`Dimension must be greater than 1.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input a positive tolerance for the power method.\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be a positive number.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of iterations for the `);
> printf(`power method.\n`);
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`The number must be a positive integer.\n`);
> fi;
> od;
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> od;
> od;
> OK := FALSE;
> for I from 1 to N do
> V[I-1] := fscanf(INP, `%f`)[1];
> if abs(V[I-1]) > 0 then
> OK := TRUE;
> fi;
> od;
> XMU := fscanf(INP, `%f`)[1];
> M := N-1;
> if OK = TRUE then
> OK := FALSE;
> for I from 1 to M do
> X[I-1] := fscanf(INP, `%f`)[1];
> if abs(X[I-1]) > 0 then
> OK := TRUE;
> fi;
> od;
> fi;
> if OK = FALSE then
> printf(`Input Error - All vectors must be nonzero.\n`);
> fi;
> fclose(INP);
> else
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> 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, `WIELANDT DEFLATION\n\n`);
> # Step 1
> I := 1;
> AMAX := abs(V[0]);
> for J from 2 to N do
> if abs(V[J-1]) > AMAX then
> I := J;
> AMAX := abs(V[J-1]);
> fi;
> od;
> # Step 2
> if I <> 1 then
> for K from 1 to I-1 do
> for J from 1 to I-1 do
> B[K-1,J-1] := A[K-1,J-1]-V[K-1]*A[I-1,J-1]/V[I-1];
> od;
> od;
> fi;
> # Step 3
> if I <> 1 and I <> N then
> for K from I to N-1 do
> for J from 1 to I-1 do
> B[K-1,J-1] := A[K,J-1]-V[K]*A[I-1,J-1]/V[I-1];
> B[J-1,K-1] := A[J-1,K]-V[J-1]*A[I-1,K]/V[I-1];
> od;
> od;
> fi;
> # Step 4
> if I <> N then
> for K from I to N-1 do
> for J from I to N-1 do
> B[K-1,J-1] := A[K,J]-V[K]*A[I-1,J]/V[I-1];
> od;
> od;
> fi;
> POWER(X, M, OK, Y, B, YMU, TOL, NN, OUP);
> if OK = TRUE then
> # Step 6
> if I <> 1 then
> for K from 1 to I-1 do
> W[K-1] := Y[K-1];
> od;
> fi;
> # Step 7
> W[I-1] := 0;
> # Step 8
> if I <>  N then
> for K from I+1 to N do
> W[K-1] := Y[K - 2];
> od;
> fi;
> # Step 9
> S := 0;
> for J from 1 to N do
> S := S + A[I-1,J-1] * W[J-1];
> od;
> S := S/V[I-1];
> for K from 1 to N do
> # Compute eigenvector
> # VV is used in place of u.
> VV[K-1] := (YMU-XMU)*W[K-1]+S*V[K-1];
> od;
> fprintf(OUP, `The reduced matrix B:\n`);
> for L1 from 1 to M do
> for L2 from 1 to M do
> fprintf(OUP, `%.10e  `, B[L1-1,L2-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> fprintf(OUP, `\nThe Eigenvalue = %12.8f`, YMU);
> fprintf(OUP, ` to Tolerance = %.10e\n\n`, TOL);
> fprintf(OUP, `Eigenvector is:\n`);
> for I from 1 to N do
> fprintf(OUP,` %11.8f`, VV[I-1]);
> od;
> fprintf(OUP, `\n`);
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg094();