--- rpl/lapack/lapack/dlaebz.f 2010/08/06 15:28:38 1.3 +++ rpl/lapack/lapack/dlaebz.f 2012/08/22 09:48:16 1.11 @@ -1,11 +1,328 @@ +*> \brief \b DLAEBZ +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLAEBZ + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, +* RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, +* NAB, WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX +* DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * ) +* DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ), +* $ WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAEBZ 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. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] IJOB +*> \verbatim +*> 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 DLAEBZ to return with INFO=-1. +*> \endverbatim +*> +*> \param[in] NITMAX +*> \verbatim +*> NITMAX is 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. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The dimension n of the tridiagonal matrix T. It must be at +*> least 1. +*> \endverbatim +*> +*> \param[in] MMAX +*> \verbatim +*> MMAX is INTEGER +*> The maximum number of intervals. If more than MMAX intervals +*> are generated, then DLAEBZ will quit with INFO=MMAX+1. +*> \endverbatim +*> +*> \param[in] MINP +*> \verbatim +*> MINP is INTEGER +*> The initial number of intervals. It may not be greater than +*> MMAX. +*> \endverbatim +*> +*> \param[in] NBMIN +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in] ABSTOL +*> \verbatim +*> ABSTOL is DOUBLE PRECISION +*> 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. +*> \endverbatim +*> +*> \param[in] RELTOL +*> \verbatim +*> RELTOL is DOUBLE PRECISION +*> 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. +*> \endverbatim +*> +*> \param[in] PIVMIN +*> \verbatim +*> PIVMIN is DOUBLE PRECISION +*> 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. +*> \endverbatim +*> +*> \param[in] D +*> \verbatim +*> D is DOUBLE PRECISION array, dimension (N) +*> The diagonal elements of the tridiagonal matrix T. +*> \endverbatim +*> +*> \param[in] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N) +*> The offdiagonal elements of the tridiagonal matrix T in +*> positions 1 through N-1. E(N) is arbitrary. +*> \endverbatim +*> +*> \param[in] E2 +*> \verbatim +*> E2 is DOUBLE PRECISION array, dimension (N) +*> The squares of the offdiagonal elements of the tridiagonal +*> matrix T. E2(N) is ignored. +*> \endverbatim +*> +*> \param[in,out] NVAL +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in,out] AB +*> \verbatim +*> AB is DOUBLE PRECISION 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. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is DOUBLE PRECISION 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. +*> \endverbatim +*> +*> \param[out] MOUT +*> \verbatim +*> 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. +*> \endverbatim +*> +*> \param[in,out] NAB +*> \verbatim +*> 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 DLAEBZ 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 DLAEBZ is called. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MMAX) +*> Workspace. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (MMAX) +*> Workspace. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: All intervals converged. +*> = 1--MMAX: The last INFO intervals did not converge. +*> = MMAX+1: More than MMAX intervals were generated. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup auxOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> 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, DLAEBZ should have one or +*> more initial intervals set up in AB, and DLAEBZ 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. DLAEBZ 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). DLAEBZ 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. +*> \endverbatim +*> +* ===================================================================== SUBROUTINE DLAEBZ( 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.2) -- +* -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* November 2011 * * .. Scalar Arguments .. INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX @@ -17,208 +334,6 @@ $ WORK( * ) * .. * -* Purpose -* ======= -* -* DLAEBZ 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 DLAEBZ 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 DLAEBZ 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) DOUBLE PRECISION -* 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) DOUBLE PRECISION -* 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) DOUBLE PRECISION -* 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) DOUBLE PRECISION array, dimension (N) -* The diagonal elements of the tridiagonal matrix T. -* -* E (input) DOUBLE PRECISION array, dimension (N) -* The offdiagonal elements of the tridiagonal matrix T in -* positions 1 through N-1. E(N) is arbitrary. -* -* E2 (input) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DLAEBZ 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 DLAEBZ is called. -* -* WORK (workspace) DOUBLE PRECISION 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, DLAEBZ should have one or -* more initial intervals set up in AB, and DLAEBZ 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. DLAEBZ 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). DLAEBZ 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 .. @@ -251,7 +366,6 @@ * Compute the number of eigenvalues in the initial intervals. * MOUT = 0 -*DIR$ NOVECTOR DO 30 JI = 1, MINP DO 20 JP = 1, 2 TMP1 = D( 1 ) - AB( JI, JP ) @@ -407,21 +521,6 @@ 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 @@ -487,8 +586,6 @@ CIBM PREFER SCALAR 100 CONTINUE KL = KLNEW * -* End of Serial Version of the loop -* END IF * * Check for convergence