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