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:\134\134name.ext`): print(`for example: A:\134\134DATA.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 :\134n`): fprintf(OUP, `Diagonal:\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, A[I1-1]): od: fprintf(OUP, `\134nOff diagonal:\134n`): for I1 from 2 to N do fprintf(OUP, ` %11.8f`, B[I1-1]): od: fprintf(OUP, `\134n`):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:\134\134name.ext`): print(`for example A:\134\134OUTPUT.DTA`): NAME := scanf(`%s`)[1]: print(`Output file is `):print(NAME): OUP := fopen(NAME,WRITE,TEXT): else OUP := default: fi: fprintf(OUP, `QR METHOD\134n\134n`):# STEP 1 SHIFT := 0: K := 1:# STEP 2 while K <= L and OK = TRUE do fprintf(OUP, `Iteration number %d N = %d\134n`, K, N): fprintf(OUP, `The array A is now as follows:\134n`): fprintf(OUP, `Diagonal:\134n`): for I1 from 1 to N do fprintf(OUP, ` %11.8f`, A[I1-1]): od: fprintf(OUP, `\134nOff diagonal:\134n`): for I1 from 2 to N do fprintf(OUP, ` %11.8f`, B[I1-1]): od: fprintf(OUP, `\134n \134n`):# 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\134n`, 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\134n`, 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\134n`, 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\134n`): for I1 from 1 to J-1 do fprintf(OUP,`%11.8f`,A[I1-1]): od: fprintf(OUP,`\134n`): for I1 from 2 to J-1 do fprintf(OUP,`%11.8f`,B[I1-1]): od: fprintf(OUP,`\134n and \134n`): for I1 from J to N do fprintf(OUP,`%11.8f`,A[I1-1]): od: fprintf(OUP,`\134n`): for I1 from J+1 to N do fprintf(OUP,`%11.8f`,B[I1-1]): od: fprintf(OUP,`\134n`): 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\134n`, 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\134n`,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.\134n`): fi:# the process is complete if OUP <> default then fclose(OUP): print(`Output file `,NAME,` created successfully`): fi:fi:JSFH