restart:

# INVERSE POWER METHOD ALGORITHM 9.3

# To approximate an eigenvalue and an associated eigenvector of the

# n by n 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.

MULTIP := proc(N,OK,NROW,Q,A) local I2, M, IMAX, J, IP, L1, L2, JJ, I1, J1, K:

# Procedure MULTIP determines the row ordering and multipliers for the

# matrix (A-Q*I)

for I1 from 1 to N do

NROW[I1-1] := I1:

od: OK := TRUE:

I1 := 1:

M := N - 1:

while I1 <= M and OK = TRUE do

IMAX := I1:

J := I1+1:

for IP from J to N do

L1 := NROW[IMAX-1]:

L2 := NROW[IP-1]:

if abs(A[L2-1,I1-1]) > abs(A[L1-1,I1-1]) then

IMAX := IP:

fi:

od:

if abs(A[NROW[IMAX-1]-1,I1-1]) <= 1.0e-20 then

OK := FALSE:

print(`A - Q * I is singular, Q is an eigenvalue, Q = `):print(Q):

else

JJ := NROW[I1-1]:

NROW[I1-1] := NROW[IMAX-1]:

NROW[IMAX-1] := JJ:

I2 := NROW[I1-1]:

for JJ from J to N do

J1 := NROW[JJ-1]:

A[J1-1,I1-1] := A[J1-1,I1-1] / A[I2-1,I1-1]:

for K from J to N do

A[J1-1,K-1] := A[J1-1,K-1] - A[J1-1,I1-1] * A[I2-1,K-1]:

od:

fi:

I1 := I1+1:

od:

if abs(A[NROW[N-1]-1,N-1]) <= 1.0e-20 then

OK := FALSE:

print(`A - Q * I is singular, Q is an eigenvalue, Q = `):print(Q):

fi:

end:

SOLVE := proc(N,B,A,Y,NROW) local M, I2, J, I1, JJ, J1, N1, L, K, N2, KK:

# Procedure SOLVE solves the linear system (A-Q*I)*Y=X given a new

# vector X and returns the solution in Y

M := N - 1:

for I2 from 1 to M do

J := I2+1:

I1 := NROW[I2-1]:

for JJ from J to N do:

J1 := NROW[JJ-1]:

B[J1-1] := B[J1-1] - A[J1-1,I2-1] * B[I1-1]:

od:

N1 := NROW[N-1]:

Y[N-1] := B[N1-1] / A[N1-1,N-1]:

L := N - 1:

for K from 1 to L do J := L - K + 1:

JJ := J + 1:

N2 := NROW[J-1]:

Y[J-1] := B[N2-1]:

for KK from JJ to N do

Y[J-1] := Y[J-1] - A[N2-1,KK-1] * Y[KK-1]:

od:

Y[J-1] := Y[J-1] / A[N2-1,J-1]:

od:

end:

print(`This is the Inverse Power Method.\134n`):

OK := FALSE:

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)\134n`):

print(`Place as many entries as desired on each line, but separate `):

print(`entries with`):

print(`at least one blank.`):

print(`The initial approximation should follow in same format.\134n`):

print(`Has the input file been created? - enter 1 for yes or 2 for no.`):

AA := scanf(`%d`)[1]: print(`Your response 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(`The file name is`): print(NAME):

INP := fopen(NAME,READ,TEXT):

OK := FALSE:

while OK = FALSE do

print(`Input the dimension n - an integer.`):

N := scanf(`%d`)[1]: print(`N is`): print(N):

if N > 0 then

for I1 from 1 to N do

for J from 1 to N do

A[I1-1,J-1] := fscanf(INP, `%f`)[1]:

od:

for I1 from 1 to N do

X[I1-1] := fscanf(INP, `%f`)[1]:

od:

OK := TRUE:

fclose(INP):

else print(`The number must be a positive integer.\134n`):

fi:

od:

else

print(`The program will end so the input file can be created.\134n`):

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

for I1 from 1 to N do

for J from 1 to N do

print(`input entry in position `,I1,J):

A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]):

od:

print(`Input initial approximation vector `):

for I1 from 1 to N do

print(`input entry in position `,I1):

X[I1-1] := scanf(`%f`)[1]:print(`Data is `):print(X[I1-1]):

od:

OK := TRUE:

else print(`The number must be a positive integer.\134n`):

fi:

od:

fi:

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 initial approximation vector is :\134n`):

for I1 from 1 to N do

fprintf(OUP, ` %11.8f`, X[I1-1]):

od:

fprintf(OUP, `\134n`):

OK := FALSE:

while OK = FALSE do

print(`Input the tolerance.\134n`):

TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):

if TOL > 0 then

OK := TRUE:

else

print(`Tolerance must be positive number.\134n`):

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input maximum number of iterations `):

print(`- integer.\134n`):

NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):

# Use NN in place of N

if NN > 0 then

OK := TRUE:

else

print(`Number must be positive integer.\134n`):

fi:

od:

fi:

if OK = TRUE then

# Step 1

# Q could be input instead of computed.

Q := 0:

S := 0:

for I1 from 1 to N do

S := S + X[I1-1] * X[I1-1]:

for J from 1 to N do

Q := Q + A[I1-1,J-1] * X[I1-1] * X[J-1]:

od:

Q := Q / S:

print(`Q = `, Q):

print(`Input new Q? Enter 1 for yes or 2 for no.`):

AA := scanf(`%d`)[1]:print(`Input is `):print(AA):

if AA = 1 then

print(`input new Q`):

Q := scanf(`%f`)[1]:print(`Q = `):print(Q):

fi:

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(`Input is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext\134n`):

print(`for example A:\134\134OUTPUT.DTA\134n`):

NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):

OUP := fopen(NAME,WRITE,TEXT):

else

OUP := default:

fi:

fprintf(OUP, `INVERSE POWER METHOD\134n\134n`):

fprintf(OUP, `Iteration Eigenvalue Eigenvector\134n`):

# Step 2

K := 1:

for I1 from 1 to N do

A[I1-1,I1-1] := A[I1-1,I1-1] - Q:

od:

# Call subroutine to compute multipliers M(I,J) and upper triangular

# matrix for the matrix (A-Q*I)

MULTIP(N, OK, NROW, Q, A):

if OK = TRUE then

# Step 3

LP := 1:

for I1 from 2 to N do

if abs(X[I1-1]) > abs(X[LP-1]) then

LP := I1:

fi:

od:

# Step 4

AMAX := X[LP-1]:

for I1 from 1 to N do

X[I1-1] := X[I1-1] / (AMAX):

od:

# Step 5

while K <= NN and OK = TRUE do

# Steps 6 and 7

for I1 from 1 to N do

B[I1-1] := X[I1-1]:

od:

# Subroutine SOLVE returns the solution of (A-Q*I)*Y=b in Y

SOLVE(N, B, A, Y, NROW):

# Step 8

YMU := Y[LP-1]:

# Steps 9 and 10

LP := 1:

for I1 from 2 to N do

if abs(Y[I1-1]) > abs(Y[LP-1]) then

LP := I1:

fi:

od:

AMAX := Y[LP-1]:

ERR := 0:

for I1 from 1 to N 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:

YMU := 1 / YMU + Q:

# Step 11

fprintf(OUP, `%3d %12.8f\134n`, K, YMU):

for I1 from 1 to N do

fprintf(OUP, ` %11.8f`, X[I1-1]):

od:

fprintf(OUP, `\134n`):

if ERR < TOL then

OK := FALSE:

fprintf(OUP, `Eigenvalue = %12.8f`, YMU):

fprintf(OUP, ` to tolerance = %.10e\134n`, TOL):

fprintf(OUP, `obtained on iteration number = %d\134n\134n`, K):

fprintf(OUP, `Unit eigenvector is :\134n`):

for I1 from 1 to N do

fprintf(OUP, ` %11.8f`, X[I1-1]):

od:

fprintf(OUP, `\134n`):

else

# Step 12

K := K+1:

fi:

od:

if K > NN then

fprintf(OUP, `No convergence in %d iterations\134n`,NN):

fi:

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created successfully`):

fi:

SUNUaGlzfmlzfnRoZX5JbnZlcnNlflBvd2Vyfk1ldGhvZC58K0c2Ig==

STdDaG9pY2V+b2Z+aW5wdXR+bWV0aG9kRzYi

SWZuMS5+aW5wdXR+ZnJvbX5rZXlib2FyZH4tfm5vdH5yZWNvbW1lbmRlZH5mb3J+bGFyZ2V+bWF0cmljZXNHNiI=

SToyLn5pbnB1dH5mcm9tfmF+dGV4dH5maWxlRzYi

STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig==

SS5Zb3VyfmlucHV0fmlzRzYi

IiIj

SVZUaGV+YXJyYXl+d2lsbH5iZX5pbnB1dH5mcm9tfmF+dGV4dH5maWxlfmlufnRoZX5vcmRlckc2Ig==

SWduQSgxLDEpLH5BKDEsMiksfi4uLix+QSgxLE4pLH5BKDIsMSksfkEoMiwyKSx+Li4uLH58K35+fkEoMixOKUc2Ig==

SUIuLi4sfkEoTiwxKSx+QShOLDIpLH4uLi4sfkEoTixOKXwrRzYi

SWduUGxhY2V+YXN+bWFueX5lbnRyaWVzfmFzfmRlc2lyZWR+b25+ZWFjaH5saW5lLH5idXR+c2VwYXJhdGV+RzYi

SS1lbnRyaWVzfndpdGhHNiI=

STRhdH5sZWFzdH5vbmV+YmxhbmsuRzYi

SVlUaGV+aW5pdGlhbH5hcHByb3hpbWF0aW9ufnNob3VsZH5mb2xsb3d+aW5+c2FtZX5mb3JtYXQufCtHNiI=

SWpuSGFzfnRoZX5pbnB1dH5maWxlfmJlZW5+Y3JlYXRlZD9+LX5lbnRlcn4xfmZvcn55ZXN+b3J+Mn5mb3J+bm8uRzYi

STFZb3VyfnJlc3BvbnNlfmlzRzYi

IiIi

SVJJbnB1dH50aGV+ZmlsZX5uYW1lfmlufnRoZX5mb3Jtfi1+ZHJpdmU6XG5hbWUuZXh0RzYi

STtmb3J+ZXhhbXBsZTp+fn5BOlxEQVRBLkRUQUc2Ig==

STFUaGV+ZmlsZX5uYW1lfmlzRzYi

US5FOlxBTEcwOTMuRFRBNiI=

SURJbnB1dH50aGV+ZGltZW5zaW9ufm5+LX5hbn5pbnRlZ2VyLkc2Ig==

SSVOfmlzRzYi

IiIk

The original matrix - output by rows: -4.00000000 14.00000000 0.00000000 -5.00000000 13.00000000 0.00000000 -1.00000000 0.00000000 2.00000000 The initial approximation vector is : 1.00000000 1.00000000 1.00000000

STZJbnB1dH50aGV+dG9sZXJhbmNlLnwrRzYi

SS1Ub2xlcmFuY2V+PX5HNiI=

JCIiIiEiJg==

SURJbnB1dH5tYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zfkc2Ig==

SSwtfmludGVnZXIufCtHNiI=

SUBNYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zfj1+RzYi

IiNE

NiRJJVF+PX5HNiIkIitMTExMaiEiKg==

SUpJbnB1dH5uZXd+UT9+RW50ZXJ+MX5mb3J+eWVzfm9yfjJ+Zm9yfm5vLkc2Ig==

SSpJbnB1dH5pc35HNiI=

IiIj

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5maWxlRzYi

STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig==

SSpJbnB1dH5pc35HNiI=

IiIi

INVERSE POWER METHOD Iteration Eigenvalue Eigenvector 1 6.18181818 1.00000000 0.72727273 -0.19580420 2 6.01724137 1.00000000 0.71551724 -0.24505203 3 6.00171526 1.00000000 0.71440823 -0.24952239 4 6.00017143 1.00000000 0.71429796 -0.24995339 5 6.00001714 1.00000000 0.71428694 -0.24999543 6 6.00000171 1.00000000 0.71428584 -0.24999955 Eigenvalue = 6.00000171 to tolerance = 1.0000000000e-05 obtained on iteration number = 6 Unit eigenvector is : 1.00000000 0.71428584 -0.24999955

JSFH