1 SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
2 $ 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, LDA, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 A( LDA, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
22 * complex Hermitian matrix A.
23 *
24 * Arguments
25 * =========
26 *
27 * JOBZ (input) CHARACTER*1
28 * = 'N': Compute eigenvalues only;
29 * = 'V': Compute eigenvalues and eigenvectors.
30 *
31 * UPLO (input) CHARACTER*1
32 * = 'U': Upper triangle of A is stored;
33 * = 'L': Lower triangle of A is stored.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
39 * On entry, the Hermitian matrix A. If UPLO = 'U', the
40 * leading N-by-N upper triangular part of A contains the
41 * upper triangular part of the matrix A. If UPLO = 'L',
42 * the leading N-by-N lower triangular part of A contains
43 * the lower triangular part of the matrix A.
44 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
45 * orthonormal eigenvectors of the matrix A.
46 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
47 * or the upper triangle (if UPLO='U') of A, including the
48 * diagonal, is destroyed.
49 *
50 * LDA (input) INTEGER
51 * The leading dimension of the array A. LDA >= max(1,N).
52 *
53 * W (output) DOUBLE PRECISION array, dimension (N)
54 * If INFO = 0, the eigenvalues in ascending order.
55 *
56 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
57 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
58 *
59 * LWORK (input) INTEGER
60 * The length of the array WORK. LWORK >= max(1,2*N-1).
61 * For optimal efficiency, LWORK >= (NB+1)*N,
62 * where NB is the blocksize for ZHETRD returned by ILAENV.
63 *
64 * If LWORK = -1, then a workspace query is assumed; the routine
65 * only calculates the optimal size of the WORK array, returns
66 * this value as the first entry of the WORK array, and no error
67 * message related to LWORK is issued by XERBLA.
68 *
69 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
70 *
71 * INFO (output) INTEGER
72 * = 0: successful exit
73 * < 0: if INFO = -i, the i-th argument had an illegal value
74 * > 0: if INFO = i, the algorithm failed to converge; i
75 * off-diagonal elements of an intermediate tridiagonal
76 * form did not converge to zero.
77 *
78 * =====================================================================
79 *
80 * .. Parameters ..
81 DOUBLE PRECISION ZERO, ONE
82 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
83 COMPLEX*16 CONE
84 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
85 * ..
86 * .. Local Scalars ..
87 LOGICAL LOWER, LQUERY, WANTZ
88 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
89 $ LLWORK, LWKOPT, NB
90 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
91 $ SMLNUM
92 * ..
93 * .. External Functions ..
94 LOGICAL LSAME
95 INTEGER ILAENV
96 DOUBLE PRECISION DLAMCH, ZLANHE
97 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE
98 * ..
99 * .. External Subroutines ..
100 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR,
101 $ ZUNGTR
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC MAX, SQRT
105 * ..
106 * .. Executable Statements ..
107 *
108 * Test the input parameters.
109 *
110 WANTZ = LSAME( JOBZ, 'V' )
111 LOWER = LSAME( UPLO, 'L' )
112 LQUERY = ( LWORK.EQ.-1 )
113 *
114 INFO = 0
115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
116 INFO = -1
117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
118 INFO = -2
119 ELSE IF( N.LT.0 ) THEN
120 INFO = -3
121 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
122 INFO = -5
123 END IF
124 *
125 IF( INFO.EQ.0 ) THEN
126 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
127 LWKOPT = MAX( 1, ( NB+1 )*N )
128 WORK( 1 ) = LWKOPT
129 *
130 IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY )
131 $ INFO = -8
132 END IF
133 *
134 IF( INFO.NE.0 ) THEN
135 CALL XERBLA( 'ZHEEV ', -INFO )
136 RETURN
137 ELSE IF( LQUERY ) THEN
138 RETURN
139 END IF
140 *
141 * Quick return if possible
142 *
143 IF( N.EQ.0 ) THEN
144 RETURN
145 END IF
146 *
147 IF( N.EQ.1 ) THEN
148 W( 1 ) = A( 1, 1 )
149 WORK( 1 ) = 1
150 IF( WANTZ )
151 $ A( 1, 1 ) = CONE
152 RETURN
153 END IF
154 *
155 * Get machine constants.
156 *
157 SAFMIN = DLAMCH( 'Safe minimum' )
158 EPS = DLAMCH( 'Precision' )
159 SMLNUM = SAFMIN / EPS
160 BIGNUM = ONE / SMLNUM
161 RMIN = SQRT( SMLNUM )
162 RMAX = SQRT( BIGNUM )
163 *
164 * Scale matrix to allowable range, if necessary.
165 *
166 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
167 ISCALE = 0
168 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
169 ISCALE = 1
170 SIGMA = RMIN / ANRM
171 ELSE IF( ANRM.GT.RMAX ) THEN
172 ISCALE = 1
173 SIGMA = RMAX / ANRM
174 END IF
175 IF( ISCALE.EQ.1 )
176 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
177 *
178 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
179 *
180 INDE = 1
181 INDTAU = 1
182 INDWRK = INDTAU + N
183 LLWORK = LWORK - INDWRK + 1
184 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
185 $ WORK( INDWRK ), LLWORK, IINFO )
186 *
187 * For eigenvalues only, call DSTERF. For eigenvectors, first call
188 * ZUNGTR to generate the unitary matrix, then call ZSTEQR.
189 *
190 IF( .NOT.WANTZ ) THEN
191 CALL DSTERF( N, W, RWORK( INDE ), INFO )
192 ELSE
193 CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
194 $ LLWORK, IINFO )
195 INDWRK = INDE + N
196 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA,
197 $ RWORK( INDWRK ), INFO )
198 END IF
199 *
200 * If matrix was scaled, then rescale eigenvalues appropriately.
201 *
202 IF( ISCALE.EQ.1 ) THEN
203 IF( INFO.EQ.0 ) THEN
204 IMAX = N
205 ELSE
206 IMAX = INFO - 1
207 END IF
208 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
209 END IF
210 *
211 * Set WORK(1) to optimal complex workspace size.
212 *
213 WORK( 1 ) = LWKOPT
214 *
215 RETURN
216 *
217 * End of ZHEEV
218 *
219 END
2 $ 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, LDA, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * ), W( * )
15 COMPLEX*16 A( LDA, * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
22 * complex Hermitian matrix A.
23 *
24 * Arguments
25 * =========
26 *
27 * JOBZ (input) CHARACTER*1
28 * = 'N': Compute eigenvalues only;
29 * = 'V': Compute eigenvalues and eigenvectors.
30 *
31 * UPLO (input) CHARACTER*1
32 * = 'U': Upper triangle of A is stored;
33 * = 'L': Lower triangle of A is stored.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
39 * On entry, the Hermitian matrix A. If UPLO = 'U', the
40 * leading N-by-N upper triangular part of A contains the
41 * upper triangular part of the matrix A. If UPLO = 'L',
42 * the leading N-by-N lower triangular part of A contains
43 * the lower triangular part of the matrix A.
44 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
45 * orthonormal eigenvectors of the matrix A.
46 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
47 * or the upper triangle (if UPLO='U') of A, including the
48 * diagonal, is destroyed.
49 *
50 * LDA (input) INTEGER
51 * The leading dimension of the array A. LDA >= max(1,N).
52 *
53 * W (output) DOUBLE PRECISION array, dimension (N)
54 * If INFO = 0, the eigenvalues in ascending order.
55 *
56 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
57 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
58 *
59 * LWORK (input) INTEGER
60 * The length of the array WORK. LWORK >= max(1,2*N-1).
61 * For optimal efficiency, LWORK >= (NB+1)*N,
62 * where NB is the blocksize for ZHETRD returned by ILAENV.
63 *
64 * If LWORK = -1, then a workspace query is assumed; the routine
65 * only calculates the optimal size of the WORK array, returns
66 * this value as the first entry of the WORK array, and no error
67 * message related to LWORK is issued by XERBLA.
68 *
69 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
70 *
71 * INFO (output) INTEGER
72 * = 0: successful exit
73 * < 0: if INFO = -i, the i-th argument had an illegal value
74 * > 0: if INFO = i, the algorithm failed to converge; i
75 * off-diagonal elements of an intermediate tridiagonal
76 * form did not converge to zero.
77 *
78 * =====================================================================
79 *
80 * .. Parameters ..
81 DOUBLE PRECISION ZERO, ONE
82 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
83 COMPLEX*16 CONE
84 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
85 * ..
86 * .. Local Scalars ..
87 LOGICAL LOWER, LQUERY, WANTZ
88 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
89 $ LLWORK, LWKOPT, NB
90 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
91 $ SMLNUM
92 * ..
93 * .. External Functions ..
94 LOGICAL LSAME
95 INTEGER ILAENV
96 DOUBLE PRECISION DLAMCH, ZLANHE
97 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE
98 * ..
99 * .. External Subroutines ..
100 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR,
101 $ ZUNGTR
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC MAX, SQRT
105 * ..
106 * .. Executable Statements ..
107 *
108 * Test the input parameters.
109 *
110 WANTZ = LSAME( JOBZ, 'V' )
111 LOWER = LSAME( UPLO, 'L' )
112 LQUERY = ( LWORK.EQ.-1 )
113 *
114 INFO = 0
115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
116 INFO = -1
117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
118 INFO = -2
119 ELSE IF( N.LT.0 ) THEN
120 INFO = -3
121 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
122 INFO = -5
123 END IF
124 *
125 IF( INFO.EQ.0 ) THEN
126 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
127 LWKOPT = MAX( 1, ( NB+1 )*N )
128 WORK( 1 ) = LWKOPT
129 *
130 IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY )
131 $ INFO = -8
132 END IF
133 *
134 IF( INFO.NE.0 ) THEN
135 CALL XERBLA( 'ZHEEV ', -INFO )
136 RETURN
137 ELSE IF( LQUERY ) THEN
138 RETURN
139 END IF
140 *
141 * Quick return if possible
142 *
143 IF( N.EQ.0 ) THEN
144 RETURN
145 END IF
146 *
147 IF( N.EQ.1 ) THEN
148 W( 1 ) = A( 1, 1 )
149 WORK( 1 ) = 1
150 IF( WANTZ )
151 $ A( 1, 1 ) = CONE
152 RETURN
153 END IF
154 *
155 * Get machine constants.
156 *
157 SAFMIN = DLAMCH( 'Safe minimum' )
158 EPS = DLAMCH( 'Precision' )
159 SMLNUM = SAFMIN / EPS
160 BIGNUM = ONE / SMLNUM
161 RMIN = SQRT( SMLNUM )
162 RMAX = SQRT( BIGNUM )
163 *
164 * Scale matrix to allowable range, if necessary.
165 *
166 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
167 ISCALE = 0
168 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
169 ISCALE = 1
170 SIGMA = RMIN / ANRM
171 ELSE IF( ANRM.GT.RMAX ) THEN
172 ISCALE = 1
173 SIGMA = RMAX / ANRM
174 END IF
175 IF( ISCALE.EQ.1 )
176 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
177 *
178 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
179 *
180 INDE = 1
181 INDTAU = 1
182 INDWRK = INDTAU + N
183 LLWORK = LWORK - INDWRK + 1
184 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
185 $ WORK( INDWRK ), LLWORK, IINFO )
186 *
187 * For eigenvalues only, call DSTERF. For eigenvectors, first call
188 * ZUNGTR to generate the unitary matrix, then call ZSTEQR.
189 *
190 IF( .NOT.WANTZ ) THEN
191 CALL DSTERF( N, W, RWORK( INDE ), INFO )
192 ELSE
193 CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
194 $ LLWORK, IINFO )
195 INDWRK = INDE + N
196 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA,
197 $ RWORK( INDWRK ), INFO )
198 END IF
199 *
200 * If matrix was scaled, then rescale eigenvalues appropriately.
201 *
202 IF( ISCALE.EQ.1 ) THEN
203 IF( INFO.EQ.0 ) THEN
204 IMAX = N
205 ELSE
206 IMAX = INFO - 1
207 END IF
208 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
209 END IF
210 *
211 * Set WORK(1) to optimal complex workspace size.
212 *
213 WORK( 1 ) = LWKOPT
214 *
215 RETURN
216 *
217 * End of ZHEEV
218 *
219 END