From:	IN%"joan@cs.albany.edu" 25-FEB-1993 15:57:17.43To:	IN%"VPAN@lcvax.lehman.cuny.edu"CC:	Subj:	Return-path: <joan@cs.albany.edu>Received: from vice1.lehman.cuny.edu by lcvax.lehman.cuny.edu (PMDF #12809) id <01GV57GYR9O0001XHQ@lcvax.lehman.cuny.edu>; Thu, 25 Feb 1993 15:57 ESTReceived: from cssun.albany.edu by vice1.lehman.cuny.edu (5.51/SMI-3.0DEV3) id AA23995; Thu, 25 Feb 93 15:52:27 ESTReceived: from steinmetz.albany.edu by cssun.albany.edu (5.65/SMI-3.2) id <AA23321@cssun.albany.edu>; Thu, 25 Feb 93 15:53:28 -0500Received: by steinmetz.albany.edu (5.65/ALBANY_CS-1.1) id AA27330; Thu, 25 Feb 93 15:53:14 -0500Date: Thu, 25 Feb 93 15:53:14 -0500From: joan@cs.albany.edu (Joan Bentley)To: VPAN@lcvax.lehman.cuny.eduMessage-id: <9302252053.AA27330@steinmetz.albany.edu>.rm CM.so refmac.AM.ll 6.5i.nr LL 6.5i.nr HM .5i.nr FM .5i.ps 11.nr PS 11.EQdelim @@gfont Rgsize 11.EN.DS C.B\s+2A New Algorithm for the Symmetric Eigenvalue Problem\s-2.R.DE.TSexpand, tab(#);c c.V. Pan#J. DemmelComputer Science Dept., SUNY#Computer Science Dept.Albany, NY 12222#Courant Instituteand#New York UniversityMathematics and Computer Science Dept.#251 Mercer St.Lehman College#New York, NY  10012CUNY, Bronx, NY 10468(Supported by NSF Grant CCR-8805782 and byPSC-CUNY Awards No. 668541 and 669210).TE.ls 2.SHSummary\u*)\d.PPWe apply a novel approach to approximate within @epsilon@ to all the.FS.sp .3m\u*)\d This paper is a revision of the earlier report.[p804.],presented at the SIAM Annual Meeting, San Diego, California, July 1989..FEeigenvalues of an @n times n@ symmetric tridiagonal matrix Ausing at most@n sup 2 (^\(lc^3 ~log sub 2 (625n sup 6 )^\(rc ~+~(83n-34)^\(lc^log sub 2 log sub 2 (( lambda sub 1 - lambda sub n )/^(2 epsilon ))/^log sub 2 (25n))~\(rc@ arithmetic operationswhere @lambda sub 1@ and @lambda sub n@ denote theextremal eigenvalues of A.  The algorithm can be modified to compute any fixednumbers of the largest and the smallesteigenvalues of A and may also be applied to the bandsymmetric matrices without their reduction to the tridiagonal form..SH1.  Introduction.PPIn this paper we present a new algorithm for the classical problem ofapproximate evaluation of the eigenvalues of an @n times n@symmetric tridiagonal matrix A.  Similarly to some other authors (see.[p418.],.[p691.]and.[p690.]),we rely on the bisection procedure, but we accelerate it by quite differentmeans, that is, by means of the techniques adapted from.[p447.],and we approximate within @epsilon@ to all the eigenvalues of Ausing at most@n sup 2 (^\(lc^3 ~log sub 2 (625n sup 6 )^\(rc ~+~(83n-34)^\(lc^log sub 2 log sub 2 (( lambda sub 1 - lambda sub n )/^(2 epsilon ))/^log sub 2 (25n))~\(rc@complex arithmetic operations (hereafter also referred toas to ops) where @lambda sub 1@ and @lambda sub n@ denote theextremal eigenvalues of A.  The algorithm may also be applied to the bandsymmetric matrices without their reduction to tridiagonal form..PPWe think that, besides the novel techniques used, theattractive feature of our algorithm is its convergence in the numbers ofsteps and arithmetic operations roughly proportional to@log sub 2 ~log sub 2 (1/^epsilon )@ rather thanto @log sub 2 (1/^epsilon ).@  The acceleration of thebisection procedure for computing the eigenvalues can betraced back to.[p890.],but the known algorithms rarely share the latter feature in thecases where several eigenvalues may form clusters.  In particular,clustered eigenvalues may slow down the convergence of QR-algorithmsand of the algorithm of.[p635.](see also.[p636.]),but the convergence of our algorithm is bothat least globally quadratic and insensitive to clustering theeigenvalues of the input matrix.  The bisection algorithmshares such an insensitivity to clustering the eigenvaluesbut has an inferior convergence rate; in particular,its convergence is much slower in the cases specifiedin the beginning of section 5.  On the other hand,our algorithm is never inferior to the bisection process,to which we may always shift back from our accelerated version if,say, we run into any numerical stability problems (see section 11 below)..PPThe algorithm is routinely extended to computing thesingular values of a bidiagonal matrix [see.[p036.],p. 427]..PPSomewhat similar techniques are used in.[p692.]in order to approximate to the absolutely largestand absolutely smallest eigenvalues of A..PPEven though we improved the classical bisection approach to thelevel of @O (log~log (1/ epsilon )),@ it is desirable to improveour construction itself, in particular, to decrease its startup cost of @n sup 2 \(lc ^3~ log sub 2 (625n sup 6 ) ^\(rc@ arithmeticoperations and to make our algorithm more robust.  Furthermore, ourauxiliary techniques of numerical integration may be eventually outperformedby some other methods, in particular, by computing Laguerre's intervals eachtime when the suspect intervals are bisected, as this was suggested by areferee.  (Such a computation outputs the number of the eigenvaluesin each suspect interval as a by-product.).PPWe will present our results in the following order.  In section 2we will recall some definitions and simple auxiliary facts.  Insection 3 we will present the simpler but slower version of ouralgorithm, with the computational cost estimates followed in section 4.  Insection 5 we will outline the acceleration procedure, leading to our mainalgorithm, specified in sections 6-8.  In section 9 wewill summarize our results on the complexityof our algorithm and will comment on its numericalstability.  In section 10 we will show an alternateversion of our algorithm, where we replace its blocks of factorization andinversion of tridiagonal matrices by the blocks for evaluating thecharacteristic polynomials of such matrices and of their leadingprincipal submatrices, as well as the derivatives ofsuch polynomials.  (We will assume that the matrices aretridiagonal where we need this.).PP\fB Acknowledgements. \fR  The authors thank the referees forhelpful comments and also Sally Goodall and Joan Bentley whoskillfully typed the manuscript..SH2.  Some Definitions and Auxiliary Facts.PPWe will recall some definitions and basic facts on matrix computations (see.[p036.]and.[p418.])..PPHereafter, @W~=~[w sub ij ]@ denotes a matrix; \fB v \fR a vector,@W sup T@ and @ bold v sup T@ their transposes; I the @n times n@ identitymatrix; @Trace~W~=~sum from i ~w sub ii ;@.EQ (2.1)|| W || sub h ~=~ max from {bold v} ~ || W bold v || sub h ~/~|| bold v || sub h.EN.LPwhere the maximum is over all the vectors @ bold v~!=~bold 0 ;@.EQ (2.2)|| W || sub 1 ~=~ || W sup T || sub inf ~=~ max from j ~ sum from i ~| w sub ij |..EN.PP.BProposition 2.1..R.IP a).IFor an @n times n@ symmetric matrix @A@, all its eigenvalues@lambda sub 1 , ..., lambda sub n@ are real;.IP b).I@|| A || sub 1 ~/~ sqrt n ~<=~|| A || sub 2 ~=~ || A sup T || sub 2 ~=~max from {1 <= j <= n } ~| lambda sub j |~<=~|| A || sub 1 ;@.IP c).Ithe matrices @(A~-~ mu I ) sup s@ have the eigenvalues@( lambda sub j - mu ) sup s@ for @j=1, ...,n@ and for all @mu@ and @s@,as long as @(A- mu I) sup s@ is defined;.IP d).I@Trace~A~=~ sum from j=1 to n lambda sub j .@.SH3.  Bisection-Exclusion Algorithm for the Eigenvalues of a Symmetric Matrix.PPFor an @n times n@ symmetric matrix A, we will seek all itseigenvalues lying between two fixed realscalars @lambda sup -@ and @lambda sup + ,~ lambda sup - ~<~ lambda sup + ,@provided that we are also given a black box subroutineperforming the following task:.LP.BSubroutine COUNT EIGENVALUES @( mu )@.LP\fB Input: \fR  a real @mu@ and a real symmetric matrix A..LP\fB Output: \fR  the numbers of the eigenvalues of A that are less than @mu ,@equal @mu@ and exceed @mu .@.PPHere and hereafter, the eigenvalues are always counted withtheir algebraic multiplicities..PPWe will compute the desired eigenvalues of A by means of the bisectionalgorithm, which may also be viewed as Weyl's exclusion algorithm.[(p190.])for computing the real zeros of the characteristic polynomial@p( lambda )~=~det( lambda I-A)@ (without computing its coefficients):.LP\fB Algorithm 3.1 \fR (compare Remarks 3.1 and 4.1 below)..LP\fB Input: \fR  an @n times n@ symmetric matrix A, two realscalars @lambda sup - @ and @lambda sup + ,~ lambda sup - ~<~ lambda sup + ,@the two numbers @ k sup - ~and~k sup + @ of the eigenvalues of @A@ that areless than @ lambda sup - @ and that are not lessthan @ lambda sup + ,@ respectively, and a positive @epsilon .@.LP\fB Output: \fR  approximations, within absolute errors@epsilon ,@ to all the @n-k sup - - k sup +@ eigenvalues of A inthe interval @"{" lambda : lambda sup - ~<=~ lambda ~<~ lambda sup + "}".@.LP\fB Initialize: \fR  The real semiclosed interval@"{" lambda :~ lambda sup - ~<=~ lambda ~<~ lambda sup + "}"@is said to be \fI suspect at the first step \fR of computations.  Callthe subroutine COUNT EIGENVALUES@( mu )@ for@mu ~=~ lambda sup -@ and @mu ~=~ lambda sup + .@.LP\fB Computations (recursive step): \fR  For every semiclosed suspect interval@"{" lambda :~ italic l ~<=~ lambda ~<~r"}",@call the subroutine COUNT EIGENVALUES(@mu@) for@mu ~=~( italic l +r)/^2.@  Count the number of the eigenvalues of A ineach of the two semiclosed intervals:  @"{" lambda :~ italic l ~<=~lambda~<~( italic l +r)/^2"}"@and @"{" lambda :~( italic l +r)/^2~<=~lambda~<~r"}"@ and call such aninterval \fI suspect at the next recursive step \fR unless no eigenvalues ofA lie in this interval..LP\fB Stopping criterion: \fR  do not apply recursive steps tothe suspect intervals that have lengths less than@2 epsilon ,@ but for each such an interval, output its midpointand the number of the eigenvalues of A lying in it..PP\fB Remark 3.1. \fR  To compute all the eigenvalues of @A,@ we shallapply Algorithm 3.1 to @ lambda sup - ~and~ lambda sup + @ such that.EQ (3.1)lambda sup - ~<=~ lambda sub n ~<=~"..."~<=~ lambda sub 1 ~<~lambda sup +.ENIf we only need to compute the g largest and the h smallesteigenvalues of A, we may apply a simplified version ofAlgorithm 3.1 (we will call it \fB Restricted Algorithm 3.1 \fR),where in each recursive step we discard all the suspect intervalsnot containing the desired eigenvalues of A and only work withthe remaining (at most g+h) suspect intervals..PPThe subroutine COUNT EIGENVALUES(@mu@) can be reduced tocomputing the factorization.EQ (3.2)A~-~mu I~=~L( mu ) D( mu )L sup T ( mu ).EN.LPwhere @L( mu )@ is a triangular matrix and @D( mu )@ is adiagonal matrix (see section 10 for an alternate computation)..PPThe factorization (3.2) is computed by usual algorithms.[(p036.],.[p418.],.[p488.]).Symmetric pivoting may improve numerical stability (see.[p675.]),unless A is tridiagonal, in which case pivoting is not needed (see.[p488.],pp. 302-303).  The preliminaryreduction of A to the tridiagonal form is recommended,as usual, in order to decrease the computational cost of thefactorization, except for the special cases where the factorization isinexpensive even without such a reduction (compare.[p418.],pp. 119, 122,.[p203.],.[p204.],.[p640.])..PP\fB Remark 3.2 \fR (compare section 10).  Hereafter, we will simplify ourstudy by always assuming that we deal with @mu@ and Athat permit such a factorization.  For some unfortunatechoices of @ mu ,@ the factorization (3.2) may failto exist even with the introduction of pivoting, and insuch cases we would have to shift @ mu ,@ which usuallygives satisfactory results according to.[p418.],Section 3.3.  In practice, one must guard againstoverflow and division by zero in computing the entries of@D( mu )@ as this is done in the EISPACKroutine based on [EISPACK] or as is described in.[p896.].  This technique is straightforwardand maintains the excellent numerical stability properties of thealgorithm..PPTo complete the evaluation of the eigenvalues of A, here is anexpression for the scalars @lambda sup -@ and @lambda sup +@satisfying (3.1) (where we assume that A is tridiagonal):.PP.BProposition 3.1..R.ILet @A@ be a symmetric tridiagonal matrix, @t~=~(1/n)Trace~A,@ @B~=~A-tI,@@b~=~|| B || sub 1 ,@ @delta@ be a positive scalar,@lambda sup - ~=~t-b,@ @lambda sup + ~=~t+b+ delta .@  Then.EQ (3.3)lambda sub 1 - lambda sub n ~<=~ lambda sup + - lambda sup - ~<=~ 2sqrt 3 ~ ( lambda sub 1 - lambda sub n )+ delta ..EN.PP\fB Proof. \fR  Deduce from Proposition 2.1 and from (2.1) and (2.2) that@|| B || sub 2 ~=~ max "{" lambda sub 1 -t,@@t- lambda sub n "}"~<=~|| B || sub 1 ~<=~sqrt 3~ || B || sub 2 .@  Combiningthese relations implies (3.3).\(sq.SH4.  Computational Cost of the Bisection-Exclusion Algorithm.PPEstimating the computational cost of Algorithm 3.1, we will use thefollowing simple result:.PP.BProposition 4.1..R.IThere are at most @n@ suspect intervals at each recursive step ofAlgorithm 3.1.  They have the same length, @L@; the midpoint of each suchan interval lies at the distance at most @L/^2@ from an eigenvalueof @A@, and every eigenvalue of @A@ lies at the distance at most @L/^2@from the midpoint of an interval suspect at this recursive step..R.PPCombining Propositions 3.1 and 4.1 and the observation that the lengthof a suspect interval is halved in each recursive step, we deduce that at most.EQn~\(lc log(2 sqrt 3 ( lambda sub 1 ~-~ lambda sub n )/^ (2 epsilon ))^\(rc ~=~ n~\(lc~log sub 2 ( sqrt 3 ( lambda sub 1 - lambda sub n )/^epsilon )~\(rc.EN.LPcalls to the subroutine COUNT EIGENVALUES(@mu@) are needed in theentire Algorithm 3.1..PPNow, we may estimate the cost of computing the eigenvalues of A.  Wewill assume that A is already in the tridiagonal formand will use the bound (3.3) and the known estimate of @3n-3@ ops forcomputing the factorization (3.2) (see.[p418.]).We arrive at the following results:.PP.BTheorem 4.1..R.IUnder the assumption of Remark 3.2, Algorithm 3.1 approximates within@epsilon@ to all the @k@ eigenvalues of an @n times n@symmetric tridiagonal matrix @A@ lying in the interval@"{" lambda :~ lambda sup - ~<=~ lambda ~<~ lambda sup + "}" @ byusing @3nk ^ \(lc log sub 2 ( sqrt 3 ( lambda sup + ~-~ lambda sup - )/^(2 epsilon ))^ \(rc @arithmetic operations,which means @3n sup 2 ~\(lc~log sub 2 ( sqrt 3 ( lambda sub 1 - lambda sub n )/^epsilon )~\(rc@under (3.1) and (3.3).  If only the @g@largest and the @h@ smallest eigenvalues of @A@ need to be approximated, then @3n(g+h)~\(lc~log sub 2 ( sqrt 3 ( lambda sub 1 - lambda sub n )/^epsilon )~\(rc~)@arithmetic operations suffice in the Restricted Algorithm 3.1..R.PP\fB Remark 4.1. \fR  Having computed a realinterval containing a single eigenvalue of A,Algorithm 3.1 continues recursive bisection steps to narrow such an interval..[p418.]recommends shifting to the secant algorithm at this point, but Newton'siteration @mu (k+1)~=~mu (k)-c( mu (k))/^c'( mu (k))@converges to a zero of @c ( mu )@ even faster than the secant algorithmdoes.  We may compute the ratio by using the well-known equation:.EQ-c( mu )/^c'( mu )~=~1/^Trace(A- mu I) sup -1.ENwhere @c ( mu ) ~=~det(A ~-~ mu I ).@  Thecost of this computation is 10n ops for a tridiagonal matrix A (see.[p673.])..SH5.  Acceleration of the Eigenvalue Computation.  An Outline of the Algorithm.PPIn this section we will accelerate our approximation to all the eigenvalues ofA by incorporating some techniques from.[p447.]into Algorithm 3.1.  To see a potential for an improvement of Algorithm 3.1,consider the case where all the eigenvalues of an @n times n@ input matrix Alie in the two semiclosed intervals @"{" lambda :~0~<=~ lambda ~<~ epsilon "}"@ and@"{" lambda :~1- epsilon ~<=~ lambda ~<~1"}".@  In this case, Algorithm 3.1requires about @log sub 2 (1/ epsilon )@ recursive steps, but the solutioncould be obtained immediately if we recognized that the interval@"{" lambda :~epsilon~<=~lambda~<~1- epsilon "}"@ contains no eigenvalues ofA.  We will improve Algorithm 3.1 in this case and, moreover, in the caseof a general symmetric matrix, decreasing the estimate@3n sup 2 ~\(lc~log sub 2 ( sqrt 3 ( lambda sub 1 - lambda sub n )/^epsilon )~\(rc ~)@of Theorem 4.1 to@n sup 2 (~\(lc~3~log sub 2 (625n sup 6 )~\(rc~+~ (83n-34)^\(lc~log sub 2 (log sub 2 (( lambda sub 1 -lambda sub n )/^(2 epsilon )) /^log sub 2 (25n))~\(rc ~).@  Towardsthis goal, we will now recall some concepts from.[p447.],where these concepts have been defined for the complex discs ratherthan for the real intervals..PP\fB Definition 5.1. \fR  Given a symmetric matrix A, two intervals ofthe real line are \fI equivalent \fR to each other if they containexactly the same eigenvalues of A..PP\fB Definition 5.2. \fR  For a real interval J, |J| denotes its length;furthermore, its \fI rigidity ratio \fR,r.r.(J), and its \fI isolation ratio \fR, i.r.(J), are defined as follows:.EQ Ir.r.(J)~=~ "inf" |J^ sup - | / | J | ,.EN.EQ Ii.r.(J)~=~ "sup" |J^ sup + |/ | J |..EN.LPwhere the infimum and the supremum are over all the real intervals@J sup - ~!subset~ J@ and @J sup + ~!supset~ J@ equivalent to Jprovided that the real intervals J and @J sup +@ areconcentric, that is, have the same midpoint.  For example, if@A~=~diag "{" 1,10,100,1000 "}",~J~=~ [10, 101],@ thenr.r (J) = 90/91, i.r. (J) = 54.5/45.5 = 109/ 91..PPNext, we will accelerate convergence of Algorithm 3.1 to the eigenvalues ofA lying in a fixed interval J, under the following assumption:.PP\fB Assumption 5.1. \fR  For threefixed real numbers \fI l \fR, r and @nu@ such that.EQ (5.1)r~>~ italic l ,~~~nu sup {^2} ~>=~ 8 sqrt n ,.EN.LPand for a positive integer @k~<=~n,@ the interval@J~=~"{" lambda :~ italic l ~<=~lambda~<~r"}"@ contains exactlyk eigenvalues of A, and furthermore,.EQ (5.2)i.r.(J)~>=~ nu sup {^2} ..EN.PPSpecifically, under the latter assumption, we will compute a subinterval@J sub 1 sup *@ of J, equivalent to J and, furthermore, such that.EQ (5.3)r.r.(J sub 1 sup * )~>=~1/~ nu sup {^2}.EN.LPunless.EQ (5.4)|J sub 1 sup * | ~<=~2 epsilon ..EN.PPThe midpoint of @J sub 1 sup *@ already approximates within @epsilon@ to allthe k eigenvalues of A lying in J if (5.4) holds.  Otherwise (5.3) holds,and then the maximum distance between two successive eigenvalues of A on@J sub 1 sup *@ is not less than @|J sub 1 sup * |/^( nu sup {^2} (k-1)).@  In this case,@\(lc~log sub 2 (k sup 2 ~nu sup {^4} )~\(rc@ recursive steps of Algorithm3.1 applied to the interval @J sub 1 sup *@ will surely suffice in order toyield at least two suspect intervals lying in @J sub 1 sup *@ and havingisolation ratios at least @nu sup {^2} .@  (This follows becausein every step of Algorithm 3.1 the length of suspect intervals ishalved.)  The same argument can be recursively applied to the new suspectintervals, until we compute the desired approximations to all theeigenvalues of A in J within @epsilon .@  We will start with the interval@J~=~"{" lambda :~lambda sup - ~=~italic l ~<=~lambda~<~lambda sup + ~=~r"}",@satisfying Assumption 5.1 for k=n and for any @nu @ satisfying (5.1), and wemay need at most @n-1@ such recursive subdivisions of J into suspect intervals(because every suspect interval contains at least one eigenvalueof A).  Next, let us specify the algorithm that computes a desired subinterval@J sub 1 sup *@ for a given interval J satisfying Assumption 5.1..LP.BAlgorithm 5.1..R.LP\fB Input: \fR  an @n times n@ symmetric matrix A,a positive @epsilon ,@ and the parameters of Assumption 5.1..LP\fB Output: \fR  real @m sup *@ and nonnegative @d sup *@ that define a realsubinterval@J sub 1 sup * ~=~"{" lambda :~m sup * -d sup * ~<=~lambda~<~m sup * +d sup * "}"@of J containing k eigenvalues of A and satisfying (5.3) and/or (5.4)..LP\fB Computations \fR (also see section 7):.IP 1.Compute an approximation @m sup *@ to m, the average of all the keigenvalues of A in J..IP 2.Compute an approximation @d sup *@ to the minimum d such that all thek eigenvalues of A in J also lie in both intervals@J sub 1 ~=~"{"lambda :~m-d~<=~lambda~<~m+d"}"@ and @J sub 1 sup * .@.PPWe observe that @r.r.(J sub 1 )~>=~n/^(2n-2)~>~\(12@ unless@|J sub 1 |~=~0.@  If @m sup *@ is close to m and @d sup *@ is close to d,then the interval @J sub 1 sup * ~>=~J sub 1@ closely approximates to@J sub 1 ,@ and we may expect that @r.r.(J sub 1 sup * )@ also exceeds \(12 or at least exceeds @1/ ^nu sup {^2}@ unless @|J sub 1 sup * |@ isvery small..PPIn the next sections, we will specify stages 1 and 2 of Algorithm 5.1and will formalize the latter argument..SH6.  An Auxiliary Result.PPIn this section we will recall the following auxiliary result from.[p120.]and.[p640.]:.PP.BProposition 6.1..R.ILet an @n times n@ matrix @B@ have eigenvalues @mu sub 1 , ..., mu sub n@(counted with their algebraic multiplicity) such that.EQ (6.1)inf~=~| mu sub 0 |~>~| mu sub 1 |~>=~ | mu sub 2 |~>=~"..."~>=~| mu sub n-k |~>~h~>~| mu sub n-k+1 |~>=~"..."~>=~| mu sub n |~>~0,.EN.LP\fI for some positive h.  Let \fR.EQ (6.2)s sub p ~=~ sum from j=n-k+1 to n ~ mu sub j sup p ~for~p~>=~0;.EN.EQ (6.3)s sub p ~=~- sum from j=1 to n-k ~ mu sub j sup p ~for~p~<~0;.EN.EQ (6.4)s sub p (q)~=~(1/q)~ sum from i=0 to q-1 ~(h~ omega sup i ) sup p+1 Trace(h ~omega sup i I-B) sup -1.EN.LP.Iwhere @omega@ is a primitive @q-th@ root of1 (that is, @omega sup q ~=~1,@ @omega sup i ~!=~1@ for @0~<~i~<~q),@and an integer @q@ exceeds @|p|.@  Finally, let.EQ (6.5)g~=~max from {1 <=~ j ~<= n} ~min~ "{"| mu sub j |~/~h,~h~/~| mu sub j |"}"..EN.LP.I Then.EQ (6.6)DELTA sub p sup * (q)~=~|s sub p (q)~-~s sub p |~<=~ DELTA sub p sup + (q) ~=~(k~g sup q+p +(n-k)g sup q-p )h sup p /(1-g sup q )..EN.PP\fB Proof. \fR  With no loss of generality, let @h~=~1.@  Let @| lambda |~=~1@and consider the Laurent expansion:.EQ (6.7)Trace~( lambda I-B) sup -1 ~=~sum from j=1 to n ~1/( lambda - lambda sub j )~=~sum from {p= - inf} to inf ~s sub p lambda sup -p-1 ~=~sum from {i= - inf} to inf ~c sub i lambda sup i ,.EN.LPwhere the second equation follows from the relations (6.2), (6.3) and thedecompositions.EQ1 over {lambda - lambda sub j} ~=~ 1 over lambda ~sum from i=0 to inf ~ left ( {lambda sub j} over lambda right ) sup {^i} ~if~j~>~n-k,.EN.EQ1 over {lambda - lambda sub j} ~=~ -~ 1 over {lambda sub j} ~sum from i=0 to inf ~ left ( {lambda sub j} over lambda right ) sup -i ~if~j~<=~n-k..EN.LPCombine the equations (6.4) and (6.7) and deduce that.EQs sub p (q)~=~ sum from {i= - inf} to inf ~c sub -p-1+iq ..EN.LPOn the other hand, the equations (6.7) imply that@s sub p+j ~=~c sub -p-1-j@ for all integers j.  Therefore,.EQ|s sub p (q)-s sub p |~=~| sum from {i= - inf} to -1 ~s sub p-iq ~+~sum from i=1 to inf ~s sub p-iq |~<=~| sum from i=1 to inf ~s sub p+iq |~+~| sum from i=1 to inf ~s sub p-iq | ~..EN.EQ|s sub p+iq |~<=~ sum from j=n-k+1 to n ~| lambda sub j | sup p+iq ~<=~k| lambda sub n-k+1 | sup p+iq ~<=~kg sup p+iq ,.EN.LPso that.EQ| sum from i=1 to inf ~s sub p+iq |~<=~k~sum from i=1 to inf ~g sup p+iq ~<=~kg sup q+p /(1-g sup q )..EN.LPSimilarly,.EQ| sum from i=1 to inf ~s sub p-iq |~<=~(n-k)~sum from i=1 to inf ~g sup p-iq ~<=~ (n-k)g sup q-p /(1-g sup q ),.EN.LPand Proposition 6.1 follows..SH7.  Contraction of the Interval Containing the Eigenvalues.PPNext, we will apply Proposition 6.1 to approximateto the average m of the k eigenvalues of A in J underAssumption 5.1.  We will initially set.EQ (7.1)B~=~B sub 0 ~=~A-m sub 0 I,~2m sub 0 ~=~ italic l +r ,~h~=~h sub 0 ~=~nu ~(r- italic l )/^2..EN.LP[so that @g~<=~1/^nu@ due to (5.2)], choose a sufficiently largeq and compute the values @s sub 1 (q),~s sub 2 (q)@ and.EQ (7.2)m(q)~=~s sub 1 (q)/^k~+~( italic l +r )/^2.EN.LPm(q) approximates.EQ (7.3)m~=~(s sub 1 /^k)~+~( italic l +r)/^2,.EN.LPthe average value of the k eigenvalues of A in J [see (6.1)-(6.6)]..PPNext, we observe that for the matrix B redefined as@B~=~A-mI,@ the value @sqrt {s sub 2}@ lies within a factor of@sqrt k@ from d, the maximum distance between m and an eigenvalueof A in J, specifically,.EQ (7.4)sqrt {s sub 2} ~/~ sqrt k ~<=~d~<=~ sqrt {s sub 2} ,.EN.LPso the pair @m sup * ~=~m@ and @d sup * ~=~ sqrt {s sub 2}@ satisfies therequirements to the output of Algorithm 5.1.  Actually, we will use theapproximations m(q) to m, @s sub 2 (q)@ to @s sub 2@ and the auxiliary value.EQ (7.5)d(q)~=~ sqrt {s sub 2 (q)+ alpha (q)} ,~~alpha (q)~=~1.001 DELTA sub 2 sup + (q),.EN.LPwhere @DELTA sub 2 sup + (q)@ is defined by (6.6).  Thelatter definition ensures that the interval@J sub 1 (q)~=~"{" lambda :~m(q)-d(q)~<=~lambda~<~m(q)+d(q)"}"@contains all the k eigenvalues of A lying in J..PP\fB Remark 7.1. \fR  In fact,we may replace 1.001 in (7.5) by any constant exceeding 1..PPApplying the bound (6.6) for p=1 and p=2, as well as therelations (7.2)-(7.5), we deduce that the values m(q) and d(q)converge to m and @d sup * ,@ respectively, as q and/or@1/^g @ grow to @inf .@  This could suggest choosing a larger q and setting.EQ (7.6)m sup * ~=~m(q), ~d sup * ~=~d(q),~J sub 1 sup * ~=~J sub 1 (q).EN.PPHowever, the cost of computing @s sub p (q)@ of (6.4) for p=1 and p=2(which are the most costly stages of our computations) includes (andactually is dominated by) the cost of computing@Trace(h omega sup i I-B) sup -1@ for @i=0,1, ...,q-1,@which is proportional to qn and is too large for large q.  Actually, thiscost is at most @8q(n-1)@ complex ops reduced to at most@20q(n-1)@ real ops for all tridiagonal matrices @h omega sup i I-B@ ifwe apply the algorithm of section 10 below..PP\fB Assumption 7.1. \fR  Hereafter, we willassume for our operation count that the cost of computing@Trace(h omega sup i I-B) sup -1@ is at most @8(n-1)@ complex ops(reduced to at most 20n real ops) for all h, @omega ,@ i and Bthat we will encounter (note that @ omega @ is a complex number)..PPThus, instead of increasing q as much as possible,we will let q take a relatively small value,.EQ (7.7)q~=~8~+~\(lc~log sub 2 (4.03n/^log sub 2 ~nu )~\(rc ,.EN.LPwhich means that @q~=~9@ if @nu~>=~4.03 n.@  [Section 8 willexplain the choice of the factor 4.03, which is due to thefactors 1.001 in (7.6) and 0.501 in (8.2) below and which canbe decreased to any constant exceeding 4 as @n ~->~ inf @ (see Remarks7.1 and 8.1)].  Some increase of @q@ and, therefore, of the complexityof the computations can be justified if we need to decrease the approximationerror bound @ DELTA sub p sup + (q)@ of (6.6)..PPIf we start with @m sub 0 ,~ h sub 0 ,~ and~ B sub 0@ of (7.1) anddefine q by (7.7), the resulting values of (7.6) do not generallysatisfy the requirements to the output of Algorithm 5.1.  We may,however, update @m sub 0 ,~h sub 0 ,~and~B sub 0@ of (7.1) by setting@m sub {0,new} ~=~m(q),~h sub {0,new} ~=~h(q),~B sub {0,new} ~=~B (q) ~=~A~-~m(q)~I,@ thenrespectively update m(q), h(q) and B(q) and recursively repeat thisprocess until (7.6) finally furnishes us with the desired outputof Algorithm 5.1..PPThis recursive process is formally defined by (7.8)-(7.11) below.  The keyobservation is that in each recursive step i [for an appropriate choice of@h~=~h(i)]@ the length of the interval @J sub 1 (q)@ decreases bythe factor of at least @ nu sup {2 sup i },@ as long as all the keigenvalues of A in @J sub 1 (q)@ actually lie in a much smallersubinterval of @J sub 1 (q).@  When the latter requirement is no longersatisfied, we arrive at (5.3) and/or (5.4) and stop the recursive process..PPLet us now give its formal description.  We will initially set@m sub old ~=~m sub 0 ,@ @h sub old ~=~h sub 0@ [see (7.1)]and will recursively recompute the matrix B and the values h,@s sub 1 (q),@ @s sub 2 (q),@ m(q) and d(q) according to the followingformulae:.EQ (7.8)B sub {new} ~=A-m sub {old} (q)I;.EN.EQ (7.9)h sub new sup 2 ~=~d sub old (q)min~"{"m sub old (q)-m sub 0 +h sub 0 nu ,~m sub 0 +h sub 0 nu -m sub old (q)"}";.EN.EQ (7.10)m sub new (q)~=~m sub old (q)+s sub 1 (q)/^k;.EN.EQ (7.11)d sub new (q)~=~ sqrt {s sub 2 (q)+ alpha (q)}.EN.LPfor @m sub 0@ and @h sub 0@ of (7.1) and for @s sub 1 (q),@ @s sub 2 (q)@and @alpha (q)@ defined by (6.4) and (7.6) for @B~=~B sub new @and @h~=~h sub new .@  (Note that exactly k eigenvalues of A lie in theinterval @"{" lambda :~ m sub 0 -h sub 0 nu ~<=~ lambda ~<~ m sub 0 +h sub 0 nu "}".)@.PPThus, we complete Algorithm 5.1 as follows:.LP\fB Algorithm 5.1 \fR (see input and output in section 5)..LP\fB Initialization: \fR  @h sub 0 ~=~ nu (r- italic l )/^2,~q@ is defined by(7.7), @ omega ~=~ exp ~(2 pi ~ sqrt -1 /^q),@ @h sub {old} ~=~h sub 0 ,~m sub {old} ~=~( italic l +r)/^2,~d sub {old} ~=~(r- italic l )/^2.@.LP\fB Computations \fR (CONTRACT and UPDATE):  Successivelycompute @g sub {new} ~=~d sub {old} /^h sub {old} ,~ DELTA sub {new} sup + ~=~ {kg sub {new} sup {q+2} ~+~(n-k)~g sub {new} sup {q-2}} over {1~-~g sub {new} sup q} @[compare (6.6)], @ alpha (q) ~=~1.001 DELTA sub {new} sup + ,@.LP@h sub {new} sup 2 @ from (7.9), @h sub {new} ~=~ sqrt {(h sub {new} sup 2 )} ,~m sub {new} (q)@ from (7.10), @B sub {new} @ from (7.8), @s sub 1 (q)~and~s sub 2 (q)@ from (6.4) [whereonly the real part of the right side of (6.4) is needed since@s sub 1 (q)~and~s sub 2 (q)@ are real] and @d sub {new} (q)@ from (7.11)..RSIf @d sub {new} ~<~ epsilon ,@ print @m sub {new} @ and stop..brElse if @d sub {new} /^d sub {old} ~>~1/^ nu sup 2 @ [then (5.3) holds], set@lambda sup - ~=~m sub {new} ~-~d sub {new} ,~ lambda sup + ~=~m sub {new} ~+~ d sub {new} @ and call Algorithm 3.1..RSElse replace @h sub {old} ,~m sub {old} ,~d sub {old} @ by @h sub {new} ,~m sub {new} ,~d sub {new} ,@ respectively, and repeat computations..RE.RE.PPThe correctness of the algorithm [in particular, the implication of (5.3) bythe inequality @d sub {new} /^d sub {old} ~>~ 1/^ nu sup 2 ]@ follows fromPropositions 7.1-7.3 below..PP\fB Proposition 7.1. \fR  @\(lc~log sub 2 (log(|J|/^(2 epsilon ))/^log~nu ~ sup 2 )~\(rc@.Irecursivesteps of the procedure (7.8)-(7.11) suffice to yield approximations@m sup * ~=~m(q)@ to @m@ and @d sup * ~=~d(q)@ to @d@ satisfying therequirements (5.3) and/or (5.4) for @J sub 1 sup * ~=~J sub 1 (q).@.R.PPWe will start our proof of Proposition 7.1 with the following simple proposition:.PP.BProposition 7.2..R.IThe inequality (5.3) is ensured for @J sub 1 sup * ~=~J sub 1 (q)@\fR [ \fIthat is, @r.r.(J sub 1 (q))~>=~1/^nu sup {^2} ]@ if@1~-~ {DELTA sub q (J sub 1 )} over 2d(q) ~>=~2 over {nu sup {^2}}@ where@DELTA sub q (J sub 1 )~=~|J sub 1 (q)|-|J sub 1 |.@.R.PP\fB Proof. \fR  Recall that @r.r.(J sub 1 )~>~\(12 @ and deduce that.EQr.r.(J sub 1 (q))~=~r.r.(J sub 1 )~{|J sub 1 |} over {|J sub 1 (q)|} ~>=~{|J sub 1 |} over {2|J sub 1 (q)|} ~>=~{|J sub 1 (q) |- DELTA sub q (J sub 1 )} over {2|J sub 1 (q)|} ~=.EN.EQ(2d(q)~-~DELTA sub q (J sub 1 ))/^(4d(q)),~since~ |J sub 1 (q)|~=~2d(q)..EN.LPThis implies Proposition 7.2..ti +5.5i\(sq.ti 0.PPIt remains to consider the case where, in every recursive step (7.7)-(7.10),@|J sub 1 (q)|~>~2 epsilon @ and.EQ (7.12)DELTA sub q (J sub 1 )/^(2d(q))~>~1~-~2/^nu sup {^2} ..EN.LPIn the next section, we will prove the following result:.PP.BProposition 7.3..R.IThe inequality (7.12) implies that.EQ|J sub 1 (q)|~<=~|J|/^nu sup {^2} ..EN.LP.Iunless @r.r.(J sub 1 (q))~>=~1/^nu sup {^2} ,@in which case (5.3) holds for @J sub 1 sup * ~=~J sub 1 (q).@.PPConsequently, since @J sub 1 (q)@ is a subinterval of J,.EQi.r.(J sub 1 (q))~>=~i.r.(J) nu sup {^2} ~>=~ nu sup 4 ..EN.PPRecursively applying the same construction with the isolationratios growing to exceed @nu sup {^4} , nu sup {^8}, nu sup {^16} , "...",@we will yield a sequence of intervals equivalent to J, having isolationratios at least @nu sup 4 ,@ @nu ~ sup 8 ,@ @nu sup 16 ,"..."@and lengths at most @|J|/^ nu sup 4 ,@@|J|/^nu ~ sup 8 ,@ @|J|/^ nu sup 16 , "...",@and this immediately implies Proposition 7.1..SH8.  Proof of Proposition 7.3..PPObserve that.EQ (8.1)DELTA sub q (J sub 1 )~<=~2( DELTA sub q (m)~+~DELTA sub q (d)).EN.LPwhere @DELTA sub q (m)@ and @DELTA sub q (d)@ denote the absolute values ofthe absolute errors of the approximation to m by m(q) and to d by d(q),respectively.  In our case, @nu ~>~3@ for @n~>~1@ [see (5.1)],@h~=~|J| nu /^2,@ @g~<=~1/^nu@ in (6.6), and we may assume that@q~>=~9@ [see (7.7)] and deduce that@g~<=~1/^3,@ @1-g sup q ~>=~1-3 sup -9 .@  Then (6.6)with @DELTA sub 1 sup * (q)/^k~=~DELTA sub q (m)@ implies that surely.EQ (8.2)DELTA sub q (m)/^|J|~<~(0.501~n/^k) nu sup {^2-q} ..EN.PP\fB Remark 8.1. \fR  The factor 0.501 can be replaced by @0.5~+~ delta (n)@where @delta (n) ~->~0@ as @n ~->~ inf @ [compare (5.1) and (6.5)]..LPWe also deduce that.EQDELTA sub q (d)~=~d(q)-d~=~(d(q)~-~ size 10 {sqrt {s sub 2}} )~+~( size 10 {sqrt {s sub 2}} -d)~<=.EN.EQ{d(q) sup 2 -s sub 2} over {d(q)+ size 10 {sqrt {s sub 2}} } ~+~( sqrt k -1) d.EN.LPsince @d sup 2 ~<=~s sub 2 ~<=~kd sup 2 .@  Furthermore,@d(q) sup 2 -s sub 2 ~=~s sub 2 (q)-s sub 2 + alpha (q)~<=~2.001 DELTA sub 2 sup * (q)@[see (6.6) and (7.6)].  Therefore,@DELTA sub q (d)~<=~2.001 DELTA sub 2 sup * (q)/^(d(q)~+~ sqrt {s sub 2} )~+~( sqrt k -1 )d.@  Applying (6.6), we deducethat @DELTA sub 2 sup * (q)~<=~n nu sup 4-q |J| sup 2 /^(4(1-3 sup -9 )),@and consequently,.EQ (8.3)DELTA sub q (d)/^|J|~<=~[0.5005 ~n nu sup 4-q |J|/^d(q)]~+~( sqrt k -1 )d/|J|..EN.LPCombining (8.1)-(8.3) implies that.EQDELTA sub q (J sub 1 )/^|J|~<=~2[ DELTA sub q (m)~+~DELTA sub q (d)]/^|J|~<=.EN.EQn nu sup 4-q [(1.001 |J|/^d(q))~+~1.002/^(k nu sup {^2} )]~+~2( sqrt k -1 )d/^ |J|..EN.LPRecall (7.12) and deduce that.EQ1 ~-~ 2 over nu sup {^2} ~<~ {DELTA sub q (J sub 1 )} over {|J sub 1 (q)|} ~=~{DELTA sub q (J sub 1 )} over {|J|} ~~ |J| over {|J sub 1 (q)|} ~<=.EN.EQ|J| over {|J sub 1 (q)|} ~~n nu sup 4-q ~left ( {1.001|J|} over d(q) ~+~ 1.002 over {k nu sup {^2}} right ) ~+~{2( sqrt k -1) d} over {|J sub 1 (q)|} ..EN.LPNote that @2d(q)~=~|J sub 1 (q)|@ and deduce that.EQf( nu ,k,q)~=~1-2/^ nu sup {^2} -2( sqrt k -1 )d/^|J sub 1 (q)|~<=.EN.EQ (8.4)|J| over {|J sub 1 (q)|} ~n nu sup 4-q ~ left ({2.002 |J|} over {|J sub 1 (q)|} ~+~ 1.002 over {k nu sup {^2}} right ) ..EN.LPWe may assume that@|J|/^|J sub 1 (q)|~<~ nu sup {^2}@ (for otherwise Proposition 7.3would follow) and that @k~>=~1;@also, we recall that @nu ~>=~3@ and deduce from (5.1) and (8.4) that.EQ (8.5)f(v, k, q)~<=~2.015 n nu ~ sup 6-q |J|/^|J sub 1 (q)|..EN.LPOn the other hand,.EQ (8.6)d/^|J sub 1 (q)|~<=~2/^nu sup {^2} ,.EN.LPfor otherwise@r.r.(J sub 1 (q))~=~r.r.(J sub 1 )d/^|J sub 1 (q)|~>=~1/ nu sup {^2} ,@and then Proposition 7.1 would followfor @J sub 1 sup * ~=~J sub 1 (q)@.  Combine (8.4) and (8.6) and deduce that.EQf( nu ,k,q)~>=~1-2/^ nu sup {^2} -4( sqrt k -1)/ ^ nu sup {^2} ~=~1-4 sqrt k /^ nu sup {^2} +2/^ nu sup {^2} ..EN.LPBy Assumption 5.1, @nu sup {^2} ~>=~8 sqrt n ~>=~8 sqrt k @[see (5.1)], so that.EQ (8.7)f( nu ,k,q)~>=~\(12 ..EN.LPCombining the inequalities (8.5) and (8.7), we obtain that.EQ|J|/^|J sub 1 (q)|~>=~ nu sup q-6 ~/~ (4.03n)..EN.LPCombining the latter inequality and (7.7) implies Proposition 7.3..SH9.  Computational Cost of the Accelerated Algorithm and Its Numerical Stabilization.PPLet us summarize our computations and estimates..PP\fB Algorithm 9.1. \fR.PP\fB Input: \fR  the interval@J sub 1 sup * ~=~"{" lambda :~ lambda sup - ~<=~ lambda~<~lambda sup + "}"@defined in Proposition 3.1 and satisfying (5.3)..PP\fB Output: \fR  the same as for Algorithm 3.1..PP\fB Computations: \fR.RS.IP 1.Fix @nu@ satisfying (5.1) by letting @nu~=~5n@ (to be certain)..IP 2.Apply Algorithm 3.1 until, in at most @\(lc~log sub 2 (n sup 2 nu sup 4 )~\(rc~=~\(lc~log sub 2 (625n sup 6 )~\(rc@recursive steps, at least two suspect intervals on @J sub 1 sup *@ withisolation ratios @nu ~ sup 2@ or greater have been defined..IP 3.Apply Algorithm 5.1 to each ofthese intervals until (5.3) and/or (5.4) have been satisfied..IP 4.If this process ends up with (5.3) [rather than with (5.4)], applyAlgorithm 3.1, and in at most @\(lc~log sub 2 (625n sup 6 )~\(rc@its steps (each involving @3n-3@ arithmetic operations),either satisfy (5.4) and print an approximation to aneigenvalue of @A@ within @ epsilon @ or else compute two or moredisjoint subintervals of the output interval ofAlgorithms 5.1, all these subintervals satisfyingAssumption 5.1.  Then recursively repeat applying Algorithm 5.1(step 3), each application followed by at most@\(lc~log sub 2 (n sup 2 nu sup 4 )~\(rc@ applications ofAlgorithm 3.1, until you approximate within @epsilon@ toall the eigenvalues of A..PPAt most @n-1@ applications of Algorithm 5.1 are needed, for in eachapplication, we strictly increase the number of disjoint intervals satisfyingAssumption 5.1, and this number is at most n, the number of theeigenvalues of A.  Each application of Algorithm 5.1 requires at most@\(lc~log sub 2 (log sub 2 (( lambda sub 1 - lambda sub n )/^(2 epsilon ))/^log sub 2 ~nu ~ sup 2 )~\(rc~@recursive steps, each step involving at most@9qn~+~2n-6q~+~20@ complex arithmetic operations, which is@83n~-~34@ if q = 9.  Summarizing, we arrive at the following result:.PP\fB Theorem 9.1 \fR [also see (10.1) below]..IUnder the assumption of Remark 3.2 and Assumption 7.1, all theeigenvalues of an @n times n@ symmetric tridiagonal matrix@A@ can be approximated within @epsilon@ by using at most@n sup 2 (83n-34)~\(lc~log sub 2 (log sub 2 (( lambda sub 1 - lambda sub n )/^(2 epsilon ))/~log sub 2 (25n))~\(rc~+~3 ^\(lc~log sub 2 (625n sup 6 )~\(rc~)@complex arithmetic operations.  This bound decreasesby @n/(g+h)@ times if we only need to compute the@g@ smallest and the @h@ largest eigenvalues of @A@..R.SH10.  Alternate Computations Based on Computing the Values of theCharacteristic Polynomials..PPIn this section we will recall how to compute thevalues @c sub i ( lambda ),~i=1, ...,n,@ and @c sub n sup ' ( lambda ),@where @c sub i ( lambda )@ is the characteristic polynomial of the@i times i @ leading principal submatrix of a given symmetric tridiagonalmatrix @A~=~[a sub ij ],~a sub ii ~=~a sub i ,~a sub i,i+1 ~=~ a sub i+1,i ~=~b sub i ,~i=1, ...,n-1.@.PPIt is well-known.[(p418.],section 7-10) that the number of sign agreements between the consecutive termsin the sequence @c sub 0 ( mu )~=~1,~c sub 1 ( mu ), ...,~c sub n ( mu )@equals the number of eigenvalues of @A@ which are less and @ mu ,@ and thisgives us an alternate way for performing the subroutine COUNT EIGENVALUES@( mu )@ of section 3, that is, we may compute @c sub i ( mu )@ for @i=1, ...,n@instead of computing the factorization (3.2)..PPFurthermore, @-p' ( lambda )/^p ( lambda ) ~=~Trace~(A~-~ lambda I) sup -1 ,@and this gives us an alternate way forcomputing @Trace ~(h omega sup i I~-~B) sup -1 @ of (6.4)..PPNow, we recall.[(p418.],section 3-4) that.EQc sub 0 ( lambda )~=~1,~c sub 1 ( lambda )~=~a sub 1 ~-~ lambda ,~c sub i+1 ( lambda )~=~(a sub i+1 ~-~ lambda )~c sub i ( lambda )~-~b sub i sup 2 ~c sub i-1 ( lambda ),~i=1, ...,n-1 ,.ENand thus we may compute the values @c sub j ( lambda )@ at @ lambda ~=~ mu @for @j=1,2, ...,n@ by using less than @4(n-1)@ arithmetic operations,which is less than we estimated for the subroutine COUNT EIGENVALUES @( mu ).@.PPFurthermore, as this has been observed in.[p882.],performing the same computation with @ lambda @ replaced by @ lambda ~+~ mu @for the variable @ lambda @ and for a fixed constant @ mu @and reducing all the polynomials modulo @ lambda sup 2 ,@ wearrive at the linearpolynomial @c sub n ( mu )~+~ lambda c sub n sup ' ( mu ),@ whichgives us @c sub 1 ( mu ), ...,c sub n ( mu )@and @c sub n sup ' ( mu )@ for the overall cost of less than@8(n-1)@ arithmetic operations, estimated in Assumption 7.1.  Some operandsinvolved in half of these operations @(a sub i+1 ~and~b sub i sup 2 )@ arereal, so that the overall cost of computing @Trace~(A~-~ lambda I) sup -1 @is less than @20(n-1)@ real arithmetic operations provided that eachcomplex addition/subtraction reduces to two real ops and each complexmultiplication reduces to 6 real ops..SH11.  Some Comments on Numerical Stability.PPWe have already pointed out that Algorithm 3.1, as a bisection algorithmfor the symmetric eigenvalue problem, is known to benumerically stable.  The only potential source of numericalinstability of Algorithm 5.1 could be the application of (6.4),and more specifically, the instability is a dangeronly to the extent to which the inversion of the shiftedinput matrices @C(q)~=~h(q)~ omega sup i I~-~B(q)@ is lessnumerically stable than their factorization.  This depends,of course, on the spectral radius of @(C(q)) sup -1 .@  On the other hand,.EQh(q) ~>=~ nu ~ | J | /^2.ENprovided that we seek the eigenvalues in the interval @J@ havingisolation ratio @ nu sup 2 ,@ and we fix a natural i and ensurethat in our computations, @ nu ~>=~(64n) sup 1/4 2 sup i@ [until this issatisfied, we will stay with the (unaccelerated) bisection algorithm(Algorithm 3.1)]; increasing @i@ by @i(0)@ means that we requireby @i(0)@ more bisection steps and thus ensure by @2 sup i(0) @ timesgreater an isolation ratio of the interval to which we apply theaccelerated procedure for the search for theeigenvalues.  Even if @ | J | ~=~2 epsilon ,@ the spectralradius of @C(q) sup -1 @ is at most@ rho (C(q) sup -1 )~<=~1/^( nu epsilon )~<=~1/^(2 sup i/2 (64n) sup 1/4 epsilon )@and, to ensure that@ rho (C(q) sup -1 )~<=~ rho (0)@ for a fixed @ rho (0),@ we onlyneed @i~=~ \(lc 3+2~log(n sup 1/4 rho sub 0 epsilon )^\(rc@ additionalbisection steps.  If this means a larger value @i,@ then we willjust stay longer with the (unaccelerated) bisection algorithm,but clearly, this approach does accelerate the bisectionalgorithm for a large subclass of symmetric tridiagonal matrices..PPWe note that the numerical stability of our algorithm depends onthe isolation ratios [which bound the spectral radii @ rho (C(q) sup -1 )@of the matrices @C(q) sup -1 ]@ but for a fixed @ rho (C(q) sup -1 )@does not depend on the multiplicity of or the clustering ofthe eigenvalues of @B.@.PPFinally, regarding the alternate approach of section 10, it is knownthat the evaluation of the characteristic polynomial @c ( lambda )@of a symmetric tridiagonal matrix does not usually lead to anysubstantial problems of numerical stability, and this shouldalso apply to the evaluation of @ c' ( lambda ).@.bp.[$LIST$.]