1: SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
2: $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER JOBZ, UPLO
11: INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: INTEGER IWORK( * )
15: DOUBLE PRECISION RWORK( * ), W( * )
16: COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZHBEVD computes all the eigenvalues and, optionally, eigenvectors of
23: * a complex Hermitian band matrix A. If eigenvectors are desired, it
24: * uses a divide and conquer algorithm.
25: *
26: * The divide and conquer algorithm makes very mild assumptions about
27: * floating point arithmetic. It will work on machines with a guard
28: * digit in add/subtract, or on those binary machines without guard
29: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
30: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
31: * without guard digits, but we know of none.
32: *
33: * Arguments
34: * =========
35: *
36: * JOBZ (input) CHARACTER*1
37: * = 'N': Compute eigenvalues only;
38: * = 'V': Compute eigenvalues and eigenvectors.
39: *
40: * UPLO (input) CHARACTER*1
41: * = 'U': Upper triangle of A is stored;
42: * = 'L': Lower triangle of A is stored.
43: *
44: * N (input) INTEGER
45: * The order of the matrix A. N >= 0.
46: *
47: * KD (input) INTEGER
48: * The number of superdiagonals of the matrix A if UPLO = 'U',
49: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
50: *
51: * AB (input/output) COMPLEX*16 array, dimension (LDAB, N)
52: * On entry, the upper or lower triangle of the Hermitian band
53: * matrix A, stored in the first KD+1 rows of the array. The
54: * j-th column of A is stored in the j-th column of the array AB
55: * as follows:
56: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
57: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
58: *
59: * On exit, AB is overwritten by values generated during the
60: * reduction to tridiagonal form. If UPLO = 'U', the first
61: * superdiagonal and the diagonal of the tridiagonal matrix T
62: * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
63: * the diagonal and first subdiagonal of T are returned in the
64: * first two rows of AB.
65: *
66: * LDAB (input) INTEGER
67: * The leading dimension of the array AB. LDAB >= KD + 1.
68: *
69: * W (output) DOUBLE PRECISION array, dimension (N)
70: * If INFO = 0, the eigenvalues in ascending order.
71: *
72: * Z (output) COMPLEX*16 array, dimension (LDZ, N)
73: * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
74: * eigenvectors of the matrix A, with the i-th column of Z
75: * holding the eigenvector associated with W(i).
76: * If JOBZ = 'N', then Z is not referenced.
77: *
78: * LDZ (input) INTEGER
79: * The leading dimension of the array Z. LDZ >= 1, and if
80: * JOBZ = 'V', LDZ >= max(1,N).
81: *
82: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
83: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84: *
85: * LWORK (input) INTEGER
86: * The dimension of the array WORK.
87: * If N <= 1, LWORK must be at least 1.
88: * If JOBZ = 'N' and N > 1, LWORK must be at least N.
89: * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2.
90: *
91: * If LWORK = -1, then a workspace query is assumed; the routine
92: * only calculates the optimal sizes of the WORK, RWORK and
93: * IWORK arrays, returns these values as the first entries of
94: * the WORK, RWORK and IWORK arrays, and no error message
95: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
96: *
97: * RWORK (workspace/output) DOUBLE PRECISION array,
98: * dimension (LRWORK)
99: * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
100: *
101: * LRWORK (input) INTEGER
102: * The dimension of array RWORK.
103: * If N <= 1, LRWORK must be at least 1.
104: * If JOBZ = 'N' and N > 1, LRWORK must be at least N.
105: * If JOBZ = 'V' and N > 1, LRWORK must be at least
106: * 1 + 5*N + 2*N**2.
107: *
108: * If LRWORK = -1, then a workspace query is assumed; the
109: * routine only calculates the optimal sizes of the WORK, RWORK
110: * and IWORK arrays, returns these values as the first entries
111: * of the WORK, RWORK and IWORK arrays, and no error message
112: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
113: *
114: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
115: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
116: *
117: * LIWORK (input) INTEGER
118: * The dimension of array IWORK.
119: * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
120: * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
121: *
122: * If LIWORK = -1, then a workspace query is assumed; the
123: * routine only calculates the optimal sizes of the WORK, RWORK
124: * and IWORK arrays, returns these values as the first entries
125: * of the WORK, RWORK and IWORK arrays, and no error message
126: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
127: *
128: * INFO (output) INTEGER
129: * = 0: successful exit.
130: * < 0: if INFO = -i, the i-th argument had an illegal value.
131: * > 0: if INFO = i, the algorithm failed to converge; i
132: * off-diagonal elements of an intermediate tridiagonal
133: * form did not converge to zero.
134: *
135: * =====================================================================
136: *
137: * .. Parameters ..
138: DOUBLE PRECISION ZERO, ONE
139: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
140: COMPLEX*16 CZERO, CONE
141: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
142: $ CONE = ( 1.0D0, 0.0D0 ) )
143: * ..
144: * .. Local Scalars ..
145: LOGICAL LOWER, LQUERY, WANTZ
146: INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
147: $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
148: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
149: $ SMLNUM
150: * ..
151: * .. External Functions ..
152: LOGICAL LSAME
153: DOUBLE PRECISION DLAMCH, ZLANHB
154: EXTERNAL LSAME, DLAMCH, ZLANHB
155: * ..
156: * .. External Subroutines ..
157: EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
158: $ ZLASCL, ZSTEDC
159: * ..
160: * .. Intrinsic Functions ..
161: INTRINSIC SQRT
162: * ..
163: * .. Executable Statements ..
164: *
165: * Test the input parameters.
166: *
167: WANTZ = LSAME( JOBZ, 'V' )
168: LOWER = LSAME( UPLO, 'L' )
169: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
170: *
171: INFO = 0
172: IF( N.LE.1 ) THEN
173: LWMIN = 1
174: LRWMIN = 1
175: LIWMIN = 1
176: ELSE
177: IF( WANTZ ) THEN
178: LWMIN = 2*N**2
179: LRWMIN = 1 + 5*N + 2*N**2
180: LIWMIN = 3 + 5*N
181: ELSE
182: LWMIN = N
183: LRWMIN = N
184: LIWMIN = 1
185: END IF
186: END IF
187: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
188: INFO = -1
189: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
190: INFO = -2
191: ELSE IF( N.LT.0 ) THEN
192: INFO = -3
193: ELSE IF( KD.LT.0 ) THEN
194: INFO = -4
195: ELSE IF( LDAB.LT.KD+1 ) THEN
196: INFO = -6
197: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
198: INFO = -9
199: END IF
200: *
201: IF( INFO.EQ.0 ) THEN
202: WORK( 1 ) = LWMIN
203: RWORK( 1 ) = LRWMIN
204: IWORK( 1 ) = LIWMIN
205: *
206: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
207: INFO = -11
208: ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
209: INFO = -13
210: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
211: INFO = -15
212: END IF
213: END IF
214: *
215: IF( INFO.NE.0 ) THEN
216: CALL XERBLA( 'ZHBEVD', -INFO )
217: RETURN
218: ELSE IF( LQUERY ) THEN
219: RETURN
220: END IF
221: *
222: * Quick return if possible
223: *
224: IF( N.EQ.0 )
225: $ RETURN
226: *
227: IF( N.EQ.1 ) THEN
228: W( 1 ) = AB( 1, 1 )
229: IF( WANTZ )
230: $ Z( 1, 1 ) = CONE
231: RETURN
232: END IF
233: *
234: * Get machine constants.
235: *
236: SAFMIN = DLAMCH( 'Safe minimum' )
237: EPS = DLAMCH( 'Precision' )
238: SMLNUM = SAFMIN / EPS
239: BIGNUM = ONE / SMLNUM
240: RMIN = SQRT( SMLNUM )
241: RMAX = SQRT( BIGNUM )
242: *
243: * Scale matrix to allowable range, if necessary.
244: *
245: ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
246: ISCALE = 0
247: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
248: ISCALE = 1
249: SIGMA = RMIN / ANRM
250: ELSE IF( ANRM.GT.RMAX ) THEN
251: ISCALE = 1
252: SIGMA = RMAX / ANRM
253: END IF
254: IF( ISCALE.EQ.1 ) THEN
255: IF( LOWER ) THEN
256: CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
257: ELSE
258: CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
259: END IF
260: END IF
261: *
262: * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
263: *
264: INDE = 1
265: INDWRK = INDE + N
266: INDWK2 = 1 + N*N
267: LLWK2 = LWORK - INDWK2 + 1
268: LLRWK = LRWORK - INDWRK + 1
269: CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
270: $ LDZ, WORK, IINFO )
271: *
272: * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
273: *
274: IF( .NOT.WANTZ ) THEN
275: CALL DSTERF( N, W, RWORK( INDE ), INFO )
276: ELSE
277: CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
278: $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
279: $ INFO )
280: CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
281: $ WORK( INDWK2 ), N )
282: CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
283: END IF
284: *
285: * If matrix was scaled, then rescale eigenvalues appropriately.
286: *
287: IF( ISCALE.EQ.1 ) THEN
288: IF( INFO.EQ.0 ) THEN
289: IMAX = N
290: ELSE
291: IMAX = INFO - 1
292: END IF
293: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
294: END IF
295: *
296: WORK( 1 ) = LWMIN
297: RWORK( 1 ) = LRWMIN
298: IWORK( 1 ) = LIWMIN
299: RETURN
300: *
301: * End of ZHBEVD
302: *
303: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhbevd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack