> 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.
> print(`This is the Symmetric 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.\n`):
>    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
>             Y[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
>          print(`Input Matrix `):
>          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):
>             Y[I1-1] := scanf(`%f`)[1]:print(`Data is `):print(Y[I1-1]):
>          od:
>          OK := TRUE:
>       else 
>          print(`The number must be a positive integer.\n`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    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
>    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`, Y[I1-1]):
>    od:
>    fprintf(OUP, `\n`):
> fi:
> if OK = TRUE then 
> for I1 from 1 to N do
> X[I1-1] := 0:
> od:
> 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`):
> print(`for example   A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]:print(`Output file is `):print(NAME):
> 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 I1 from 1 to N do
> XL := XL+Y[I1-1]*Y[I1-1]:
> od:
> # 2-Norm of Y
> XL := sqrt(XL):
> ERR := 0:
> if XL > 0 then
> for I1 from 1 to N do
> T := Y[I1-1]/XL:
> ERR := ERR+(X[I1-1]-T)*(X[I1-1]-T):
> X[I1-1] := T:
> od:
> # X has a 2-Norm of 1.0
> ERR := sqrt(ERR):
> else
> print(`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 I1 from 1 to N do
> Y[I1-1] := 0:
> for J from 1 to N do
> Y[I1-1] := Y[I1-1]+A[I1-1,J-1]*X[J-1]:
> od:
> YMU := YMU+X[I1-1]*Y[I1-1]:
> od:
> # Steps 5 and 6
> XL := 0:
> for I1 from 1 to N do
> XL := XL+Y[I1-1]*Y[I1-1]:
> od:
> # 2-Norm of Y
> XL := sqrt(XL):
> ERR := 0:
> if XL > 0 then
> for I1 from 1 to N do
> T := Y[I1-1]/XL:
> ERR := ERR+(X[I1-1]-T)*(X[I1-1]-T):
> X[I1-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 I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I1-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 I1 from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I1-1]):
> od:
> fprintf(OUP, `\n`):
> OK := FALSE:
> else
> # Step 8
> K := K+1:
> fi:
> fi:
> od:
> # Step 9
> if K > NN then
> print(`No convergence within `,NN,` iterations\n`):
> fi:
> fi:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: