> 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:
> 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:
> 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.\n`):
> 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)\n`):
>    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.\n`):
>    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:\\name.ext`):
>       print(`for example:   A:\\DATA.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:
>             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.\n`):
>          fi:
>       od:
>    else 
>       print(`The program will end so the input file can be created.\n`):
>    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:
>          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.\n`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OUP := default:
>    fprintf(OUP, `The original matrix - output by rows:\n`):
>    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, `\n`):
>    od:
>    fprintf(OUP, `The initial approximation vector is :\n`):
>    for I1 from 1 to N do
>       fprintf(OUP, ` %11.8f`, X[I1-1]):
>    od:
>    fprintf(OUP, `\n`):
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input the tolerance.\n`):
>       TOL := scanf(`%f`)[1]:print(`Tolerance = `):print(TOL):
>       if TOL > 0 then
>          OK := TRUE:
>       else
>          print(`Tolerance must be positive number.\n`):
>       fi:
>    od:
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input maximum number of iterations `):
>       print(`- integer.\n`):
>       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.\n`):
>       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:
> 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:\\name.ext\n`):
> print(`for example   A:\\OUTPUT.DTA\n`):
> NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `INVERSE POWER METHOD\n\n`):
> fprintf(OUP, `Iteration  Eigenvalue  Eigenvector\n`):
> # 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\n`, K, YMU):
> for I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> if ERR < TOL then 
> OK := FALSE:
> fprintf(OUP, `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`):
> for I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> else
> # Step 12
> K := K+1:
> fi:
> od:
> if K > NN then
> fprintf(OUP,  `No convergence in %d iterations\n`,NN):
> fi:
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: