> 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\n \n`, 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:\\name.ext`):
>       print(`for example:   A:\\DATA.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:\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 approximate eigenvector:\n`):
>       for I1 from 1 to N do
>          fprintf(OUP, ` %11.8f`, V[I1-1]):
>       od:
>       fprintf(OUP, `\n`):
>       fprintf(OUP, `The approximate eigenvalue:\n`):
>       fprintf(OUP, ` %11.8f \n`, XMU):
>       fprintf(OUP, `The initial vector of dimension n-1:\n`):
>       for I1 from 1 to N-1 do
>          fprintf(OUP, ` %11.8f`, X[I1-1]):
>       od:
>       fprintf(OUP, `\n`):
>    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:\\name.ext`):
>          print(`for example   A:\\OUTPUT.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\n\n`):
> # 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:\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 I1 from 1 to N do
>             fprintf(OUP,` %11.8f`, VV[I1-1]):
>          od:
>          fprintf(OUP, `\n`):
>       fi:
>       if OUP <> default then
>          fclose(OUP):
>          print(`Output file `,NAME,` created successfully`):
>       fi:
>    fi:
> fi: