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.POWER := proc(X,M,OK,Y,B,YMU,TOL,NN,OUP) local K, LP, AMAX, I1, DONE, J, ERR, T:K := 1:LP := 1:AMAX := abs(X[0]):for I1 from 2 to M do if abs(X[I1-1]) > AMAX then AMAX := abs(X[I1-1]): LP := I1: fi:od:DONE := FALSE:for I1 from 1 to M do X[I1-1] := X[I1-1] / AMAX:od:while K <= NN and OK = TRUE and DONE = FALSE do for I1 from 1 to M do Y[I1-1] := 0: for J from 1 to M do Y[I1-1] := Y[I1-1] + B[I1-1,J-1] * X[J-1]: od: od: YMU := Y[LP-1]: LP := 1: AMAX := abs(Y[0]): for I1 from 2 to M do if abs(Y[I1-1]) > AMAX then AMAX := abs(Y[I1-1]): LP := I1: fi: od: if AMAX <= 1.0e-20 then print(`Zero eigenvalue - B is singular`): OK := FALSE: else ERR := 0: for I1 from 1 to M do T := Y[I1-1]/AMAX: if abs(X[I1-1]-T) > ERR then ERR := abs(X[I1-1]-T): fi: X[I1-1] := T: od: if ERR < TOL then for I1 from 1 to M do Y[I1-1] := X[I1-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 NN iterations. NN = `):print(NN): OK := FALSE:else fprintf(OUP, `Number Iterations for Power Method = %d\134n \134n`, K):fi:end:print(`This is Wielandt Deflation.`):print(`Choice of input method`):print(`1. input from keyboard - not recommended for large matrices`):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), A(2,1), A(2,2), ..., A(2,n),`): print(`..., A(n,1), A(n,2), ..., A(n,n)`): print(`Next place the approximate eigenvector V(1), ..., V(n) and follow it`): print(`by the approximate eigenvalue. Finally, an initial approximate `): print(`eigenvector of dimension n-1: X(1), ..., X(n-1) should follow.`): print(`Place as many entries as desired on each line, but separate `): print(`entries with 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 input is`): print(AA): if AA = 1 then print(`Input the file name in the form - drive:\134\134name.ext`): print(`for example: A:\134\134DATA.DTA`): NAME := scanf(`%s`)[1]: print(`Your input is`): print(NAME): INP := fopen(NAME,READ,TEXT): OK := FALSE: while OK = FALSE do print(`Input the dimension n.`): N := scanf(`%d`)[1]: if N > 1 then OK := TRUE else print(`Dimension must be greater than 1.`): fi: od: for I1 from 1 to N do for J from 1 to N do A[I1-1,J-1] := fscanf(INP, `%f`)[1]: od: od: OK := FALSE: for I1 from 1 to N do V[I1-1] := fscanf(INP, `%f`)[1]: if abs(V[I1-1]) > 0 then OK := TRUE: fi: od: XMU := fscanf(INP, `%f`)[1]: M := N-1: if OK = TRUE then OK := FALSE: for I1 from 1 to M do X[I1-1] := fscanf(INP, `%f`)[1]: if abs(X[I1-1]) > 0 then OK := TRUE: fi: od: fi: if OK = FALSE then print(`Input Error - All vectors must be nonzero.`): fi: fclose(INP): else print(`The program will end so the input file can be created.`): fi:else OK := FALSE: while OK = FALSE do print(`Input the dimension n - an integer.`): N := scanf(`%d`)[1]: print(`N= `): print(N): if N > 0 then print(`Input of the matrix A follows:`): for I1 from 1 to N do for J from 1 to N do print(`input entry in position `): print(I1,J): A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]): od: od: print(`Input of the approximate eigenvector V follows:`): OK := FALSE: while OK = FALSE do for I1 from 1 to N do print(`input entry in position `): print(I1): V[I1-1] := scanf(`%f`)[1]:print(`Data is `):print(V[I1-1]): if abs(V[I1-1]) > 0 then OK := TRUE: fi: od: if OK = TRUE then print(`Input the approximate eigenvalue`): XMU := scanf(`%f`)[1]:print(`Data is `):print(XMU): print(`Input of the initial approximate eigenvector X of dimension n-1 follows`): OK := FALSE: M:=N-1: for I1 from 1 to M do print(`input entry in position `): print(I1): X[I1-1] := scanf(`%f`)[1]:print(`Data is `):print(X[I1-1]): if abs(X[I1-1]) > 0 then OK := TRUE: fi: od: fi: if OK = FALSE then print(`Input Error - All vectors must be nonzero.`): fi: od: OK := TRUE: else print(`The number must be a positive integer.`): fi: od:fi:if OK = TRUE then OK := FALSE: while OK = FALSE do print(`Input a positive tolerance for the power method.`): TOL := scanf(`%f`)[1]: if TOL > 0 then OK := TRUE: else print(`Tolerance must be a positive number.`): fi: od: OK := FALSE: while OK = FALSE do print(`Input the maximum number of iterations for the `): print(`power method.`): NN := scanf(`%d`)[1]: if NN > 0 then OK := TRUE: else print(`The number must be a positive integer.`): fi: od: if OK = TRUE then OUP := default: fprintf(OUP, `The original matrix - output by rows:\134n`): for I1 from 1 to N do for J from 1 to N do fprintf(OUP, ` %11.8f`, A[I1-1,J-1]): od: fprintf(OUP, `\134n`): od: fprintf(OUP, `The approximate eigenvector:\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, V[I1-1]): od: fprintf(OUP, `\134n`): fprintf(OUP, `The approximate eigenvalue:\134n`): fprintf(OUP, ` %11.8f \134n`, XMU): fprintf(OUP, `The initial vector of dimension n-1:\134n`): for I1 from 1 to N-1 do fprintf(OUP, ` %11.8f`, X[I1-1]): od: fprintf(OUP, `\134n`): fi: if OK = TRUE then print(`Choice of output method:`): print(`1. Output to screen`): print(`2. Output to text file`): print(`Please enter 1 or 2.`): FLAG := scanf(`%d`)[1]: print(`Your input is`): print(FLAG): if FLAG = 2 then print(`Input the file name in the form - drive:\134\134name.ext`): print(`for example A:\134\134OUTPUT.DTA`): NAME := scanf(`%s`)[1]: print(`The output file is`): print(NAME): OUP := fopen(NAME,WRITE,TEXT): else OUP := default: fi: fprintf(OUP, `WIELANDT DEFLATION\134n\134n`):# Step 1 I1 := 1: AMAX := abs(V[0]): for J from 2 to N do if abs(V[J-1]) > AMAX then I1 := J: AMAX := abs(V[J-1]): fi: od:# Step 2 if I1 <> 1 then for K from 1 to I1-1 do for J from 1 to I1-1 do B[K-1,J-1] := A[K-1,J-1]-V[K-1]*A[I1-1,J-1]/V[I1-1]: od: od: fi:# Step 3 if I1 <> 1 and I1 <> N then for K from I1 to N-1 do for J from 1 to I1-1 do B[K-1,J-1] := A[K,J-1]-V[K]*A[I1-1,J-1]/V[I1-1]: B[J-1,K-1] := A[J-1,K]-V[J-1]*A[I1-1,K]/V[I1-1]: od: od: fi:# Step 4 if I1 <> N then for K from I1 to N-1 do for J from I1 to N-1 do B[K-1,J-1] := A[K,J]-V[K]*A[I1-1,J]/V[I1-1]: od: od: fi: POWER(X, M, OK, Y, B, YMU, TOL, NN, OUP): if OK = TRUE then# Step 6 if I1 <> 1 then for K from 1 to I1-1 do W[K-1] := Y[K-1]: od: fi:# Step 7 W[I1-1] := 0:# Step 8 if I1 <> N then for K from I1+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[I1-1,J-1] * W[J-1]: od: S := S/V[I1-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:\134n`): 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, `\134n`): od: fprintf(OUP, `\134nThe Eigenvalue = %12.8f`, YMU): fprintf(OUP, ` to Tolerance = %.10e\134n\134n`, TOL): fprintf(OUP, `Eigenvector is:\134n`): for I1 from 1 to N do fprintf(OUP,` %11.8f`, VV[I1-1]): od: fprintf(OUP, `\134n`): fi: if OUP <> default then fclose(OUP): print(`Output file `,NAME,` created successfully`): fi: fi:fi:JSFH