LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slaebz()

subroutine slaebz ( integer ijob,
integer nitmax,
integer n,
integer mmax,
integer minp,
integer nbmin,
real abstol,
real reltol,
real pivmin,
real, dimension( * ) d,
real, dimension( * ) e,
real, dimension( * ) e2,
integer, dimension( * ) nval,
real, dimension( mmax, * ) ab,
real, dimension( * ) c,
integer mout,
integer, dimension( mmax, * ) nab,
real, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz.

Download SLAEBZ + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SLAEBZ contains the iteration loops which compute and use the
!> function N(w), which is the count of eigenvalues of a symmetric
!> tridiagonal matrix T less than or equal to its argument  w.  It
!> performs a choice of two types of loops:
!>
!> IJOB=1, followed by
!> IJOB=2: It takes as input a list of intervals and returns a list of
!>         sufficiently small intervals whose union contains the same
!>         eigenvalues as the union of the original intervals.
!>         The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
!>         The output interval (AB(j,1),AB(j,2)] will contain
!>         eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
!>
!> IJOB=3: It performs a binary search in each input interval
!>         (AB(j,1),AB(j,2)] for a point  w(j)  such that
!>         N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
!>         the search.  If such a w(j) is found, then on output
!>         AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
!>         (AB(j,1),AB(j,2)] will be a small interval containing the
!>         point where N(w) jumps through NVAL(j), unless that point
!>         lies outside the initial interval.
!>
!> Note that the intervals are in all cases half-open intervals,
!> i.e., of the form  (a,b] , which includes  b  but not  a .
!>
!> To avoid underflow, the matrix should be scaled so that its largest
!> element is no greater than  overflow**(1/2) * underflow**(1/4)
!> in absolute value.  To assure the most accurate computation
!> of small eigenvalues, the matrix should be scaled to be
!> not much smaller than that, either.
!>
!> See W. Kahan , Report CS41, Computer Science Dept., Stanford
!> University, July 21, 1966
!>
!> Note: the arguments are, in general, *not* checked for unreasonable
!> values.
!> 
Parameters
[in]IJOB
!>          IJOB is INTEGER
!>          Specifies what is to be done:
!>          = 1:  Compute NAB for the initial intervals.
!>          = 2:  Perform bisection iteration to find eigenvalues of T.
!>          = 3:  Perform bisection iteration to invert N(w), i.e.,
!>                to find a point which has a specified number of
!>                eigenvalues of T to its left.
!>          Other values will cause SLAEBZ to return with INFO=-1.
!> 
[in]NITMAX
!>          NITMAX is INTEGER
!>          The maximum number of  of bisection to be
!>          performed, i.e., an interval of width W will not be made
!>          smaller than 2^(-NITMAX) * W.  If not all intervals
!>          have converged after NITMAX iterations, then INFO is set
!>          to the number of non-converged intervals.
!> 
[in]N
!>          N is INTEGER
!>          The dimension n of the tridiagonal matrix T.  It must be at
!>          least 1.
!> 
[in]MMAX
!>          MMAX is INTEGER
!>          The maximum number of intervals.  If more than MMAX intervals
!>          are generated, then SLAEBZ will quit with INFO=MMAX+1.
!> 
[in]MINP
!>          MINP is INTEGER
!>          The initial number of intervals.  It may not be greater than
!>          MMAX.
!> 
[in]NBMIN
!>          NBMIN is INTEGER
!>          The smallest number of intervals that should be processed
!>          using a vector loop.  If zero, then only the scalar loop
!>          will be used.
!> 
[in]ABSTOL
!>          ABSTOL is REAL
!>          The minimum (absolute) width of an interval.  When an
!>          interval is narrower than ABSTOL, or than RELTOL times the
!>          larger (in magnitude) endpoint, then it is considered to be
!>          sufficiently small, i.e., converged.  This must be at least
!>          zero.
!> 
[in]RELTOL
!>          RELTOL is REAL
!>          The minimum relative width of an interval.  When an interval
!>          is narrower than ABSTOL, or than RELTOL times the larger (in
!>          magnitude) endpoint, then it is considered to be
!>          sufficiently small, i.e., converged.  Note: this should
!>          always be at least radix*machine epsilon.
!> 
[in]PIVMIN
!>          PIVMIN is REAL
!>          The minimum absolute value of a  in the Sturm
!>          sequence loop.
!>          This must be at least  max |e(j)**2|*safe_min  and at
!>          least safe_min, where safe_min is at least
!>          the smallest number that can divide one without overflow.
!> 
[in]D
!>          D is REAL array, dimension (N)
!>          The diagonal elements of the tridiagonal matrix T.
!> 
[in]E
!>          E is REAL array, dimension (N)
!>          The offdiagonal elements of the tridiagonal matrix T in
!>          positions 1 through N-1.  E(N) is arbitrary.
!> 
[in]E2
!>          E2 is REAL array, dimension (N)
!>          The squares of the offdiagonal elements of the tridiagonal
!>          matrix T.  E2(N) is ignored.
!> 
[in,out]NVAL
!>          NVAL is INTEGER array, dimension (MINP)
!>          If IJOB=1 or 2, not referenced.
!>          If IJOB=3, the desired values of N(w).  The elements of NVAL
!>          will be reordered to correspond with the intervals in AB.
!>          Thus, NVAL(j) on output will not, in general be the same as
!>          NVAL(j) on input, but it will correspond with the interval
!>          (AB(j,1),AB(j,2)] on output.
!> 
[in,out]AB
!>          AB is REAL array, dimension (MMAX,2)
!>          The endpoints of the intervals.  AB(j,1) is  a(j), the left
!>          endpoint of the j-th interval, and AB(j,2) is b(j), the
!>          right endpoint of the j-th interval.  The input intervals
!>          will, in general, be modified, split, and reordered by the
!>          calculation.
!> 
[in,out]C
!>          C is REAL array, dimension (MMAX)
!>          If IJOB=1, ignored.
!>          If IJOB=2, workspace.
!>          If IJOB=3, then on input C(j) should be initialized to the
!>          first search point in the binary search.
!> 
[out]MOUT
!>          MOUT is INTEGER
!>          If IJOB=1, the number of eigenvalues in the intervals.
!>          If IJOB=2 or 3, the number of intervals output.
!>          If IJOB=3, MOUT will equal MINP.
!> 
[in,out]NAB
!>          NAB is INTEGER array, dimension (MMAX,2)
!>          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
!>          If IJOB=2, then on input, NAB(i,j) should be set.  It must
!>             satisfy the condition:
!>             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
!>             which means that in interval i only eigenvalues
!>             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
!>             NAB(i,j)=N(AB(i,j)), from a previous call to SLAEBZ with
!>             IJOB=1.
!>             On output, NAB(i,j) will contain
!>             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
!>             the input interval that the output interval
!>             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
!>             the input values of NAB(k,1) and NAB(k,2).
!>          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
!>             unless N(w) > NVAL(i) for all search points  w , in which
!>             case NAB(i,1) will not be modified, i.e., the output
!>             value will be the same as the input value (modulo
!>             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
!>             for all search points  w , in which case NAB(i,2) will
!>             not be modified.  Normally, NAB should be set to some
!>             distinctive value(s) before SLAEBZ is called.
!> 
[out]WORK
!>          WORK is REAL array, dimension (MMAX)
!>          Workspace.
!> 
[out]IWORK
!>          IWORK is INTEGER array, dimension (MMAX)
!>          Workspace.
!> 
[out]INFO
!>          INFO is INTEGER
!>          = 0:       All intervals converged.
!>          = 1--MMAX: The last INFO intervals did not converge.
!>          = MMAX+1:  More than MMAX intervals were generated.
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>      This routine is intended to be called only by other LAPACK
!>  routines, thus the interface is less user-friendly.  It is intended
!>  for two purposes:
!>
!>  (a) finding eigenvalues.  In this case, SLAEBZ should have one or
!>      more initial intervals set up in AB, and SLAEBZ should be called
!>      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
!>      Intervals with no eigenvalues would usually be thrown out at
!>      this point.  Also, if not all the eigenvalues in an interval i
!>      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
!>      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
!>      eigenvalue.  SLAEBZ is then called with IJOB=2 and MMAX
!>      no smaller than the value of MOUT returned by the call with
!>      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
!>      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
!>      tolerance specified by ABSTOL and RELTOL.
!>
!>  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
!>      In this case, start with a Gershgorin interval  (a,b).  Set up
!>      AB to contain 2 search intervals, both initially (a,b).  One
!>      NVAL element should contain  f-1  and the other should contain  l
!>      , while C should contain a and b, resp.  NAB(i,1) should be -1
!>      and NAB(i,2) should be N+1, to flag an error if the desired
!>      interval does not lie in (a,b).  SLAEBZ is then called with
!>      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
!>      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
!>      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
!>      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
!>      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
!>      w(l-r)=...=w(l+k) are handled similarly.
!> 

Definition at line 314 of file slaebz.f.

317*
318* -- LAPACK auxiliary routine --
319* -- LAPACK is a software package provided by Univ. of Tennessee, --
320* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
321*
322* .. Scalar Arguments ..
323 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
324 REAL ABSTOL, PIVMIN, RELTOL
325* ..
326* .. Array Arguments ..
327 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
328 REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
329 $ WORK( * )
330* ..
331*
332* =====================================================================
333*
334* .. Parameters ..
335 REAL ZERO, TWO, HALF
336 parameter( zero = 0.0e0, two = 2.0e0,
337 $ half = 1.0e0 / two )
338* ..
339* .. Local Scalars ..
340 INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
341 $ KLNEW
342 REAL TMP1, TMP2
343* ..
344* .. Intrinsic Functions ..
345 INTRINSIC abs, max, min
346* ..
347* .. Executable Statements ..
348*
349* Check for Errors
350*
351 info = 0
352 IF( ijob.LT.1 .OR. ijob.GT.3 ) THEN
353 info = -1
354 RETURN
355 END IF
356*
357* Initialize NAB
358*
359 IF( ijob.EQ.1 ) THEN
360*
361* Compute the number of eigenvalues in the initial intervals.
362*
363 mout = 0
364 DO 30 ji = 1, minp
365 DO 20 jp = 1, 2
366 tmp1 = d( 1 ) - ab( ji, jp )
367 IF( abs( tmp1 ).LT.pivmin )
368 $ tmp1 = -pivmin
369 nab( ji, jp ) = 0
370 IF( tmp1.LE.zero )
371 $ nab( ji, jp ) = 1
372*
373 DO 10 j = 2, n
374 tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp )
375 IF( abs( tmp1 ).LT.pivmin )
376 $ tmp1 = -pivmin
377 IF( tmp1.LE.zero )
378 $ nab( ji, jp ) = nab( ji, jp ) + 1
379 10 CONTINUE
380 20 CONTINUE
381 mout = mout + nab( ji, 2 ) - nab( ji, 1 )
382 30 CONTINUE
383 RETURN
384 END IF
385*
386* Initialize for loop
387*
388* KF and KL have the following meaning:
389* Intervals 1,...,KF-1 have converged.
390* Intervals KF,...,KL still need to be refined.
391*
392 kf = 1
393 kl = minp
394*
395* If IJOB=2, initialize C.
396* If IJOB=3, use the user-supplied starting point.
397*
398 IF( ijob.EQ.2 ) THEN
399 DO 40 ji = 1, minp
400 c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
401 40 CONTINUE
402 END IF
403*
404* Iteration loop
405*
406 DO 130 jit = 1, nitmax
407*
408* Loop over intervals
409*
410 IF( kl-kf+1.GE.nbmin .AND. nbmin.GT.0 ) THEN
411*
412* Begin of Parallel Version of the loop
413*
414 DO 60 ji = kf, kl
415*
416* Compute N(c), the number of eigenvalues less than c
417*
418 work( ji ) = d( 1 ) - c( ji )
419 iwork( ji ) = 0
420 IF( work( ji ).LE.pivmin ) THEN
421 iwork( ji ) = 1
422 work( ji ) = min( work( ji ), -pivmin )
423 END IF
424*
425 DO 50 j = 2, n
426 work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji )
427 IF( work( ji ).LE.pivmin ) THEN
428 iwork( ji ) = iwork( ji ) + 1
429 work( ji ) = min( work( ji ), -pivmin )
430 END IF
431 50 CONTINUE
432 60 CONTINUE
433*
434 IF( ijob.LE.2 ) THEN
435*
436* IJOB=2: Choose all intervals containing eigenvalues.
437*
438 klnew = kl
439 DO 70 ji = kf, kl
440*
441* Insure that N(w) is monotone
442*
443 iwork( ji ) = min( nab( ji, 2 ),
444 $ max( nab( ji, 1 ), iwork( ji ) ) )
445*
446* Update the Queue -- add intervals if both halves
447* contain eigenvalues.
448*
449 IF( iwork( ji ).EQ.nab( ji, 2 ) ) THEN
450*
451* No eigenvalue in the upper interval:
452* just use the lower interval.
453*
454 ab( ji, 2 ) = c( ji )
455*
456 ELSE IF( iwork( ji ).EQ.nab( ji, 1 ) ) THEN
457*
458* No eigenvalue in the lower interval:
459* just use the upper interval.
460*
461 ab( ji, 1 ) = c( ji )
462 ELSE
463 klnew = klnew + 1
464 IF( klnew.LE.mmax ) THEN
465*
466* Eigenvalue in both intervals -- add upper to
467* queue.
468*
469 ab( klnew, 2 ) = ab( ji, 2 )
470 nab( klnew, 2 ) = nab( ji, 2 )
471 ab( klnew, 1 ) = c( ji )
472 nab( klnew, 1 ) = iwork( ji )
473 ab( ji, 2 ) = c( ji )
474 nab( ji, 2 ) = iwork( ji )
475 ELSE
476 info = mmax + 1
477 END IF
478 END IF
479 70 CONTINUE
480 IF( info.NE.0 )
481 $ RETURN
482 kl = klnew
483 ELSE
484*
485* IJOB=3: Binary search. Keep only the interval containing
486* w s.t. N(w) = NVAL
487*
488 DO 80 ji = kf, kl
489 IF( iwork( ji ).LE.nval( ji ) ) THEN
490 ab( ji, 1 ) = c( ji )
491 nab( ji, 1 ) = iwork( ji )
492 END IF
493 IF( iwork( ji ).GE.nval( ji ) ) THEN
494 ab( ji, 2 ) = c( ji )
495 nab( ji, 2 ) = iwork( ji )
496 END IF
497 80 CONTINUE
498 END IF
499*
500 ELSE
501*
502* End of Parallel Version of the loop
503*
504* Begin of Serial Version of the loop
505*
506 klnew = kl
507 DO 100 ji = kf, kl
508*
509* Compute N(w), the number of eigenvalues less than w
510*
511 tmp1 = c( ji )
512 tmp2 = d( 1 ) - tmp1
513 itmp1 = 0
514 IF( tmp2.LE.pivmin ) THEN
515 itmp1 = 1
516 tmp2 = min( tmp2, -pivmin )
517 END IF
518*
519 DO 90 j = 2, n
520 tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1
521 IF( tmp2.LE.pivmin ) THEN
522 itmp1 = itmp1 + 1
523 tmp2 = min( tmp2, -pivmin )
524 END IF
525 90 CONTINUE
526*
527 IF( ijob.LE.2 ) THEN
528*
529* IJOB=2: Choose all intervals containing eigenvalues.
530*
531* Insure that N(w) is monotone
532*
533 itmp1 = min( nab( ji, 2 ),
534 $ max( nab( ji, 1 ), itmp1 ) )
535*
536* Update the Queue -- add intervals if both halves
537* contain eigenvalues.
538*
539 IF( itmp1.EQ.nab( ji, 2 ) ) THEN
540*
541* No eigenvalue in the upper interval:
542* just use the lower interval.
543*
544 ab( ji, 2 ) = tmp1
545*
546 ELSE IF( itmp1.EQ.nab( ji, 1 ) ) THEN
547*
548* No eigenvalue in the lower interval:
549* just use the upper interval.
550*
551 ab( ji, 1 ) = tmp1
552 ELSE IF( klnew.LT.mmax ) THEN
553*
554* Eigenvalue in both intervals -- add upper to queue.
555*
556 klnew = klnew + 1
557 ab( klnew, 2 ) = ab( ji, 2 )
558 nab( klnew, 2 ) = nab( ji, 2 )
559 ab( klnew, 1 ) = tmp1
560 nab( klnew, 1 ) = itmp1
561 ab( ji, 2 ) = tmp1
562 nab( ji, 2 ) = itmp1
563 ELSE
564 info = mmax + 1
565 RETURN
566 END IF
567 ELSE
568*
569* IJOB=3: Binary search. Keep only the interval
570* containing w s.t. N(w) = NVAL
571*
572 IF( itmp1.LE.nval( ji ) ) THEN
573 ab( ji, 1 ) = tmp1
574 nab( ji, 1 ) = itmp1
575 END IF
576 IF( itmp1.GE.nval( ji ) ) THEN
577 ab( ji, 2 ) = tmp1
578 nab( ji, 2 ) = itmp1
579 END IF
580 END IF
581 100 CONTINUE
582 kl = klnew
583*
584 END IF
585*
586* Check for convergence
587*
588 kfnew = kf
589 DO 110 ji = kf, kl
590 tmp1 = abs( ab( ji, 2 )-ab( ji, 1 ) )
591 tmp2 = max( abs( ab( ji, 2 ) ), abs( ab( ji, 1 ) ) )
592 IF( tmp1.LT.max( abstol, pivmin, reltol*tmp2 ) .OR.
593 $ nab( ji, 1 ).GE.nab( ji, 2 ) ) THEN
594*
595* Converged -- Swap with position KFNEW,
596* then increment KFNEW
597*
598 IF( ji.GT.kfnew ) THEN
599 tmp1 = ab( ji, 1 )
600 tmp2 = ab( ji, 2 )
601 itmp1 = nab( ji, 1 )
602 itmp2 = nab( ji, 2 )
603 ab( ji, 1 ) = ab( kfnew, 1 )
604 ab( ji, 2 ) = ab( kfnew, 2 )
605 nab( ji, 1 ) = nab( kfnew, 1 )
606 nab( ji, 2 ) = nab( kfnew, 2 )
607 ab( kfnew, 1 ) = tmp1
608 ab( kfnew, 2 ) = tmp2
609 nab( kfnew, 1 ) = itmp1
610 nab( kfnew, 2 ) = itmp2
611 IF( ijob.EQ.3 ) THEN
612 itmp1 = nval( ji )
613 nval( ji ) = nval( kfnew )
614 nval( kfnew ) = itmp1
615 END IF
616 END IF
617 kfnew = kfnew + 1
618 END IF
619 110 CONTINUE
620 kf = kfnew
621*
622* Choose Midpoints
623*
624 DO 120 ji = kf, kl
625 c( ji ) = half*( ab( ji, 1 )+ab( ji, 2 ) )
626 120 CONTINUE
627*
628* If no more intervals to refine, quit.
629*
630 IF( kf.GT.kl )
631 $ GO TO 140
632 130 CONTINUE
633*
634* Converged
635*
636 140 CONTINUE
637 info = max( kl+1-kf, 0 )
638 mout = kl
639*
640 RETURN
641*
642* End of SLAEBZ
643*
Here is the caller graph for this function: