Annotation of rpl/lapack/lapack/zhbgvd.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZHBGST
! 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 2011
! 242: *
! 243: *> \ingroup complex16OTHEReigen
! 244: *
! 245: *> \par Contributors:
! 246: * ==================
! 247: *>
! 248: *> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
! 249: *
! 250: * =====================================================================
1.1 bertrand 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: *
1.9 ! bertrand 255: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 256: * -- LAPACK is a software package provided by Univ. of Tennessee, --
257: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 258: * November 2011
1.1 bertrand 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
1.8 bertrand 303: LWMIN = 1+N
304: LRWMIN = 1+N
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>