001:       SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
002:      $                   RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
003:      $                   NAB, WORK, IWORK, INFO )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
011:       DOUBLE PRECISION   ABSTOL, PIVMIN, RELTOL
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
015:       DOUBLE PRECISION   AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
016:      $                   WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DLAEBZ contains the iteration loops which compute and use the
023: *  function N(w), which is the count of eigenvalues of a symmetric
024: *  tridiagonal matrix T less than or equal to its argument  w.  It
025: *  performs a choice of two types of loops:
026: *
027: *  IJOB=1, followed by
028: *  IJOB=2: It takes as input a list of intervals and returns a list of
029: *          sufficiently small intervals whose union contains the same
030: *          eigenvalues as the union of the original intervals.
031: *          The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
032: *          The output interval (AB(j,1),AB(j,2)] will contain
033: *          eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
034: *
035: *  IJOB=3: It performs a binary search in each input interval
036: *          (AB(j,1),AB(j,2)] for a point  w(j)  such that
037: *          N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
038: *          the search.  If such a w(j) is found, then on output
039: *          AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
040: *          (AB(j,1),AB(j,2)] will be a small interval containing the
041: *          point where N(w) jumps through NVAL(j), unless that point
042: *          lies outside the initial interval.
043: *
044: *  Note that the intervals are in all cases half-open intervals,
045: *  i.e., of the form  (a,b] , which includes  b  but not  a .
046: *
047: *  To avoid underflow, the matrix should be scaled so that its largest
048: *  element is no greater than  overflow**(1/2) * underflow**(1/4)
049: *  in absolute value.  To assure the most accurate computation
050: *  of small eigenvalues, the matrix should be scaled to be
051: *  not much smaller than that, either.
052: *
053: *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
054: *  Matrix", Report CS41, Computer Science Dept., Stanford
055: *  University, July 21, 1966
056: *
057: *  Note: the arguments are, in general, *not* checked for unreasonable
058: *  values.
059: *
060: *  Arguments
061: *  =========
062: *
063: *  IJOB    (input) INTEGER
064: *          Specifies what is to be done:
065: *          = 1:  Compute NAB for the initial intervals.
066: *          = 2:  Perform bisection iteration to find eigenvalues of T.
067: *          = 3:  Perform bisection iteration to invert N(w), i.e.,
068: *                to find a point which has a specified number of
069: *                eigenvalues of T to its left.
070: *          Other values will cause DLAEBZ to return with INFO=-1.
071: *
072: *  NITMAX  (input) INTEGER
073: *          The maximum number of "levels" of bisection to be
074: *          performed, i.e., an interval of width W will not be made
075: *          smaller than 2^(-NITMAX) * W.  If not all intervals
076: *          have converged after NITMAX iterations, then INFO is set
077: *          to the number of non-converged intervals.
078: *
079: *  N       (input) INTEGER
080: *          The dimension n of the tridiagonal matrix T.  It must be at
081: *          least 1.
082: *
083: *  MMAX    (input) INTEGER
084: *          The maximum number of intervals.  If more than MMAX intervals
085: *          are generated, then DLAEBZ will quit with INFO=MMAX+1.
086: *
087: *  MINP    (input) INTEGER
088: *          The initial number of intervals.  It may not be greater than
089: *          MMAX.
090: *
091: *  NBMIN   (input) INTEGER
092: *          The smallest number of intervals that should be processed
093: *          using a vector loop.  If zero, then only the scalar loop
094: *          will be used.
095: *
096: *  ABSTOL  (input) DOUBLE PRECISION
097: *          The minimum (absolute) width of an interval.  When an
098: *          interval is narrower than ABSTOL, or than RELTOL times the
099: *          larger (in magnitude) endpoint, then it is considered to be
100: *          sufficiently small, i.e., converged.  This must be at least
101: *          zero.
102: *
103: *  RELTOL  (input) DOUBLE PRECISION
104: *          The minimum relative width of an interval.  When an interval
105: *          is narrower than ABSTOL, or than RELTOL times the larger (in
106: *          magnitude) endpoint, then it is considered to be
107: *          sufficiently small, i.e., converged.  Note: this should
108: *          always be at least radix*machine epsilon.
109: *
110: *  PIVMIN  (input) DOUBLE PRECISION
111: *          The minimum absolute value of a "pivot" in the Sturm
112: *          sequence loop.  This *must* be at least  max |e(j)**2| *
113: *          safe_min  and at least safe_min, where safe_min is at least
114: *          the smallest number that can divide one without overflow.
115: *
116: *  D       (input) DOUBLE PRECISION array, dimension (N)
117: *          The diagonal elements of the tridiagonal matrix T.
118: *
119: *  E       (input) DOUBLE PRECISION array, dimension (N)
120: *          The offdiagonal elements of the tridiagonal matrix T in
121: *          positions 1 through N-1.  E(N) is arbitrary.
122: *
123: *  E2      (input) DOUBLE PRECISION array, dimension (N)
124: *          The squares of the offdiagonal elements of the tridiagonal
125: *          matrix T.  E2(N) is ignored.
126: *
127: *  NVAL    (input/output) INTEGER array, dimension (MINP)
128: *          If IJOB=1 or 2, not referenced.
129: *          If IJOB=3, the desired values of N(w).  The elements of NVAL
130: *          will be reordered to correspond with the intervals in AB.
131: *          Thus, NVAL(j) on output will not, in general be the same as
132: *          NVAL(j) on input, but it will correspond with the interval
133: *          (AB(j,1),AB(j,2)] on output.
134: *
135: *  AB      (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
136: *          The endpoints of the intervals.  AB(j,1) is  a(j), the left
137: *          endpoint of the j-th interval, and AB(j,2) is b(j), the
138: *          right endpoint of the j-th interval.  The input intervals
139: *          will, in general, be modified, split, and reordered by the
140: *          calculation.
141: *
142: *  C       (input/output) DOUBLE PRECISION array, dimension (MMAX)
143: *          If IJOB=1, ignored.
144: *          If IJOB=2, workspace.
145: *          If IJOB=3, then on input C(j) should be initialized to the
146: *          first search point in the binary search.
147: *
148: *  MOUT    (output) INTEGER
149: *          If IJOB=1, the number of eigenvalues in the intervals.
150: *          If IJOB=2 or 3, the number of intervals output.
151: *          If IJOB=3, MOUT will equal MINP.
152: *
153: *  NAB     (input/output) INTEGER array, dimension (MMAX,2)
154: *          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
155: *          If IJOB=2, then on input, NAB(i,j) should be set.  It must
156: *             satisfy the condition:
157: *             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
158: *             which means that in interval i only eigenvalues
159: *             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
160: *             NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
161: *             IJOB=1.
162: *             On output, NAB(i,j) will contain
163: *             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
164: *             the input interval that the output interval
165: *             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
166: *             the input values of NAB(k,1) and NAB(k,2).
167: *          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
168: *             unless N(w) > NVAL(i) for all search points  w , in which
169: *             case NAB(i,1) will not be modified, i.e., the output
170: *             value will be the same as the input value (modulo
171: *             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
172: *             for all search points  w , in which case NAB(i,2) will
173: *             not be modified.  Normally, NAB should be set to some
174: *             distinctive value(s) before DLAEBZ is called.
175: *
176: *  WORK    (workspace) DOUBLE PRECISION array, dimension (MMAX)
177: *          Workspace.
178: *
179: *  IWORK   (workspace) INTEGER array, dimension (MMAX)
180: *          Workspace.
181: *
182: *  INFO    (output) INTEGER
183: *          = 0:       All intervals converged.
184: *          = 1--MMAX: The last INFO intervals did not converge.
185: *          = MMAX+1:  More than MMAX intervals were generated.
186: *
187: *  Further Details
188: *  ===============
189: *
190: *      This routine is intended to be called only by other LAPACK
191: *  routines, thus the interface is less user-friendly.  It is intended
192: *  for two purposes:
193: *
194: *  (a) finding eigenvalues.  In this case, DLAEBZ should have one or
195: *      more initial intervals set up in AB, and DLAEBZ should be called
196: *      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
197: *      Intervals with no eigenvalues would usually be thrown out at
198: *      this point.  Also, if not all the eigenvalues in an interval i
199: *      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
200: *      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
201: *      eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX
202: *      no smaller than the value of MOUT returned by the call with
203: *      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
204: *      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
205: *      tolerance specified by ABSTOL and RELTOL.
206: *
207: *  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
208: *      In this case, start with a Gershgorin interval  (a,b).  Set up
209: *      AB to contain 2 search intervals, both initially (a,b).  One
210: *      NVAL element should contain  f-1  and the other should contain  l
211: *      , while C should contain a and b, resp.  NAB(i,1) should be -1
212: *      and NAB(i,2) should be N+1, to flag an error if the desired
213: *      interval does not lie in (a,b).  DLAEBZ is then called with
214: *      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
215: *      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
216: *      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
217: *      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
218: *      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
219: *      w(l-r)=...=w(l+k) are handled similarly.
220: *
221: *  =====================================================================
222: *
223: *     .. Parameters ..
224:       DOUBLE PRECISION   ZERO, TWO, HALF
225:       PARAMETER          ( ZERO = 0.0D0, TWO = 2.0D0,
226:      $                   HALF = 1.0D0 / TWO )
227: *     ..
228: *     .. Local Scalars ..
229:       INTEGER            ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
230:      $                   KLNEW
231:       DOUBLE PRECISION   TMP1, TMP2
232: *     ..
233: *     .. Intrinsic Functions ..
234:       INTRINSIC          ABS, MAX, MIN
235: *     ..
236: *     .. Executable Statements ..
237: *
238: *     Check for Errors
239: *
240:       INFO = 0
241:       IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
242:          INFO = -1
243:          RETURN
244:       END IF
245: *
246: *     Initialize NAB
247: *
248:       IF( IJOB.EQ.1 ) THEN
249: *
250: *        Compute the number of eigenvalues in the initial intervals.
251: *
252:          MOUT = 0
253: *DIR$ NOVECTOR
254:          DO 30 JI = 1, MINP
255:             DO 20 JP = 1, 2
256:                TMP1 = D( 1 ) - AB( JI, JP )
257:                IF( ABS( TMP1 ).LT.PIVMIN )
258:      $            TMP1 = -PIVMIN
259:                NAB( JI, JP ) = 0
260:                IF( TMP1.LE.ZERO )
261:      $            NAB( JI, JP ) = 1
262: *
263:                DO 10 J = 2, N
264:                   TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
265:                   IF( ABS( TMP1 ).LT.PIVMIN )
266:      $               TMP1 = -PIVMIN
267:                   IF( TMP1.LE.ZERO )
268:      $               NAB( JI, JP ) = NAB( JI, JP ) + 1
269:    10          CONTINUE
270:    20       CONTINUE
271:             MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
272:    30    CONTINUE
273:          RETURN
274:       END IF
275: *
276: *     Initialize for loop
277: *
278: *     KF and KL have the following meaning:
279: *        Intervals 1,...,KF-1 have converged.
280: *        Intervals KF,...,KL  still need to be refined.
281: *
282:       KF = 1
283:       KL = MINP
284: *
285: *     If IJOB=2, initialize C.
286: *     If IJOB=3, use the user-supplied starting point.
287: *
288:       IF( IJOB.EQ.2 ) THEN
289:          DO 40 JI = 1, MINP
290:             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
291:    40    CONTINUE
292:       END IF
293: *
294: *     Iteration loop
295: *
296:       DO 130 JIT = 1, NITMAX
297: *
298: *        Loop over intervals
299: *
300:          IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
301: *
302: *           Begin of Parallel Version of the loop
303: *
304:             DO 60 JI = KF, KL
305: *
306: *              Compute N(c), the number of eigenvalues less than c
307: *
308:                WORK( JI ) = D( 1 ) - C( JI )
309:                IWORK( JI ) = 0
310:                IF( WORK( JI ).LE.PIVMIN ) THEN
311:                   IWORK( JI ) = 1
312:                   WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
313:                END IF
314: *
315:                DO 50 J = 2, N
316:                   WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
317:                   IF( WORK( JI ).LE.PIVMIN ) THEN
318:                      IWORK( JI ) = IWORK( JI ) + 1
319:                      WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
320:                   END IF
321:    50          CONTINUE
322:    60       CONTINUE
323: *
324:             IF( IJOB.LE.2 ) THEN
325: *
326: *              IJOB=2: Choose all intervals containing eigenvalues.
327: *
328:                KLNEW = KL
329:                DO 70 JI = KF, KL
330: *
331: *                 Insure that N(w) is monotone
332: *
333:                   IWORK( JI ) = MIN( NAB( JI, 2 ),
334:      $                          MAX( NAB( JI, 1 ), IWORK( JI ) ) )
335: *
336: *                 Update the Queue -- add intervals if both halves
337: *                 contain eigenvalues.
338: *
339:                   IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
340: *
341: *                    No eigenvalue in the upper interval:
342: *                    just use the lower interval.
343: *
344:                      AB( JI, 2 ) = C( JI )
345: *
346:                   ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
347: *
348: *                    No eigenvalue in the lower interval:
349: *                    just use the upper interval.
350: *
351:                      AB( JI, 1 ) = C( JI )
352:                   ELSE
353:                      KLNEW = KLNEW + 1
354:                      IF( KLNEW.LE.MMAX ) THEN
355: *
356: *                       Eigenvalue in both intervals -- add upper to
357: *                       queue.
358: *
359:                         AB( KLNEW, 2 ) = AB( JI, 2 )
360:                         NAB( KLNEW, 2 ) = NAB( JI, 2 )
361:                         AB( KLNEW, 1 ) = C( JI )
362:                         NAB( KLNEW, 1 ) = IWORK( JI )
363:                         AB( JI, 2 ) = C( JI )
364:                         NAB( JI, 2 ) = IWORK( JI )
365:                      ELSE
366:                         INFO = MMAX + 1
367:                      END IF
368:                   END IF
369:    70          CONTINUE
370:                IF( INFO.NE.0 )
371:      $            RETURN
372:                KL = KLNEW
373:             ELSE
374: *
375: *              IJOB=3: Binary search.  Keep only the interval containing
376: *                      w   s.t. N(w) = NVAL
377: *
378:                DO 80 JI = KF, KL
379:                   IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
380:                      AB( JI, 1 ) = C( JI )
381:                      NAB( JI, 1 ) = IWORK( JI )
382:                   END IF
383:                   IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
384:                      AB( JI, 2 ) = C( JI )
385:                      NAB( JI, 2 ) = IWORK( JI )
386:                   END IF
387:    80          CONTINUE
388:             END IF
389: *
390:          ELSE
391: *
392: *           End of Parallel Version of the loop
393: *
394: *           Begin of Serial Version of the loop
395: *
396:             KLNEW = KL
397:             DO 100 JI = KF, KL
398: *
399: *              Compute N(w), the number of eigenvalues less than w
400: *
401:                TMP1 = C( JI )
402:                TMP2 = D( 1 ) - TMP1
403:                ITMP1 = 0
404:                IF( TMP2.LE.PIVMIN ) THEN
405:                   ITMP1 = 1
406:                   TMP2 = MIN( TMP2, -PIVMIN )
407:                END IF
408: *
409: *              A series of compiler directives to defeat vectorization
410: *              for the next loop
411: *
412: *$PL$ CMCHAR=' '
413: CDIR$          NEXTSCALAR
414: C$DIR          SCALAR
415: CDIR$          NEXT SCALAR
416: CVD$L          NOVECTOR
417: CDEC$          NOVECTOR
418: CVD$           NOVECTOR
419: *VDIR          NOVECTOR
420: *VOCL          LOOP,SCALAR
421: CIBM           PREFER SCALAR
422: *$PL$ CMCHAR='*'
423: *
424:                DO 90 J = 2, N
425:                   TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
426:                   IF( TMP2.LE.PIVMIN ) THEN
427:                      ITMP1 = ITMP1 + 1
428:                      TMP2 = MIN( TMP2, -PIVMIN )
429:                   END IF
430:    90          CONTINUE
431: *
432:                IF( IJOB.LE.2 ) THEN
433: *
434: *                 IJOB=2: Choose all intervals containing eigenvalues.
435: *
436: *                 Insure that N(w) is monotone
437: *
438:                   ITMP1 = MIN( NAB( JI, 2 ),
439:      $                    MAX( NAB( JI, 1 ), ITMP1 ) )
440: *
441: *                 Update the Queue -- add intervals if both halves
442: *                 contain eigenvalues.
443: *
444:                   IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
445: *
446: *                    No eigenvalue in the upper interval:
447: *                    just use the lower interval.
448: *
449:                      AB( JI, 2 ) = TMP1
450: *
451:                   ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
452: *
453: *                    No eigenvalue in the lower interval:
454: *                    just use the upper interval.
455: *
456:                      AB( JI, 1 ) = TMP1
457:                   ELSE IF( KLNEW.LT.MMAX ) THEN
458: *
459: *                    Eigenvalue in both intervals -- add upper to queue.
460: *
461:                      KLNEW = KLNEW + 1
462:                      AB( KLNEW, 2 ) = AB( JI, 2 )
463:                      NAB( KLNEW, 2 ) = NAB( JI, 2 )
464:                      AB( KLNEW, 1 ) = TMP1
465:                      NAB( KLNEW, 1 ) = ITMP1
466:                      AB( JI, 2 ) = TMP1
467:                      NAB( JI, 2 ) = ITMP1
468:                   ELSE
469:                      INFO = MMAX + 1
470:                      RETURN
471:                   END IF
472:                ELSE
473: *
474: *                 IJOB=3: Binary search.  Keep only the interval
475: *                         containing  w  s.t. N(w) = NVAL
476: *
477:                   IF( ITMP1.LE.NVAL( JI ) ) THEN
478:                      AB( JI, 1 ) = TMP1
479:                      NAB( JI, 1 ) = ITMP1
480:                   END IF
481:                   IF( ITMP1.GE.NVAL( JI ) ) THEN
482:                      AB( JI, 2 ) = TMP1
483:                      NAB( JI, 2 ) = ITMP1
484:                   END IF
485:                END IF
486:   100       CONTINUE
487:             KL = KLNEW
488: *
489: *           End of Serial Version of the loop
490: *
491:          END IF
492: *
493: *        Check for convergence
494: *
495:          KFNEW = KF
496:          DO 110 JI = KF, KL
497:             TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
498:             TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
499:             IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
500:      $          NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
501: *
502: *              Converged -- Swap with position KFNEW,
503: *                           then increment KFNEW
504: *
505:                IF( JI.GT.KFNEW ) THEN
506:                   TMP1 = AB( JI, 1 )
507:                   TMP2 = AB( JI, 2 )
508:                   ITMP1 = NAB( JI, 1 )
509:                   ITMP2 = NAB( JI, 2 )
510:                   AB( JI, 1 ) = AB( KFNEW, 1 )
511:                   AB( JI, 2 ) = AB( KFNEW, 2 )
512:                   NAB( JI, 1 ) = NAB( KFNEW, 1 )
513:                   NAB( JI, 2 ) = NAB( KFNEW, 2 )
514:                   AB( KFNEW, 1 ) = TMP1
515:                   AB( KFNEW, 2 ) = TMP2
516:                   NAB( KFNEW, 1 ) = ITMP1
517:                   NAB( KFNEW, 2 ) = ITMP2
518:                   IF( IJOB.EQ.3 ) THEN
519:                      ITMP1 = NVAL( JI )
520:                      NVAL( JI ) = NVAL( KFNEW )
521:                      NVAL( KFNEW ) = ITMP1
522:                   END IF
523:                END IF
524:                KFNEW = KFNEW + 1
525:             END IF
526:   110    CONTINUE
527:          KF = KFNEW
528: *
529: *        Choose Midpoints
530: *
531:          DO 120 JI = KF, KL
532:             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
533:   120    CONTINUE
534: *
535: *        If no more intervals to refine, quit.
536: *
537:          IF( KF.GT.KL )
538:      $      GO TO 140
539:   130 CONTINUE
540: *
541: *     Converged
542: *
543:   140 CONTINUE
544:       INFO = MAX( KL+1-KF, 0 )
545:       MOUT = KL
546: *
547:       RETURN
548: *
549: *     End of DLAEBZ
550: *
551:       END
552: