> restart:
> # QR ALGORITHM 9.6
> #
> # To obtain the eigenvalues of a symmetric, tridiagonal n by n matrix
> #
> #          a(1)   b(2)
> #          b(2)   a(2)   b(3)
> #             .      .      .
> #               .      .      .
> #                 .      .      .
> #                 b(n-1)  a(n-1)  b(n)
> #                            b(n)   a(n)
> #
> # INPUT:   n: A(1),...,A(n) (diagonal of A): B(2),...,B(n)
> #          (off-diagonal of A): maximum number of iterations M, 
> #          tolerance TOL.
> #
> # OUTPUT:  Eigenvalues of A or recommended splitting of A, or a 
> #          message that the maximum number of iterations was 
> #          exceeded.
> print(`This is the QR Method.`):
> 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 response is`): print(FLAG):
> if FLAG = 2 then
>    print(`The tridiagonal symmetric array A will be input from `):
>    print(`a text file in the order:`):
>    print(` (diagonal): A(1), A(2), ..., A(n),`):
>    print(` (subdiagonal): B(2), B(3), ..., B(n).`):
>    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 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(`Input file name is `): print(NAME):
>       INP := fopen(NAME,READ,TEXT):
>       OK := FALSE:
>       while OK = FALSE do
>          print(`Input the dimension n.`):
>          N := scanf(`%d`)[1]: print(`Dimension n = `): print(N):
>          if N > 1 then
>             for I1 from 1 to N do
>                A[I1-1] := fscanf(INP, `%f`)[1]:
>             od:
>             for I1 from 2 to N do
>                B[I1-1] := fscanf(INP, `%f`)[1]:
>             od:
>             OK := TRUE:
>             fclose(INP):
>          else print(`The number must be a positive integer.`):
>          fi:
>       od:  
>    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(`The tridiagonal symmetric array A will be input in the order `):
>          print(` (diagonal): A(1), A(2), ..., A(n),`):
>          print(` (subdiagonal): B(2), B(3), ..., B(n).`):          
>          for I1 from 1 to N do
>             print(`Input diagonal entry A(I) for I = `):  print(I1): 
>             A[I1-1] := scanf( `%f`)[1]: print(`Data is `): print(A[I1-1]):
>          od:
>          for I1 from 2 to N do
>             print(`Input subdiagonal entry B[I] for I = `): print(I1):
>             B[I1-1] := scanf(`%f`)[1]: print(`Data is `): print(B[I1-1]):
>          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 the tolerance.`):
>       TOL := scanf(`%f`)[1]: print(`Tolerance = `): print(TOL):
>       if TOL > 0 then
>          OK := TRUE:
>       else
>          print(`Tolerance must be a positive real number.`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input the maximum number of iterations.`):
>       L := scanf(`%d`)[1]: print(`Maximun number of iterations = `): print(L):
>       if L > 0 then
>          OK := TRUE:
>       else
>          print(`The number must be a positive integer.`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OUP := default:
>    fprintf(OUP, `The original matrix :\n`):
>       fprintf(OUP, `Diagonal:\n`):
>       for I1 from 1 to N do
>          fprintf(OUP, ` %11.8f`, A[I1-1]):
>       od:
>       fprintf(OUP, `\nOff diagonal:\n`):
>       for I1 from 2 to N do
>          fprintf(OUP, ` %11.8f`, B[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 response 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, `QR METHOD\n\n`):
> # STEP 1 
>    SHIFT := 0:
>    K := 1:
> # STEP 2 
>    while K <= L and OK = TRUE do
>       fprintf(OUP, `Iteration number %d N = %d\n`, K, N):
>       fprintf(OUP, `The array A is now as follows:\n`):
>       fprintf(OUP, `Diagonal:\n`):
>       for I1 from 1 to N do
>          fprintf(OUP, ` %11.8f`, A[I1-1]):
>       od:
>       fprintf(OUP, `\nOff diagonal:\n`):
>       for I1 from 2 to N do
>          fprintf(OUP, ` %11.8f`, B[I1-1]):
>       od:
>       fprintf(OUP, `\n \n`):
> # Steps 3-7 test for success  
> # STEP 3  
>       if abs(B[N-1]) <= TOL then
>          A[N-1] := A[N-1] + SHIFT:
>          fprintf(OUP, `Eigenvalue = %12.8f\n`, A[N-1]):
>          N := N-1:
>       fi:
> # STEP 4
>       if abs(B[1]) <= TOL then
>          A[0] := A[0]+SHIFT:
>          fprintf(OUP, `Eigenvalue = %12.8f\n`, A[0]):
>          N := N-1:
>          A[0] := A[1]:
>          for J from 2 to N do
>             A[J-1] := A[J]:
>             B[J-1] := B[J]:
>          od:
>       fi:
> # STEP 5  
>       if N = 0 then
>          OK := FALSE:
>       fi:
> # STEP 6
>       if N = 1 then
>          A[0] := A[0] + SHIFT:
>          fprintf(OUP,`Eigenvalue = %12.8f\n`, A[0]):
>          OK := FALSE:
>       fi:
> # STEP 7  
>       if OK = TRUE then
>          M := N-1:
>          if M >= 2 then
>             for I1 from 2 to M do
>                if abs(B[I1-1]) <= TOL then
>                   OK := FALSE:
>                   J := I1:
>                fi:
>             od:
>             if OK = FALSE then
>                fprintf(OUP, `Split the matrix into\n`):
>                for I1 from 1 to J-1 do
>                   fprintf(OUP,`%11.8f`,A[I1-1]):
>                od:
>                fprintf(OUP,`\n`):
>                for I1 from 2 to J-1 do
>                   fprintf(OUP,`%11.8f`,B[I1-1]):
>                od:
>                fprintf(OUP,`\n and \n`):
>                for I1 from J to N do
>                   fprintf(OUP,`%11.8f`,A[I1-1]):
>                od:
>                fprintf(OUP,`\n`):
>                for I1 from J+1 to N do
>                   fprintf(OUP,`%11.8f`,B[I1-1]):
>                od:
>                fprintf(OUP,`\n`):
>             fi:
>          fi:
>       fi:
>       if OK = TRUE then
> # STEP 8 
> # compute shift 
>          B1 := -(A[N-1]+A[N-2]):
>          C1 := A[N-1]*A[N-2]-B[N-1]*B[N-1]:
>          D1 := B1*B1-4*C1:
>          if D1 < 0 then
>             fprintf(OUP, `Problem with complex roots, D1 = %.8e\n`, D1):
>             OK := FALSE:
>          else
>             D1 := sqrt(D1):
> #  STEP 9
>             if B1 > 0 then
>                X1 := -2*C1/(B1+D1):
>                X2 := -(B1+D1)/2:
>             else
>                X1 := (D1-B1)/2:
>                X2 := 2*C1/(D1-B1):
>             fi:
> # if N = 2 then the 2 eigenvalues have been computed 
> #  STEP 10  
>             if N = 2 then
>                X1 := X1+SHIFT:
>                X2 := X2+SHIFT:
>                fprintf(OUP, `The last two eigenvalues are: %12.8f%11.8f\n`,X1, X2):
>                OK := FALSE:
>             else
> # STEP 11 
>                if abs(A[N-1]-X1) > abs(A[N-1]-X2) then
>                   X1 := X2:
>                fi:
> # STEP 12 
> # accumulate shift 
>                SHIFT := SHIFT+X1:
> # STEP 13
> # perform shift 
>                for I1 from 1 to N do
>                   DI[I1-1] := A[I1-1]-X1:
>                od:
> # STEP 14 and 15 compute R(K) 
> # STEP 14 
>                X[0] := DI[0]:
>                Y[0] := B[1]:
> # STEP 15 
>                for J from 2 to N do
>                   Z[J-2] := sqrt((X[J-2]*X[J-2])+(B[J-1]*B[J-1])):
>                   C[J-1] := X[J-2]/Z[J-2]:
>                   S[J-1] := B[J-1]/Z[J-2]:
>                   Q[J-2] := C[J-1]*Y[J-2]+S[J-1]*DI[J-1]:
>                   X[J-1] := C[J-1]*DI[J-1]-S[J-1]*Y[J-2]:
>                   if J <> N then
>                      R[J-2] := S[J-1]*B[J]:
>                      Y[J-1] := C[J-1]*B[J]:
>                   fi:
>                od:
>                M := N-1:
>                MM := N-2:
> # Steps 16-18 compute A(K+1) 
> # STEP 16 
>                Z[N-1] := X[N-1]:
>                A[0] := C[1]*Z[0]+S[1]*Q[0]:
>                B[1] := S[1]*Z[1]:
> # STEP 17 
>                if N > 2 then
>                   for J from 2 to M do
>                      A[J-1] := C[J]*C[J-1]*Z[J-1]+S[J]*Q[J-1]:
>                      B[J] := S[J]*Z[J]:
>                   od:
>                fi:
> # STEP 18 
>                A[N-1] := C[N-1]*Z[N-1]:
>             fi:
>          fi:
>       fi:
> # STEP 19 
>       K := K+1:
>    od:
> # STEP 20 
>    if OK = TRUE and K > L then
>       fprintf(OUP, `Maximum Number of Iterations Exceeded.\n`):
>    fi:
> # the process is complete 
>    if OUP <> default then
>       fclose(OUP):
>       print(`Output file `,NAME,` created successfully`):
>    fi:
> fi: