001:       SUBROUTINE ZHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
002:      $                   Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK,
003:      $                   LIWORK, INFO )
004: *
005: *  -- LAPACK driver routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          JOBZ, UPLO
012:       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK,
013:      $                   LWORK, N
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IWORK( * )
017:       DOUBLE PRECISION   RWORK( * ), W( * )
018:       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
019:      $                   Z( LDZ, * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZHBGVD computes all the eigenvalues, and optionally, the eigenvectors
026: *  of a complex generalized Hermitian-definite banded eigenproblem, of
027: *  the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
028: *  and banded, and B is also positive definite.  If eigenvectors are
029: *  desired, it uses a divide and conquer algorithm.
030: *
031: *  The divide and conquer algorithm makes very mild assumptions about
032: *  floating point arithmetic. It will work on machines with a guard
033: *  digit in add/subtract, or on those binary machines without guard
034: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
035: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
036: *  without guard digits, but we know of none.
037: *
038: *  Arguments
039: *  =========
040: *
041: *  JOBZ    (input) CHARACTER*1
042: *          = 'N':  Compute eigenvalues only;
043: *          = 'V':  Compute eigenvalues and eigenvectors.
044: *
045: *  UPLO    (input) CHARACTER*1
046: *          = 'U':  Upper triangles of A and B are stored;
047: *          = 'L':  Lower triangles of A and B are stored.
048: *
049: *  N       (input) INTEGER
050: *          The order of the matrices A and B.  N >= 0.
051: *
052: *  KA      (input) INTEGER
053: *          The number of superdiagonals of the matrix A if UPLO = 'U',
054: *          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
055: *
056: *  KB      (input) INTEGER
057: *          The number of superdiagonals of the matrix B if UPLO = 'U',
058: *          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
059: *
060: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB, N)
061: *          On entry, the upper or lower triangle of the Hermitian band
062: *          matrix A, stored in the first ka+1 rows of the array.  The
063: *          j-th column of A is stored in the j-th column of the array AB
064: *          as follows:
065: *          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
066: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
067: *
068: *          On exit, the contents of AB are destroyed.
069: *
070: *  LDAB    (input) INTEGER
071: *          The leading dimension of the array AB.  LDAB >= KA+1.
072: *
073: *  BB      (input/output) COMPLEX*16 array, dimension (LDBB, N)
074: *          On entry, the upper or lower triangle of the Hermitian band
075: *          matrix B, stored in the first kb+1 rows of the array.  The
076: *          j-th column of B is stored in the j-th column of the array BB
077: *          as follows:
078: *          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
079: *          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
080: *
081: *          On exit, the factor S from the split Cholesky factorization
082: *          B = S**H*S, as returned by ZPBSTF.
083: *
084: *  LDBB    (input) INTEGER
085: *          The leading dimension of the array BB.  LDBB >= KB+1.
086: *
087: *  W       (output) DOUBLE PRECISION array, dimension (N)
088: *          If INFO = 0, the eigenvalues in ascending order.
089: *
090: *  Z       (output) COMPLEX*16 array, dimension (LDZ, N)
091: *          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
092: *          eigenvectors, with the i-th column of Z holding the
093: *          eigenvector associated with W(i). The eigenvectors are
094: *          normalized so that Z**H*B*Z = I.
095: *          If JOBZ = 'N', then Z is not referenced.
096: *
097: *  LDZ     (input) INTEGER
098: *          The leading dimension of the array Z.  LDZ >= 1, and if
099: *          JOBZ = 'V', LDZ >= N.
100: *
101: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
102: *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
103: *
104: *  LWORK   (input) INTEGER
105: *          The dimension of the array WORK.
106: *          If N <= 1,               LWORK >= 1.
107: *          If JOBZ = 'N' and N > 1, LWORK >= N.
108: *          If JOBZ = 'V' and N > 1, LWORK >= 2*N**2.
109: *
110: *          If LWORK = -1, then a workspace query is assumed; the routine
111: *          only calculates the optimal sizes of the WORK, RWORK and
112: *          IWORK arrays, returns these values as the first entries of
113: *          the WORK, RWORK and IWORK arrays, and no error message
114: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
115: *
116: *  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
117: *          On exit, if INFO=0, RWORK(1) returns the optimal LRWORK.
118: *
119: *  LRWORK  (input) INTEGER
120: *          The dimension of array RWORK.
121: *          If N <= 1,               LRWORK >= 1.
122: *          If JOBZ = 'N' and N > 1, LRWORK >= N.
123: *          If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
124: *
125: *          If LRWORK = -1, then a workspace query is assumed; the
126: *          routine only calculates the optimal sizes of the WORK, RWORK
127: *          and IWORK arrays, returns these values as the first entries
128: *          of the WORK, RWORK and IWORK arrays, and no error message
129: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
130: *
131: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
132: *          On exit, if INFO=0, IWORK(1) returns the optimal LIWORK.
133: *
134: *  LIWORK  (input) INTEGER
135: *          The dimension of array IWORK.
136: *          If JOBZ = 'N' or N <= 1, LIWORK >= 1.
137: *          If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
138: *
139: *          If LIWORK = -1, then a workspace query is assumed; the
140: *          routine only calculates the optimal sizes of the WORK, RWORK
141: *          and IWORK arrays, returns these values as the first entries
142: *          of the WORK, RWORK and IWORK arrays, and no error message
143: *          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
144: *
145: *  INFO    (output) INTEGER
146: *          = 0:  successful exit
147: *          < 0:  if INFO = -i, the i-th argument had an illegal value
148: *          > 0:  if INFO = i, and i is:
149: *             <= N:  the algorithm failed to converge:
150: *                    i off-diagonal elements of an intermediate
151: *                    tridiagonal form did not converge to zero;
152: *             > N:   if INFO = N + i, for 1 <= i <= N, then ZPBSTF
153: *                    returned INFO = i: B is not positive definite.
154: *                    The factorization of B could not be completed and
155: *                    no eigenvalues or eigenvectors were computed.
156: *
157: *  Further Details
158: *  ===============
159: *
160: *  Based on contributions by
161: *     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
162: *
163: *  =====================================================================
164: *
165: *     .. Parameters ..
166:       COMPLEX*16         CONE, CZERO
167:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ),
168:      $                   CZERO = ( 0.0D+0, 0.0D+0 ) )
169: *     ..
170: *     .. Local Scalars ..
171:       LOGICAL            LQUERY, UPPER, WANTZ
172:       CHARACTER          VECT
173:       INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK,
174:      $                   LLWK2, LRWMIN, LWMIN
175: *     ..
176: *     .. External Functions ..
177:       LOGICAL            LSAME
178:       EXTERNAL           LSAME
179: *     ..
180: *     .. External Subroutines ..
181:       EXTERNAL           DSTERF, XERBLA, ZGEMM, ZHBGST, ZHBTRD, ZLACPY,
182:      $                   ZPBSTF, ZSTEDC
183: *     ..
184: *     .. Executable Statements ..
185: *
186: *     Test the input parameters.
187: *
188:       WANTZ = LSAME( JOBZ, 'V' )
189:       UPPER = LSAME( UPLO, 'U' )
190:       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
191: *
192:       INFO = 0
193:       IF( N.LE.1 ) THEN
194:          LWMIN = 1
195:          LRWMIN = 1
196:          LIWMIN = 1
197:       ELSE IF( WANTZ ) THEN
198:          LWMIN = 2*N**2
199:          LRWMIN = 1 + 5*N + 2*N**2
200:          LIWMIN = 3 + 5*N
201:       ELSE
202:          LWMIN = N
203:          LRWMIN = N
204:          LIWMIN = 1
205:       END IF
206:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
207:          INFO = -1
208:       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
209:          INFO = -2
210:       ELSE IF( N.LT.0 ) THEN
211:          INFO = -3
212:       ELSE IF( KA.LT.0 ) THEN
213:          INFO = -4
214:       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
215:          INFO = -5
216:       ELSE IF( LDAB.LT.KA+1 ) THEN
217:          INFO = -7
218:       ELSE IF( LDBB.LT.KB+1 ) THEN
219:          INFO = -9
220:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
221:          INFO = -12
222:       END IF
223: *
224:       IF( INFO.EQ.0 ) THEN
225:          WORK( 1 ) = LWMIN
226:          RWORK( 1 ) = LRWMIN
227:          IWORK( 1 ) = LIWMIN
228: *
229:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
230:             INFO = -14
231:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
232:             INFO = -16
233:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
234:             INFO = -18
235:          END IF
236:       END IF
237: *
238:       IF( INFO.NE.0 ) THEN
239:          CALL XERBLA( 'ZHBGVD', -INFO )
240:          RETURN
241:       ELSE IF( LQUERY ) THEN
242:          RETURN
243:       END IF
244: *
245: *     Quick return if possible
246: *
247:       IF( N.EQ.0 )
248:      $   RETURN
249: *
250: *     Form a split Cholesky factorization of B.
251: *
252:       CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO )
253:       IF( INFO.NE.0 ) THEN
254:          INFO = N + INFO
255:          RETURN
256:       END IF
257: *
258: *     Transform problem to standard eigenvalue problem.
259: *
260:       INDE = 1
261:       INDWRK = INDE + N
262:       INDWK2 = 1 + N*N
263:       LLWK2 = LWORK - INDWK2 + 2
264:       LLRWK = LRWORK - INDWRK + 2
265:       CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
266:      $             WORK, RWORK( INDWRK ), IINFO )
267: *
268: *     Reduce Hermitian band matrix to tridiagonal form.
269: *
270:       IF( WANTZ ) THEN
271:          VECT = 'U'
272:       ELSE
273:          VECT = 'N'
274:       END IF
275:       CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z,
276:      $             LDZ, WORK, IINFO )
277: *
278: *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
279: *
280:       IF( .NOT.WANTZ ) THEN
281:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
282:       ELSE
283:          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
284:      $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
285:      $                INFO )
286:          CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
287:      $               WORK( INDWK2 ), N )
288:          CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
289:       END IF
290: *
291:       WORK( 1 ) = LWMIN
292:       RWORK( 1 ) = LRWMIN
293:       IWORK( 1 ) = LIWMIN
294:       RETURN
295: *
296: *     End of ZHBGVD
297: *
298:       END
299: