alg092.txt

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

> restart;
> # SYMMETRIC POWER METHOD ALGORITHM 9.2
> #
> # To approximate the dominant eigenvalue and an associated
> # eigenvector of the n by n symmetric matrix A given a nonzero vector x:
> #
> # INPUT:   Dimension n; matrix A; vector x; tolerance TOL;
> #          maximum number of iterations N.
> #
> # OUTPUT:  Approximate eigenvalue MU; approximate eigenvector x or
> #          a message that the maximum number of iterations was
> #          exceeded.
> alg092 := proc() local OK, AA, NAME, INP, N, I, J, A, Y, X, TOL, NN, FLAG, OUP, K, XL, ERR, T, YMU;
> printf(`This is the Symmetric Power Method.\n`);
> OK := FALSE;
> 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(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\n`);
> printf(`The initial approximation should follow in same format.\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 dimension n.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> 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;
> # The initial input is into Y and X is initialized as the zero vector.
> for I from 1 to N do
> Y[I-1] := fscanf(INP, `%f`)[1];
> od;
> for I from 1 to N do
> X[I-1] := 0;
> od;
> fclose(INP);
> while OK = FALSE do
> printf(`Input the tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be positive number.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations `);
> printf(`- integer.\n`);
> # NN is used in place of N.
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Number must be positive integer.\n`);
> fi;
> od;
> else
> printf(`The dimension 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 
> 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, `SYMMETRIC POWER METHOD\n\n`);
> fprintf(OUP, `iter  approx        approx eigenvector\n`);
> fprintf(OUP, `     eigenvalue\n`);   
> # Step 1
> K := 1;
> XL := 0;
> for I from 1 to N do
> XL := XL+Y[I-1]*Y[I-1];
> od;
> # 2-Norm of Y
> XL := sqrt(XL);
> ERR := 0;
> if XL > 0 then
> for I from 1 to N do
> T := Y[I-1]/XL;
> ERR := ERR+(X[I-1]-T)*(X[I-1]-T);
> X[I-1] := T;
> od;
> # X has a 2-Norm of 1.0
> ERR := sqrt(ERR);
> else
> printf(`A has a zero eigenvalue - select new vector and begin again\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> # Step 2
> while K <= NN and OK = TRUE do
> # Steps 3 and 4
> YMU := 0;
> for I from 1 to N do
> Y[I-1] := 0;
> for J from 1 to N do
> Y[I-1] := Y[I-1]+A[I-1,J-1]*X[J-1];
> od;
> YMU := YMU+X[I-1]*Y[I-1];
> od;
> # Steps 5 and 6
> XL := 0;
> for I from 1 to N do
> XL := XL+Y[I-1]*Y[I-1];
> od;
> # 2-Norm of Y
> XL := sqrt(XL);
> ERR := 0;
> if XL > 0 then
> for I from 1 to N do
> T := Y[I-1]/XL;
> ERR := ERR+(X[I-1]-T)*(X[I-1]-T);
> X[I-1] := T;
> od;
> # X has a 2-Norm of 1.0
> ERR := sqrt(ERR);
> else
> printf(`A has a zero eigenvalue - select new vector and begin again\n`);
> OK := FALSE;
> fi;
> fprintf(OUP, `%d %12.8f`, K, YMU);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> if OK = TRUE then
> # Step 7
> if  ERR < TOL then
> # Procedure completed successfully.
> fprintf(OUP, `\n\nThe eigenvalue = %12.8f`,YMU);
> fprintf(OUP, ` to tolerance = %.10e\n`, TOL);
> fprintf(OUP, `obtained on iteration number = %d\n\n`, K);
> fprintf(OUP, `Unit eigenvector is :\n\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> OK := FALSE;
> else
> # Step 8
> K := K+1;
> fi;
> fi;
> od;
> # Step 9
> if K > NN then
> printf(`No convergence within %d iterations\n`, NN);
> fi;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg092();