diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/slaebz.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/slaebz.f')
-rw-r--r-- | SRC/slaebz.f | 551 |
1 files changed, 551 insertions, 0 deletions
diff --git a/SRC/slaebz.f b/SRC/slaebz.f new file mode 100644 index 00000000..82af82af --- /dev/null +++ b/SRC/slaebz.f @@ -0,0 +1,551 @@ + SUBROUTINE SLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, + $ RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, + $ NAB, WORK, IWORK, INFO ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX + REAL ABSTOL, PIVMIN, RELTOL +* .. +* .. Array Arguments .. + INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) + REAL AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), + $ WORK( * ) +* .. +* +* 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 "Accurate Eigenvalues of a Symmetric Tridiagonal +* Matrix", Report CS41, Computer Science Dept., Stanford +* University, July 21, 1966 +* +* Note: the arguments are, in general, *not* checked for unreasonable +* values. +* +* Arguments +* ========= +* +* IJOB (input) 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. +* +* NITMAX (input) INTEGER +* The maximum number of "levels" 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. +* +* N (input) INTEGER +* The dimension n of the tridiagonal matrix T. It must be at +* least 1. +* +* MMAX (input) INTEGER +* The maximum number of intervals. If more than MMAX intervals +* are generated, then SLAEBZ will quit with INFO=MMAX+1. +* +* MINP (input) INTEGER +* The initial number of intervals. It may not be greater than +* MMAX. +* +* NBMIN (input) INTEGER +* The smallest number of intervals that should be processed +* using a vector loop. If zero, then only the scalar loop +* will be used. +* +* ABSTOL (input) 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. +* +* RELTOL (input) 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. +* +* PIVMIN (input) REAL +* The minimum absolute value of a "pivot" 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. +* +* D (input) REAL array, dimension (N) +* The diagonal elements of the tridiagonal matrix T. +* +* E (input) REAL array, dimension (N) +* The offdiagonal elements of the tridiagonal matrix T in +* positions 1 through N-1. E(N) is arbitrary. +* +* E2 (input) REAL array, dimension (N) +* The squares of the offdiagonal elements of the tridiagonal +* matrix T. E2(N) is ignored. +* +* NVAL (input/output) 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. +* +* AB (input/output) 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. +* +* C (input/output) 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. +* +* MOUT (output) 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. +* +* NAB (input/output) 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. +* +* WORK (workspace) REAL array, dimension (MMAX) +* Workspace. +* +* IWORK (workspace) INTEGER array, dimension (MMAX) +* Workspace. +* +* INFO (output) INTEGER +* = 0: All intervals converged. +* = 1--MMAX: The last INFO intervals did not converge. +* = MMAX+1: More than MMAX intervals were generated. +* +* 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. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, TWO, HALF + PARAMETER ( ZERO = 0.0E0, TWO = 2.0E0, + $ HALF = 1.0E0 / TWO ) +* .. +* .. Local Scalars .. + INTEGER ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL, + $ KLNEW + REAL TMP1, TMP2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Executable Statements .. +* +* Check for Errors +* + INFO = 0 + IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN + INFO = -1 + RETURN + END IF +* +* Initialize NAB +* + IF( IJOB.EQ.1 ) THEN +* +* Compute the number of eigenvalues in the initial intervals. +* + MOUT = 0 +CDIR$ NOVECTOR + DO 30 JI = 1, MINP + DO 20 JP = 1, 2 + TMP1 = D( 1 ) - AB( JI, JP ) + IF( ABS( TMP1 ).LT.PIVMIN ) + $ TMP1 = -PIVMIN + NAB( JI, JP ) = 0 + IF( TMP1.LE.ZERO ) + $ NAB( JI, JP ) = 1 +* + DO 10 J = 2, N + TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP ) + IF( ABS( TMP1 ).LT.PIVMIN ) + $ TMP1 = -PIVMIN + IF( TMP1.LE.ZERO ) + $ NAB( JI, JP ) = NAB( JI, JP ) + 1 + 10 CONTINUE + 20 CONTINUE + MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 ) + 30 CONTINUE + RETURN + END IF +* +* Initialize for loop +* +* KF and KL have the following meaning: +* Intervals 1,...,KF-1 have converged. +* Intervals KF,...,KL still need to be refined. +* + KF = 1 + KL = MINP +* +* If IJOB=2, initialize C. +* If IJOB=3, use the user-supplied starting point. +* + IF( IJOB.EQ.2 ) THEN + DO 40 JI = 1, MINP + C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) + 40 CONTINUE + END IF +* +* Iteration loop +* + DO 130 JIT = 1, NITMAX +* +* Loop over intervals +* + IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN +* +* Begin of Parallel Version of the loop +* + DO 60 JI = KF, KL +* +* Compute N(c), the number of eigenvalues less than c +* + WORK( JI ) = D( 1 ) - C( JI ) + IWORK( JI ) = 0 + IF( WORK( JI ).LE.PIVMIN ) THEN + IWORK( JI ) = 1 + WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) + END IF +* + DO 50 J = 2, N + WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI ) + IF( WORK( JI ).LE.PIVMIN ) THEN + IWORK( JI ) = IWORK( JI ) + 1 + WORK( JI ) = MIN( WORK( JI ), -PIVMIN ) + END IF + 50 CONTINUE + 60 CONTINUE +* + IF( IJOB.LE.2 ) THEN +* +* IJOB=2: Choose all intervals containing eigenvalues. +* + KLNEW = KL + DO 70 JI = KF, KL +* +* Insure that N(w) is monotone +* + IWORK( JI ) = MIN( NAB( JI, 2 ), + $ MAX( NAB( JI, 1 ), IWORK( JI ) ) ) +* +* Update the Queue -- add intervals if both halves +* contain eigenvalues. +* + IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN +* +* No eigenvalue in the upper interval: +* just use the lower interval. +* + AB( JI, 2 ) = C( JI ) +* + ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN +* +* No eigenvalue in the lower interval: +* just use the upper interval. +* + AB( JI, 1 ) = C( JI ) + ELSE + KLNEW = KLNEW + 1 + IF( KLNEW.LE.MMAX ) THEN +* +* Eigenvalue in both intervals -- add upper to +* queue. +* + AB( KLNEW, 2 ) = AB( JI, 2 ) + NAB( KLNEW, 2 ) = NAB( JI, 2 ) + AB( KLNEW, 1 ) = C( JI ) + NAB( KLNEW, 1 ) = IWORK( JI ) + AB( JI, 2 ) = C( JI ) + NAB( JI, 2 ) = IWORK( JI ) + ELSE + INFO = MMAX + 1 + END IF + END IF + 70 CONTINUE + IF( INFO.NE.0 ) + $ RETURN + KL = KLNEW + ELSE +* +* IJOB=3: Binary search. Keep only the interval containing +* w s.t. N(w) = NVAL +* + DO 80 JI = KF, KL + IF( IWORK( JI ).LE.NVAL( JI ) ) THEN + AB( JI, 1 ) = C( JI ) + NAB( JI, 1 ) = IWORK( JI ) + END IF + IF( IWORK( JI ).GE.NVAL( JI ) ) THEN + AB( JI, 2 ) = C( JI ) + NAB( JI, 2 ) = IWORK( JI ) + END IF + 80 CONTINUE + END IF +* + ELSE +* +* End of Parallel Version of the loop +* +* Begin of Serial Version of the loop +* + KLNEW = KL + DO 100 JI = KF, KL +* +* Compute N(w), the number of eigenvalues less than w +* + TMP1 = C( JI ) + TMP2 = D( 1 ) - TMP1 + ITMP1 = 0 + IF( TMP2.LE.PIVMIN ) THEN + ITMP1 = 1 + TMP2 = MIN( TMP2, -PIVMIN ) + END IF +* +* A series of compiler directives to defeat vectorization +* for the next loop +* +*$PL$ CMCHAR=' ' +CDIR$ NEXTSCALAR +C$DIR SCALAR +CDIR$ NEXT SCALAR +CVD$L NOVECTOR +CDEC$ NOVECTOR +CVD$ NOVECTOR +*VDIR NOVECTOR +*VOCL LOOP,SCALAR +CIBM PREFER SCALAR +*$PL$ CMCHAR='*' +* + DO 90 J = 2, N + TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1 + IF( TMP2.LE.PIVMIN ) THEN + ITMP1 = ITMP1 + 1 + TMP2 = MIN( TMP2, -PIVMIN ) + END IF + 90 CONTINUE +* + IF( IJOB.LE.2 ) THEN +* +* IJOB=2: Choose all intervals containing eigenvalues. +* +* Insure that N(w) is monotone +* + ITMP1 = MIN( NAB( JI, 2 ), + $ MAX( NAB( JI, 1 ), ITMP1 ) ) +* +* Update the Queue -- add intervals if both halves +* contain eigenvalues. +* + IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN +* +* No eigenvalue in the upper interval: +* just use the lower interval. +* + AB( JI, 2 ) = TMP1 +* + ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN +* +* No eigenvalue in the lower interval: +* just use the upper interval. +* + AB( JI, 1 ) = TMP1 + ELSE IF( KLNEW.LT.MMAX ) THEN +* +* Eigenvalue in both intervals -- add upper to queue. +* + KLNEW = KLNEW + 1 + AB( KLNEW, 2 ) = AB( JI, 2 ) + NAB( KLNEW, 2 ) = NAB( JI, 2 ) + AB( KLNEW, 1 ) = TMP1 + NAB( KLNEW, 1 ) = ITMP1 + AB( JI, 2 ) = TMP1 + NAB( JI, 2 ) = ITMP1 + ELSE + INFO = MMAX + 1 + RETURN + END IF + ELSE +* +* IJOB=3: Binary search. Keep only the interval +* containing w s.t. N(w) = NVAL +* + IF( ITMP1.LE.NVAL( JI ) ) THEN + AB( JI, 1 ) = TMP1 + NAB( JI, 1 ) = ITMP1 + END IF + IF( ITMP1.GE.NVAL( JI ) ) THEN + AB( JI, 2 ) = TMP1 + NAB( JI, 2 ) = ITMP1 + END IF + END IF + 100 CONTINUE + KL = KLNEW +* +* End of Serial Version of the loop +* + END IF +* +* Check for convergence +* + KFNEW = KF + DO 110 JI = KF, KL + TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) ) + TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) ) + IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR. + $ NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN +* +* Converged -- Swap with position KFNEW, +* then increment KFNEW +* + IF( JI.GT.KFNEW ) THEN + TMP1 = AB( JI, 1 ) + TMP2 = AB( JI, 2 ) + ITMP1 = NAB( JI, 1 ) + ITMP2 = NAB( JI, 2 ) + AB( JI, 1 ) = AB( KFNEW, 1 ) + AB( JI, 2 ) = AB( KFNEW, 2 ) + NAB( JI, 1 ) = NAB( KFNEW, 1 ) + NAB( JI, 2 ) = NAB( KFNEW, 2 ) + AB( KFNEW, 1 ) = TMP1 + AB( KFNEW, 2 ) = TMP2 + NAB( KFNEW, 1 ) = ITMP1 + NAB( KFNEW, 2 ) = ITMP2 + IF( IJOB.EQ.3 ) THEN + ITMP1 = NVAL( JI ) + NVAL( JI ) = NVAL( KFNEW ) + NVAL( KFNEW ) = ITMP1 + END IF + END IF + KFNEW = KFNEW + 1 + END IF + 110 CONTINUE + KF = KFNEW +* +* Choose Midpoints +* + DO 120 JI = KF, KL + C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) ) + 120 CONTINUE +* +* If no more intervals to refine, quit. +* + IF( KF.GT.KL ) + $ GO TO 140 + 130 CONTINUE +* +* Converged +* + 140 CONTINUE + INFO = MAX( KL+1-KF, 0 ) + MOUT = KL +* + RETURN +* +* End of SLAEBZ +* + END |