1*> \brief <b> DSBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSBEVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbevx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbevx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbevx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
22*                          VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
23*                          IFAIL, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBZ, RANGE, UPLO
27*       INTEGER            IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
28*       DOUBLE PRECISION   ABSTOL, VL, VU
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IFAIL( * ), IWORK( * )
32*       DOUBLE PRECISION   AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
33*      $                   Z( LDZ, * )
34*       ..
35*
36*
37*> \par Purpose:
38*  =============
39*>
40*> \verbatim
41*>
42*> DSBEVX computes selected eigenvalues and, optionally, eigenvectors
43*> of a real symmetric band matrix A.  Eigenvalues and eigenvectors can
44*> be selected by specifying either a range of values or a range of
45*> indices for the desired eigenvalues.
46*> \endverbatim
47*
48*  Arguments:
49*  ==========
50*
51*> \param[in] JOBZ
52*> \verbatim
53*>          JOBZ is CHARACTER*1
54*>          = 'N':  Compute eigenvalues only;
55*>          = 'V':  Compute eigenvalues and eigenvectors.
56*> \endverbatim
57*>
58*> \param[in] RANGE
59*> \verbatim
60*>          RANGE is CHARACTER*1
61*>          = 'A': all eigenvalues will be found;
62*>          = 'V': all eigenvalues in the half-open interval (VL,VU]
63*>                 will be found;
64*>          = 'I': the IL-th through IU-th eigenvalues will be found.
65*> \endverbatim
66*>
67*> \param[in] UPLO
68*> \verbatim
69*>          UPLO is CHARACTER*1
70*>          = 'U':  Upper triangle of A is stored;
71*>          = 'L':  Lower triangle of A is stored.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*>          N is INTEGER
77*>          The order of the matrix A.  N >= 0.
78*> \endverbatim
79*>
80*> \param[in] KD
81*> \verbatim
82*>          KD is INTEGER
83*>          The number of superdiagonals of the matrix A if UPLO = 'U',
84*>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
85*> \endverbatim
86*>
87*> \param[in,out] AB
88*> \verbatim
89*>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
90*>          On entry, the upper or lower triangle of the symmetric band
91*>          matrix A, stored in the first KD+1 rows of the array.  The
92*>          j-th column of A is stored in the j-th column of the array AB
93*>          as follows:
94*>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
95*>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
96*>
97*>          On exit, AB is overwritten by values generated during the
98*>          reduction to tridiagonal form.  If UPLO = 'U', the first
99*>          superdiagonal and the diagonal of the tridiagonal matrix T
100*>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
101*>          the diagonal and first subdiagonal of T are returned in the
102*>          first two rows of AB.
103*> \endverbatim
104*>
105*> \param[in] LDAB
106*> \verbatim
107*>          LDAB is INTEGER
108*>          The leading dimension of the array AB.  LDAB >= KD + 1.
109*> \endverbatim
110*>
111*> \param[out] Q
112*> \verbatim
113*>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
114*>          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
115*>                         reduction to tridiagonal form.
116*>          If JOBZ = 'N', the array Q is not referenced.
117*> \endverbatim
118*>
119*> \param[in] LDQ
120*> \verbatim
121*>          LDQ is INTEGER
122*>          The leading dimension of the array Q.  If JOBZ = 'V', then
123*>          LDQ >= max(1,N).
124*> \endverbatim
125*>
126*> \param[in] VL
127*> \verbatim
128*>          VL is DOUBLE PRECISION
129*> \endverbatim
130*>
131*> \param[in] VU
132*> \verbatim
133*>          VU is DOUBLE PRECISION
134*>          If RANGE='V', the lower and upper bounds of the interval to
135*>          be searched for eigenvalues. VL < VU.
136*>          Not referenced if RANGE = 'A' or 'I'.
137*> \endverbatim
138*>
139*> \param[in] IL
140*> \verbatim
141*>          IL is INTEGER
142*> \endverbatim
143*>
144*> \param[in] IU
145*> \verbatim
146*>          IU is INTEGER
147*>          If RANGE='I', the indices (in ascending order) of the
148*>          smallest and largest eigenvalues to be returned.
149*>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
150*>          Not referenced if RANGE = 'A' or 'V'.
151*> \endverbatim
152*>
153*> \param[in] ABSTOL
154*> \verbatim
155*>          ABSTOL is DOUBLE PRECISION
156*>          The absolute error tolerance for the eigenvalues.
157*>          An approximate eigenvalue is accepted as converged
158*>          when it is determined to lie in an interval [a,b]
159*>          of width less than or equal to
160*>
161*>                  ABSTOL + EPS *   max( |a|,|b| ) ,
162*>
163*>          where EPS is the machine precision.  If ABSTOL is less than
164*>          or equal to zero, then  EPS*|T|  will be used in its place,
165*>          where |T| is the 1-norm of the tridiagonal matrix obtained
166*>          by reducing AB to tridiagonal form.
167*>
168*>          Eigenvalues will be computed most accurately when ABSTOL is
169*>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
170*>          If this routine returns with INFO>0, indicating that some
171*>          eigenvectors did not converge, try setting ABSTOL to
172*>          2*DLAMCH('S').
173*>
174*>          See "Computing Small Singular Values of Bidiagonal Matrices
175*>          with Guaranteed High Relative Accuracy," by Demmel and
176*>          Kahan, LAPACK Working Note #3.
177*> \endverbatim
178*>
179*> \param[out] M
180*> \verbatim
181*>          M is INTEGER
182*>          The total number of eigenvalues found.  0 <= M <= N.
183*>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
184*> \endverbatim
185*>
186*> \param[out] W
187*> \verbatim
188*>          W is DOUBLE PRECISION array, dimension (N)
189*>          The first M elements contain the selected eigenvalues in
190*>          ascending order.
191*> \endverbatim
192*>
193*> \param[out] Z
194*> \verbatim
195*>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
196*>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
197*>          contain the orthonormal eigenvectors of the matrix A
198*>          corresponding to the selected eigenvalues, with the i-th
199*>          column of Z holding the eigenvector associated with W(i).
200*>          If an eigenvector fails to converge, then that column of Z
201*>          contains the latest approximation to the eigenvector, and the
202*>          index of the eigenvector is returned in IFAIL.
203*>          If JOBZ = 'N', then Z is not referenced.
204*>          Note: the user must ensure that at least max(1,M) columns are
205*>          supplied in the array Z; if RANGE = 'V', the exact value of M
206*>          is not known in advance and an upper bound must be used.
207*> \endverbatim
208*>
209*> \param[in] LDZ
210*> \verbatim
211*>          LDZ is INTEGER
212*>          The leading dimension of the array Z.  LDZ >= 1, and if
213*>          JOBZ = 'V', LDZ >= max(1,N).
214*> \endverbatim
215*>
216*> \param[out] WORK
217*> \verbatim
218*>          WORK is DOUBLE PRECISION array, dimension (7*N)
219*> \endverbatim
220*>
221*> \param[out] IWORK
222*> \verbatim
223*>          IWORK is INTEGER array, dimension (5*N)
224*> \endverbatim
225*>
226*> \param[out] IFAIL
227*> \verbatim
228*>          IFAIL is INTEGER array, dimension (N)
229*>          If JOBZ = 'V', then if INFO = 0, the first M elements of
230*>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
231*>          indices of the eigenvectors that failed to converge.
232*>          If JOBZ = 'N', then IFAIL is not referenced.
233*> \endverbatim
234*>
235*> \param[out] INFO
236*> \verbatim
237*>          INFO is INTEGER
238*>          = 0:  successful exit.
239*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
240*>          > 0:  if INFO = i, then i eigenvectors failed to converge.
241*>                Their indices are stored in array IFAIL.
242*> \endverbatim
243*
244*  Authors:
245*  ========
246*
247*> \author Univ. of Tennessee
248*> \author Univ. of California Berkeley
249*> \author Univ. of Colorado Denver
250*> \author NAG Ltd.
251*
252*> \date November 2011
253*
254*> \ingroup doubleOTHEReigen
255*
256*  =====================================================================
257      SUBROUTINE DSBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL,
258     $                   VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
259     $                   IFAIL, INFO )
260*
261*  -- LAPACK driver routine (version 3.4.0) --
262*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
263*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
264*     November 2011
265*
266*     .. Scalar Arguments ..
267      CHARACTER          JOBZ, RANGE, UPLO
268      INTEGER            IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
269      DOUBLE PRECISION   ABSTOL, VL, VU
270*     ..
271*     .. Array Arguments ..
272      INTEGER            IFAIL( * ), IWORK( * )
273      DOUBLE PRECISION   AB( LDAB, * ), Q( LDQ, * ), W( * ), WORK( * ),
274     $                   Z( LDZ, * )
275*     ..
276*
277*  =====================================================================
278*
279*     .. Parameters ..
280      DOUBLE PRECISION   ZERO, ONE
281      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
282*     ..
283*     .. Local Scalars ..
284      LOGICAL            ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
285      CHARACTER          ORDER
286      INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
287     $                   INDISP, INDIWO, INDWRK, ISCALE, ITMP1, J, JJ,
288     $                   NSPLIT
289      DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
290     $                   SIGMA, SMLNUM, TMP1, VLL, VUU
291*     ..
292*     .. External Functions ..
293      LOGICAL            LSAME
294      DOUBLE PRECISION   DLAMCH, DLANSB
295      EXTERNAL           LSAME, DLAMCH, DLANSB
296*     ..
297*     .. External Subroutines ..
298      EXTERNAL           DCOPY, DGEMV, DLACPY, DLASCL, DSBTRD, DSCAL,
299     $                   DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA
300*     ..
301*     .. Intrinsic Functions ..
302      INTRINSIC          MAX, MIN, SQRT
303*     ..
304*     .. Executable Statements ..
305*
306*     Test the input parameters.
307*
308      WANTZ = LSAME( JOBZ, 'V' )
309      ALLEIG = LSAME( RANGE, 'A' )
310      VALEIG = LSAME( RANGE, 'V' )
311      INDEIG = LSAME( RANGE, 'I' )
312      LOWER = LSAME( UPLO, 'L' )
313*
314      INFO = 0
315      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
316         INFO = -1
317      ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
318         INFO = -2
319      ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
320         INFO = -3
321      ELSE IF( N.LT.0 ) THEN
322         INFO = -4
323      ELSE IF( KD.LT.0 ) THEN
324         INFO = -5
325      ELSE IF( LDAB.LT.KD+1 ) THEN
326         INFO = -7
327      ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN
328         INFO = -9
329      ELSE
330         IF( VALEIG ) THEN
331            IF( N.GT.0 .AND. VU.LE.VL )
332     $         INFO = -11
333         ELSE IF( INDEIG ) THEN
334            IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
335               INFO = -12
336            ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
337               INFO = -13
338            END IF
339         END IF
340      END IF
341      IF( INFO.EQ.0 ) THEN
342         IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
343     $      INFO = -18
344      END IF
345*
346      IF( INFO.NE.0 ) THEN
347         CALL XERBLA( 'DSBEVX', -INFO )
348         RETURN
349      END IF
350*
351*     Quick return if possible
352*
353      M = 0
354      IF( N.EQ.0 )
355     $   RETURN
356*
357      IF( N.EQ.1 ) THEN
358         M = 1
359         IF( LOWER ) THEN
360            TMP1 = AB( 1, 1 )
361         ELSE
362            TMP1 = AB( KD+1, 1 )
363         END IF
364         IF( VALEIG ) THEN
365            IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) )
366     $         M = 0
367         END IF
368         IF( M.EQ.1 ) THEN
369            W( 1 ) = TMP1
370            IF( WANTZ )
371     $         Z( 1, 1 ) = ONE
372         END IF
373         RETURN
374      END IF
375*
376*     Get machine constants.
377*
378      SAFMIN = DLAMCH( 'Safe minimum' )
379      EPS = DLAMCH( 'Precision' )
380      SMLNUM = SAFMIN / EPS
381      BIGNUM = ONE / SMLNUM
382      RMIN = SQRT( SMLNUM )
383      RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
384*
385*     Scale matrix to allowable range, if necessary.
386*
387      ISCALE = 0
388      ABSTLL = ABSTOL
389      IF( VALEIG ) THEN
390         VLL = VL
391         VUU = VU
392      ELSE
393         VLL = ZERO
394         VUU = ZERO
395      END IF
396      ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
397      IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
398         ISCALE = 1
399         SIGMA = RMIN / ANRM
400      ELSE IF( ANRM.GT.RMAX ) THEN
401         ISCALE = 1
402         SIGMA = RMAX / ANRM
403      END IF
404      IF( ISCALE.EQ.1 ) THEN
405         IF( LOWER ) THEN
406            CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
407         ELSE
408            CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
409         END IF
410         IF( ABSTOL.GT.0 )
411     $      ABSTLL = ABSTOL*SIGMA
412         IF( VALEIG ) THEN
413            VLL = VL*SIGMA
414            VUU = VU*SIGMA
415         END IF
416      END IF
417*
418*     Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
419*
420      INDD = 1
421      INDE = INDD + N
422      INDWRK = INDE + N
423      CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, WORK( INDD ),
424     $             WORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO )
425*
426*     If all eigenvalues are desired and ABSTOL is less than or equal
427*     to zero, then call DSTERF or SSTEQR.  If this fails for some
428*     eigenvalue, then try DSTEBZ.
429*
430      TEST = .FALSE.
431      IF (INDEIG) THEN
432         IF (IL.EQ.1 .AND. IU.EQ.N) THEN
433            TEST = .TRUE.
434         END IF
435      END IF
436      IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN
437         CALL DCOPY( N, WORK( INDD ), 1, W, 1 )
438         INDEE = INDWRK + 2*N
439         IF( .NOT.WANTZ ) THEN
440            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
441            CALL DSTERF( N, W, WORK( INDEE ), INFO )
442         ELSE
443            CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
444            CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 )
445            CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ,
446     $                   WORK( INDWRK ), INFO )
447            IF( INFO.EQ.0 ) THEN
448               DO 10 I = 1, N
449                  IFAIL( I ) = 0
450   10          CONTINUE
451            END IF
452         END IF
453         IF( INFO.EQ.0 ) THEN
454            M = N
455            GO TO 30
456         END IF
457         INFO = 0
458      END IF
459*
460*     Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
461*
462      IF( WANTZ ) THEN
463         ORDER = 'B'
464      ELSE
465         ORDER = 'E'
466      END IF
467      INDIBL = 1
468      INDISP = INDIBL + N
469      INDIWO = INDISP + N
470      CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
471     $             WORK( INDD ), WORK( INDE ), M, NSPLIT, W,
472     $             IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ),
473     $             IWORK( INDIWO ), INFO )
474*
475      IF( WANTZ ) THEN
476         CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W,
477     $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
478     $                WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO )
479*
480*        Apply orthogonal matrix used in reduction to tridiagonal
481*        form to eigenvectors returned by DSTEIN.
482*
483         DO 20 J = 1, M
484            CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
485            CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO,
486     $                  Z( 1, J ), 1 )
487   20    CONTINUE
488      END IF
489*
490*     If matrix was scaled, then rescale eigenvalues appropriately.
491*
492   30 CONTINUE
493      IF( ISCALE.EQ.1 ) THEN
494         IF( INFO.EQ.0 ) THEN
495            IMAX = M
496         ELSE
497            IMAX = INFO - 1
498         END IF
499         CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
500      END IF
501*
502*     If eigenvalues are not in order, then sort them, along with
503*     eigenvectors.
504*
505      IF( WANTZ ) THEN
506         DO 50 J = 1, M - 1
507            I = 0
508            TMP1 = W( J )
509            DO 40 JJ = J + 1, M
510               IF( W( JJ ).LT.TMP1 ) THEN
511                  I = JJ
512                  TMP1 = W( JJ )
513               END IF
514   40       CONTINUE
515*
516            IF( I.NE.0 ) THEN
517               ITMP1 = IWORK( INDIBL+I-1 )
518               W( I ) = W( J )
519               IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
520               W( J ) = TMP1
521               IWORK( INDIBL+J-1 ) = ITMP1
522               CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
523               IF( INFO.NE.0 ) THEN
524                  ITMP1 = IFAIL( I )
525                  IFAIL( I ) = IFAIL( J )
526                  IFAIL( J ) = ITMP1
527               END IF
528            END IF
529   50    CONTINUE
530      END IF
531*
532      RETURN
533*
534*     End of DSBEVX
535*
536      END
537