1      SUBROUTINE DLARRB( N, D, L, LD, LLD, IFIRST, ILAST, SIGMA, RELTOL,
2     $                   W, WGAP, WERR, WORK, IWORK, INFO )
3*
4*  -- LAPACK auxiliary routine (instru to count ops, version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     June 30, 1999
8*
9*     .. Scalar Arguments ..
10      INTEGER            IFIRST, ILAST, INFO, N
11      DOUBLE PRECISION   RELTOL, SIGMA
12*     ..
13*     .. Array Arguments ..
14      INTEGER            IWORK( * )
15      DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ), W( * ),
16     $                   WERR( * ), WGAP( * ), WORK( * )
17*     ..
18*     Common block to return operation count
19*     .. Common blocks ..
20      COMMON             / LATIME / OPS, ITCNT
21*     ..
22*     .. Scalars in Common ..
23      DOUBLE PRECISION   ITCNT, OPS
24*     ..
25*
26*  Purpose
27*  =======
28*
29*  Given the relatively robust representation(RRR) L D L^T, DLARRB
30*  does ``limited'' bisection to locate the eigenvalues of L D L^T,
31*  W( IFIRST ) thru' W( ILAST ), to more accuracy. Intervals
32*  [left, right] are maintained by storing their mid-points and
33*  semi-widths in the arrays W and WERR respectively.
34*
35*  Arguments
36*  =========
37*
38*  N       (input) INTEGER
39*          The order of the matrix.
40*
41*  D       (input) DOUBLE PRECISION array, dimension (N)
42*          The n diagonal elements of the diagonal matrix D.
43*
44*  L       (input) DOUBLE PRECISION array, dimension (N-1)
45*          The n-1 subdiagonal elements of the unit bidiagonal matrix L.
46*
47*  LD      (input) DOUBLE PRECISION array, dimension (N-1)
48*          The n-1 elements L(i)*D(i).
49*
50*  LLD     (input) DOUBLE PRECISION array, dimension (N-1)
51*          The n-1 elements L(i)*L(i)*D(i).
52*
53*  IFIRST  (input) INTEGER
54*          The index of the first eigenvalue in the cluster.
55*
56*  ILAST   (input) INTEGER
57*          The index of the last eigenvalue in the cluster.
58*
59*  SIGMA   (input) DOUBLE PRECISION
60*          The shift used to form L D L^T (see DLARRF).
61*
62*  RELTOL  (input) DOUBLE PRECISION
63*          The relative tolerance.
64*
65*  W       (input/output) DOUBLE PRECISION array, dimension (N)
66*          On input, W( IFIRST ) thru' W( ILAST ) are estimates of the
67*          corresponding eigenvalues of L D L^T.
68*          On output, these estimates are ``refined''.
69*
70*  WGAP    (input/output) DOUBLE PRECISION array, dimension (N)
71*          The gaps between the eigenvalues of L D L^T. Very small
72*          gaps are changed on output.
73*
74*  WERR    (input/output) DOUBLE PRECISION array, dimension (N)
75*          On input, WERR( IFIRST ) thru' WERR( ILAST ) are the errors
76*          in the estimates W( IFIRST ) thru' W( ILAST ).
77*          On output, these are the ``refined'' errors.
78*
79*****Reminder to Inder --- WORK is never used in this subroutine *****
80*  WORK    (input) DOUBLE PRECISION array, dimension (???)
81*          Workspace.
82*
83*  IWORK   (input) INTEGER array, dimension (2*N)
84*          Workspace.
85*
86*****Reminder to Inder --- INFO is never set in this subroutine ******
87*  INFO    (output) INTEGER
88*          Error flag.
89*
90*  Further Details
91*  ===============
92*
93*  Based on contributions by
94*     Inderjit Dhillon, IBM Almaden, USA
95*     Osni Marques, LBNL/NERSC, USA
96*
97*  =====================================================================
98*
99*     .. Parameters ..
100      DOUBLE PRECISION   ZERO, TWO, HALF
101      PARAMETER          ( ZERO = 0.0D0, TWO = 2.0D0, HALF = 0.5D0 )
102*     ..
103*     .. Local Scalars ..
104      INTEGER            CNT, I, I1, I2, INITI1, INITI2, J, K, NCNVRG,
105     $                   NEIG, NINT, NRIGHT, OLNINT
106      DOUBLE PRECISION   DELTA, EPS, GAP, LEFT, MID, PERT, RIGHT, S,
107     $                   THRESH, TMP, WIDTH
108*     ..
109*     .. External Functions ..
110      DOUBLE PRECISION   DLAMCH
111      EXTERNAL           DLAMCH
112*     ..
113*     .. Intrinsic Functions ..
114      INTRINSIC          ABS, DBLE, MAX, MIN
115*     ..
116*     .. Executable Statements ..
117*
118      EPS = DLAMCH( 'Precision' )
119      I1 = IFIRST
120      I2 = IFIRST
121      NEIG = ILAST - IFIRST + 1
122      NCNVRG = 0
123      THRESH = RELTOL
124      DO 10 I = IFIRST, ILAST
125         OPS = OPS + DBLE( 3 )
126         IWORK( I ) = 0
127         PERT = EPS*( ABS( SIGMA )+ABS( W( I ) ) )
128         WERR( I ) = WERR( I ) + PERT
129         IF( WGAP( I ).LT.PERT )
130     $      WGAP( I ) = PERT
131   10 CONTINUE
132      DO 20 I = I1, ILAST
133         IF( I.EQ.1 ) THEN
134            GAP = WGAP( I )
135         ELSE IF( I.EQ.N ) THEN
136            GAP = WGAP( I-1 )
137         ELSE
138            GAP = MIN( WGAP( I-1 ), WGAP( I ) )
139         END IF
140         OPS = OPS + DBLE( 1 )
141         IF( WERR( I ).LT.THRESH*GAP ) THEN
142            NCNVRG = NCNVRG + 1
143            IWORK( I ) = 1
144            IF( I1.EQ.I )
145     $         I1 = I1 + 1
146         ELSE
147            I2 = I
148         END IF
149   20 CONTINUE
150*
151*     Initialize the unconverged intervals.
152*
153      I = I1
154      NINT = 0
155      RIGHT = ZERO
156   30 CONTINUE
157      IF( I.LE.I2 ) THEN
158         IF( IWORK( I ).EQ.0 ) THEN
159            DELTA = EPS
160            OPS = OPS + DBLE( 1 )
161            LEFT = W( I ) - WERR( I )
162*
163*           Do while( CNT(LEFT).GT.I-1 )
164*
165   40       CONTINUE
166            IF( I.GT.I1 .AND. LEFT.LE.RIGHT ) THEN
167               LEFT = RIGHT
168               CNT = I - 1
169            ELSE
170               S = -LEFT
171               CNT = 0
172               DO 50 J = 1, N - 1
173                  OPS = OPS + DBLE( 5 )
174                  TMP = D( J ) + S
175                  S = S*( LD( J ) / TMP )*L( J ) - LEFT
176                  IF( TMP.LT.ZERO )
177     $               CNT = CNT + 1
178   50          CONTINUE
179               TMP = D( N ) + S
180               IF( TMP.LT.ZERO )
181     $            CNT = CNT + 1
182               IF( CNT.GT.I-1 ) THEN
183                  OPS = OPS + DBLE( 4 )
184                  DELTA = TWO*DELTA
185                  LEFT = LEFT - ( ABS( SIGMA )+ABS( LEFT ) )*DELTA
186                  GO TO 40
187               END IF
188            END IF
189            OPS = OPS + DBLE( 1 )
190            DELTA = EPS
191            RIGHT = W( I ) + WERR( I )
192*
193*           Do while( CNT(RIGHT).LT.I )
194*
195   60       CONTINUE
196            S = -RIGHT
197            CNT = 0
198            OPS = OPS + DBLE( 5*( N-1 )+1 )
199            DO 70 J = 1, N - 1
200               TMP = D( J ) + S
201               S = S*( LD( J ) / TMP )*L( J ) - RIGHT
202               IF( TMP.LT.ZERO )
203     $            CNT = CNT + 1
204   70       CONTINUE
205            TMP = D( N ) + S
206            IF( TMP.LT.ZERO )
207     $         CNT = CNT + 1
208            IF( CNT.LT.I ) THEN
209               OPS = OPS + DBLE( 4 )
210               DELTA = TWO*DELTA
211               RIGHT = RIGHT + ( ABS( SIGMA )+ABS( RIGHT ) )*DELTA
212               GO TO 60
213            END IF
214            WERR( I ) = LEFT
215            W( I ) = RIGHT
216            IWORK( N+I ) = CNT
217            NINT = NINT + 1
218            I = CNT + 1
219         ELSE
220            I = I + 1
221         END IF
222         GO TO 30
223      END IF
224*
225*     While( NCNVRG.LT.NEIG )
226*
227      INITI1 = I1
228      INITI2 = I2
229   80 CONTINUE
230      IF( NCNVRG.LT.NEIG ) THEN
231         OLNINT = NINT
232         I = I1
233         DO 100 K = 1, OLNINT
234            NRIGHT = IWORK( N+I )
235            IF( IWORK( I ).EQ.0 ) THEN
236               OPS = OPS + DBLE( 2 )
237               MID = HALF*( WERR( I )+W( I ) )
238               S = -MID
239               CNT = 0
240               OPS = OPS + DBLE( 5*( N-1 )+1 )
241               DO 90 J = 1, N - 1
242                  TMP = D( J ) + S
243                  S = S*( LD( J ) / TMP )*L( J ) - MID
244                  IF( TMP.LT.ZERO )
245     $               CNT = CNT + 1
246   90          CONTINUE
247               TMP = D( N ) + S
248               IF( TMP.LT.ZERO )
249     $            CNT = CNT + 1
250               CNT = MAX( I-1, MIN( NRIGHT, CNT ) )
251               IF( I.EQ.NRIGHT ) THEN
252                  IF( I.EQ.IFIRST ) THEN
253                     OPS = OPS + DBLE( 1 )
254                     GAP = WERR( I+1 ) - W( I )
255                  ELSE IF( I.EQ.ILAST ) THEN
256                     OPS = OPS + DBLE( 1 )
257                     GAP = WERR( I ) - W( I-1 )
258                  ELSE
259                     OPS = OPS + DBLE( 2 )
260                     GAP = MIN( WERR( I+1 )-W( I ), WERR( I )-W( I-1 ) )
261                  END IF
262                  OPS = OPS + DBLE( 2 )
263                  WIDTH = W( I ) - MID
264                  IF( WIDTH.LT.THRESH*GAP ) THEN
265                     NCNVRG = NCNVRG + 1
266                     IWORK( I ) = 1
267                     IF( I1.EQ.I ) THEN
268                        I1 = I1 + 1
269                        NINT = NINT - 1
270                     END IF
271                  END IF
272               END IF
273               IF( IWORK( I ).EQ.0 )
274     $            I2 = K
275               IF( CNT.EQ.I-1 ) THEN
276                  WERR( I ) = MID
277               ELSE IF( CNT.EQ.NRIGHT ) THEN
278                  W( I ) = MID
279               ELSE
280                  IWORK( N+I ) = CNT
281                  NINT = NINT + 1
282                  WERR( CNT+1 ) = MID
283                  W( CNT+1 ) = W( I )
284                  W( I ) = MID
285                  I = CNT + 1
286                  IWORK( N+I ) = NRIGHT
287               END IF
288            END IF
289            I = NRIGHT + 1
290  100    CONTINUE
291         NINT = NINT - OLNINT + I2
292         GO TO 80
293      END IF
294      DO 110 I = INITI1, INITI2
295         OPS = OPS + DBLE( 3 )
296         W( I ) = HALF*( WERR( I )+W( I ) )
297         WERR( I ) = W( I ) - WERR( I )
298  110 CONTINUE
299*
300      RETURN
301*
302*     End of DLARRB
303*
304      END
305