1*> \brief \b ZHBGVD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHBGVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgvd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgvd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgvd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
22*                          Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
23*                          LIWORK, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBZ, UPLO
27*       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
28*      $                   LWORK, N
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IWORK( * )
32*       DOUBLE PRECISION   RWORK( * ), W( * )
33*       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
34*      $                   Z( LDZ, * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> ZHBGVD computes all the eigenvalues, and optionally, the eigenvectors
44*> of a complex generalized Hermitian-definite banded eigenproblem, of
45*> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
46*> and banded, and B is also positive definite.  If eigenvectors are
47*> desired, it uses a divide and conquer algorithm.
48*>
49*> The divide and conquer algorithm makes very mild assumptions about
50*> floating point arithmetic. It will work on machines with a guard
51*> digit in add/subtract, or on those binary machines without guard
52*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
53*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
54*> without guard digits, but we know of none.
55*> \endverbatim
56*
57*  Arguments:
58*  ==========
59*
60*> \param[in] JOBZ
61*> \verbatim
62*>          JOBZ is CHARACTER*1
63*>          = 'N':  Compute eigenvalues only;
64*>          = 'V':  Compute eigenvalues and eigenvectors.
65*> \endverbatim
66*>
67*> \param[in] UPLO
68*> \verbatim
69*>          UPLO is CHARACTER*1
70*>          = 'U':  Upper triangles of A and B are stored;
71*>          = 'L':  Lower triangles of A and B are stored.
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*>          N is INTEGER
77*>          The order of the matrices A and B.  N >= 0.
78*> \endverbatim
79*>
80*> \param[in] KA
81*> \verbatim
82*>          KA is INTEGER
83*>          The number of superdiagonals of the matrix A if UPLO = 'U',
84*>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
85*> \endverbatim
86*>
87*> \param[in] KB
88*> \verbatim
89*>          KB is INTEGER
90*>          The number of superdiagonals of the matrix B if UPLO = 'U',
91*>          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
92*> \endverbatim
93*>
94*> \param[in,out] AB
95*> \verbatim
96*>          AB is COMPLEX*16 array, dimension (LDAB, N)
97*>          On entry, the upper or lower triangle of the Hermitian band
98*>          matrix A, stored in the first ka+1 rows of the array.  The
99*>          j-th column of A is stored in the j-th column of the array AB
100*>          as follows:
101*>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
102*>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
103*>
104*>          On exit, the contents of AB are destroyed.
105*> \endverbatim
106*>
107*> \param[in] LDAB
108*> \verbatim
109*>          LDAB is INTEGER
110*>          The leading dimension of the array AB.  LDAB >= KA+1.
111*> \endverbatim
112*>
113*> \param[in,out] BB
114*> \verbatim
115*>          BB is COMPLEX*16 array, dimension (LDBB, N)
116*>          On entry, the upper or lower triangle of the Hermitian band
117*>          matrix B, stored in the first kb+1 rows of the array.  The
118*>          j-th column of B is stored in the j-th column of the array BB
119*>          as follows:
120*>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
121*>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
122*>
123*>          On exit, the factor S from the split Cholesky factorization
124*>          B = S**H*S, as returned by ZPBSTF.
125*> \endverbatim
126*>
127*> \param[in] LDBB
128*> \verbatim
129*>          LDBB is INTEGER
130*>          The leading dimension of the array BB.  LDBB >= KB+1.
131*> \endverbatim
132*>
133*> \param[out] W
134*> \verbatim
135*>          W is DOUBLE PRECISION array, dimension (N)
136*>          If INFO = 0, the eigenvalues in ascending order.
137*> \endverbatim
138*>
139*> \param[out] Z
140*> \verbatim
141*>          Z is COMPLEX*16 array, dimension (LDZ, N)
142*>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
143*>          eigenvectors, with the i-th column of Z holding the
144*>          eigenvector associated with W(i). The eigenvectors are
145*>          normalized so that Z**H*B*Z = I.
146*>          If JOBZ = 'N', then Z is not referenced.
147*> \endverbatim
148*>
149*> \param[in] LDZ
150*> \verbatim
151*>          LDZ is INTEGER
152*>          The leading dimension of the array Z.  LDZ >= 1, and if
153*>          JOBZ = 'V', LDZ >= N.
154*> \endverbatim
155*>
156*> \param[out] WORK
157*> \verbatim
158*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
159*>          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
160*> \endverbatim
161*>
162*> \param[in] LWORK
163*> \verbatim
164*>          LWORK is INTEGER
165*>          The dimension of the array WORK.
166*>          If N <= 1,               LWORK >= 1.
167*>          If JOBZ = 'N' and N > 1, LWORK >= N.
168*>          If JOBZ = 'V' and N > 1, LWORK >= 2*N**2.
169*>
170*>          If LWORK = -1, then a workspace query is assumed; the routine
171*>          only calculates the optimal sizes of the WORK, RWORK and
172*>          IWORK arrays, returns these values as the first entries of
173*>          the WORK, RWORK and IWORK arrays, and no error message
174*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
175*> \endverbatim
176*>
177*> \param[out] RWORK
178*> \verbatim
179*>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
180*>          On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
181*> \endverbatim
182*>
183*> \param[in] LRWORK
184*> \verbatim
185*>          LRWORK is INTEGER
186*>          The dimension of array RWORK.
187*>          If N <= 1,               LRWORK >= 1.
188*>          If JOBZ = 'N' and N > 1, LRWORK >= N.
189*>          If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
190*>
191*>          If LRWORK = -1, then a workspace query is assumed; the
192*>          routine only calculates the optimal sizes of the WORK, RWORK
193*>          and IWORK arrays, returns these values as the first entries
194*>          of the WORK, RWORK and IWORK arrays, and no error message
195*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
196*> \endverbatim
197*>
198*> \param[out] IWORK
199*> \verbatim
200*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
201*>          On exit, if INFO=0, IWORK(1) returns the optimal LIWORK.
202*> \endverbatim
203*>
204*> \param[in] LIWORK
205*> \verbatim
206*>          LIWORK is INTEGER
207*>          The dimension of array IWORK.
208*>          If JOBZ = 'N' or N <= 1, LIWORK >= 1.
209*>          If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
210*>
211*>          If LIWORK = -1, then a workspace query is assumed; the
212*>          routine only calculates the optimal sizes of the WORK, RWORK
213*>          and IWORK arrays, returns these values as the first entries
214*>          of the WORK, RWORK and IWORK arrays, and no error message
215*>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
216*> \endverbatim
217*>
218*> \param[out] INFO
219*> \verbatim
220*>          INFO is INTEGER
221*>          = 0:  successful exit
222*>          < 0:  if INFO = -i, the i-th argument had an illegal value
223*>          > 0:  if INFO = i, and i is:
224*>             <= N:  the algorithm failed to converge:
225*>                    i off-diagonal elements of an intermediate
226*>                    tridiagonal form did not converge to zero;
227*>             > N:   if INFO = N + i, for 1 <= i <= N, then ZPBSTF
228*>                    returned INFO = i: B is not positive definite.
229*>                    The factorization of B could not be completed and
230*>                    no eigenvalues or eigenvectors were computed.
231*> \endverbatim
232*
233*  Authors:
234*  ========
235*
236*> \author Univ. of Tennessee
237*> \author Univ. of California Berkeley
238*> \author Univ. of Colorado Denver
239*> \author NAG Ltd.
240*
241*> \date November 2015
242*
243*> \ingroup complex16OTHEReigen
244*
245*> \par Contributors:
246*  ==================
247*>
248*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
249*
250*  =====================================================================
251      SUBROUTINE ZHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
252     $                   Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
253     $                   LIWORK, INFO )
254*
255*  -- LAPACK driver routine (version 3.6.0) --
256*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
257*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
258*     November 2015
259*
260*     .. Scalar Arguments ..
261      CHARACTER          JOBZ, UPLO
262      INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
263     $                   LWORK, N
264*     ..
265*     .. Array Arguments ..
266      INTEGER            IWORK( * )
267      DOUBLE PRECISION   RWORK( * ), W( * )
268      COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
269     $                   Z( LDZ, * )
270*     ..
271*
272*  =====================================================================
273*
274*     .. Parameters ..
275      COMPLEX*16         CONE, CZERO
276      PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
277     $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
278*     ..
279*     .. Local Scalars ..
280      LOGICAL            LQUERY, UPPER, WANTZ
281      CHARACTER          VECT
282      INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK,
283     $                   LLWK2, LRWMIN, LWMIN
284*     ..
285*     .. External Functions ..
286      LOGICAL            LSAME
287      EXTERNAL           LSAME
288*     ..
289*     .. External Subroutines ..
290      EXTERNAL           DSTERF, XERBLA, ZGEMM, ZHBGST, ZHBTRD, ZLACPY,
291     $                   ZPBSTF, ZSTEDC
292*     ..
293*     .. Executable Statements ..
294*
295*     Test the input parameters.
296*
297      WANTZ = LSAME( JOBZ, 'V' )
298      UPPER = LSAME( UPLO, 'U' )
299      LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
300*
301      INFO = 0
302      IF( N.LE.1 ) THEN
303         LWMIN = 1+N
304         LRWMIN = 1+N
305         LIWMIN = 1
306      ELSE IF( WANTZ ) THEN
307         LWMIN = 2*N**2
308         LRWMIN = 1 + 5*N + 2*N**2
309         LIWMIN = 3 + 5*N
310      ELSE
311         LWMIN = N
312         LRWMIN = N
313         LIWMIN = 1
314      END IF
315      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
316         INFO = -1
317      ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
318         INFO = -2
319      ELSE IF( N.LT.0 ) THEN
320         INFO = -3
321      ELSE IF( KA.LT.0 ) THEN
322         INFO = -4
323      ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
324         INFO = -5
325      ELSE IF( LDAB.LT.KA+1 ) THEN
326         INFO = -7
327      ELSE IF( LDBB.LT.KB+1 ) THEN
328         INFO = -9
329      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
330         INFO = -12
331      END IF
332*
333      IF( INFO.EQ.0 ) THEN
334         WORK( 1 ) = LWMIN
335         RWORK( 1 ) = LRWMIN
336         IWORK( 1 ) = LIWMIN
337*
338         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
339            INFO = -14
340         ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
341            INFO = -16
342         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
343            INFO = -18
344         END IF
345      END IF
346*
347      IF( INFO.NE.0 ) THEN
348         CALL XERBLA( 'ZHBGVD', -INFO )
349         RETURN
350      ELSE IF( LQUERY ) THEN
351         RETURN
352      END IF
353*
354*     Quick return if possible
355*
356      IF( N.EQ.0 )
357     $   RETURN
358*
359*     Form a split Cholesky factorization of B.
360*
361      CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO )
362      IF( INFO.NE.0 ) THEN
363         INFO = N + INFO
364         RETURN
365      END IF
366*
367*     Transform problem to standard eigenvalue problem.
368*
369      INDE = 1
370      INDWRK = INDE + N
371      INDWK2 = 1 + N*N
372      LLWK2 = LWORK - INDWK2 + 2
373      LLRWK = LRWORK - INDWRK + 2
374      CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
375     $             WORK, RWORK( INDWRK ), IINFO )
376*
377*     Reduce Hermitian band matrix to tridiagonal form.
378*
379      IF( WANTZ ) THEN
380         VECT = 'U'
381      ELSE
382         VECT = 'N'
383      END IF
384      CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z,
385     $             LDZ, WORK, IINFO )
386*
387*     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
388*
389      IF( .NOT.WANTZ ) THEN
390         CALL DSTERF( N, W, RWORK( INDE ), INFO )
391      ELSE
392         CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
393     $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
394     $                INFO )
395         CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
396     $               WORK( INDWK2 ), N )
397         CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
398      END IF
399*
400      WORK( 1 ) = LWMIN
401      RWORK( 1 ) = LRWMIN
402      IWORK( 1 ) = LIWMIN
403      RETURN
404*
405*     End of ZHBGVD
406*
407      END
408