alg096.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg096 := proc() local OK, AA, NAME, INP, N, I, A, B, TOL, L, FLAG, OUP, SHIFT, K, J, M, B1, C1, D1, X1, X2, D, X, Y, Z, C, S, Q, R, MM;
> printf(`This is the QR Method.\n`);
> OK := FALSE;
> printf(`The tridiagonal symmetric array A will be input from `);
> printf(`a text file in the order:\n`);
> printf(` (diagonal): A(1), A(2), ..., A(n),\n`);
> printf(` (subdiagonal): B(2), B(3), ..., B(n).\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> AA := scanf(`%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:  A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the dimension n.\n`);
> N := scanf(`%d`)[1];
> if N > 1 then
> for I from 1 to N do
> A[I-1] := fscanf(INP, `%f`)[1];
> od;
> for I from 2 to N do
> B[I-1] := fscanf(INP, `%f`)[1];
> od;
> OK := TRUE;
> else
> printf(`Dimension must be greater then 1.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be a positive real number.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of iterations.\n`);
> L := scanf(`%d`)[1];
> if L > 0 then
> OK := TRUE;
> else
> printf(`The number must be a positive integer.\n`);
> fi;
> od;
> fclose(INP);
> else
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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 I from 1 to N do
> fprintf(OUP, ` %11.8f`, A[I-1]);
> od;
> fprintf(OUP, `\nOff diagonal:\n`);
> for I from 2 to N do
> fprintf(OUP, ` %11.8f`, B[I-1]);
> od;
> fprintf(OUP, `\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 I from 2 to M do
> if abs(B[I-1]) <= TOL then
> OK := FALSE;
> J := I;
> fi;
> od;
> if OK = FALSE then
> fprintf(OUP, `Split the matrix into\n`);
> for I from 1 to J-1 do
> fprintf(OUP,`%11.8f`,A[I-1]);
> od;
> fprintf(OUP,`\n`);
> for I from 2 to J-1 do
> fprintf(OUP,`%11.8f`,B[I-1]);
> od;
> fprintf(OUP,`\n and \n`);
> for I from J to N do
> fprintf(OUP,`%11.8f`,A[I-1]);
> od;
> fprintf(OUP,`\n`);
> for I from J+1 to N do
> fprintf(OUP,`%11.8f`,B[I-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 I from 1 to N do
> D[I-1] := A[I-1]-X1;
> od;
> # STEP 14 and 15 compute R(K) 
> # STEP 14 
> X[0] := D[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]*D[J-1];
> X[J-1] := C[J-1]*D[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):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg096();