/* WR-251 */ /* EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 251-252. */ /* c is the diagonal, b the sub-diagonal and beta the squared sub-diagonal of a symmetric tridiagonal matrix of order n. The eigenvalues lambda[m1], ..., lambda[m2], where m2 is not less than m1 and lambda[i+1] is not less than lambda[i], are calculated by the method of bisection and stored in the vector x. Bisection is continued until the upper and lower bounds for an eigenvalue differ by less than eps1, unless at some earlier stage, the upper and lower bounds differ only in the least significant digits. eps2 gives an extreme upper bound for the error in any eigenvalue, but for certain types of matrices the small eigenvalues are determined to a very much higher accuracy. In this case, eps1 should be set equal to the error to be tolerated in the smallest eigenvalue. It must not be set equal to zero. */ /* Input: c = an n x 1 array giving the diagonal elements of a tridiagonal matrix. b = an n x 1 array giving the sub-diagonal elements. b[1] may be arbitrary but is replaced by zero in the procedure. beta = an n x 1 array giving the squares of the sub-diaonal elements. beta[1] may be arbitrary but is replaced by zero in the procedure. Both b[i] and beta[i] are given since the squares may be the primary data. If storage economy is important then one of these can be dispensed with in an obvious way. n = the order of the tri-diagonal matrix. m1, m2 = the eigenvalues labda[m1], lambda[m1+1], ..., lambda[m2] are calculated (lambda[1] is the smallest eigenvalue). m1 <= m2 must hold otherwise no eigenvalues are computed. eps1 = a quantity affecting the precision to which the eigenvaues are computed. relfeh = the smallest number for which 1 + relfeh > 1 on the computer. [in PL/I, this is EPSILON]. Output: eps2 gives information concerning the accuracy of the results. z = total number of bisections to find all required eigenvalues. x = array x[m1:m2] contains the computed eigenvalues. */ BISECT:PROC(C,B,BETA,N,M1,M2,EPSA,RELFEH,EPS2,Z,X); /*PAGE 251*/ /*00000010*/ DCL (N,M1,M2,Z,I,A,K,LBN,UBN,LBM,UBM) BIN FIXED (31,0), /*00000020*/ EPSA BIN FLOAT(21), /*00000030*/ (EPS1,EPS2,RELFEH,H,XMIN,XMAX,Q,X1,XU,X0) BIN FLOAT(21), /*00000040*/ ((C,B,BETA)(*), /*00000050*/ (X(*),WU(LBM:UBM) CTL)) BIN FLOAT(21); /*00000060*/ LBN=LBOUND(C,1); UBN=LBN+N-1; /*00000070*/ LBM=LBOUND(X,1); UBM=LBM+M2-M1; /*00000080*/ ALLOC WU; /*00000090*/ /*CALCULATION OF XMIN,XMAX */ /*00000100*/ EPS1=EPSA; /*00000110*/ BETA(LBN),B(LBN)=0E0; XMIN=C(UBN)-ABS(B(UBN)); XMAX=C(UBN)+ABS(B(UBN));/*00000120*/ DO I=UBN-1 TO LBN BY -1; /*00000130*/ H=ABS(B(I))+ABS(B(I+1)); /*00000140*/ IF C(I)+H>XMAX THEN XMAX=C(I)+H; /*00000150*/ IF C(I)-H0E0 THEN EPS2=XMAX; ELSE EPS2=-XMIN; /*00000180*/ EPS2=RELFEH*EPS2; /*00000190*/ IF EPS1<=0E0 THEN EPS1=EPS2; /*00000200*/ EPS2=.5E0*EPS1+7E0*EPS2; /*00000210*/ /* INNER BLOCK */ /*00000220*/ X0=XMAX; /*00000230*/ DO I=LBM TO UBM; /*00000240*/ X(I)=XMAX; WU(I)=XMIN; /*00000250*/ END; /*00000260*/ Z=0E0; /*00000270*/ /*LOOP FOR KTH EIGENVALUE */ /*00000280*/ DO K=UBM TO LBM BY -1; /*00000290*/ XU=XMIN; /*00000300*/ DO I=K TO LBM BY -1; /*00000310*/ IF XUX(K) THEN X0=X(K); /*00000370*/ X1=(XU+X0)/2E0; /*00000380*/ DO WHILE(X0-XU>2E0*RELFEH*(ABS(XU)+ABS(X0))+EPS1); /*00000390*/ Z=Z+1; /*00000400*/ /*STURMS SEQUENCE */ /*00000410*/ A=LBM-M1; Q=1E0; /*00000420*/ DO I=LBN TO UBN; /*00000430*/ IF Q^=0E0 THEN Q=BETA(I)/Q; /*00000440*/ ELSE Q=ABS(B(I))/RELFEH; /*00000450*/ Q=C(I)-X1-Q; /*00000460*/ IF Q<0E0 THEN A=A+1; /*00000470*/ END; /*00000480*/ IF AX1 THEN X(A)=X1; /*00000520*/ END; /*00000530*/ ELSE X0=X1; /*00000540*/ X1=(XU+X0)/2E0; /*00000550*/ END; /*00000560*/ X(K)=(X0+XU)/2E0; /*00000570*/ END; /*00000580*/ FREE WU; /*00000590*/ END BISECT; /*00000600*/