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