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