LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cheevd.f
Go to the documentation of this file.
1*> \brief <b> CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHEEVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
22* LRWORK, IWORK, LIWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBZ, UPLO
26* INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
27* ..
28* .. Array Arguments ..
29* INTEGER IWORK( * )
30* REAL RWORK( * ), W( * )
31* COMPLEX A( LDA, * ), WORK( * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> CHEEVD computes all eigenvalues and, optionally, eigenvectors of a
41*> complex Hermitian matrix A. If eigenvectors are desired, it uses a
42*> divide and conquer algorithm.
43*>
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] JOBZ
50*> \verbatim
51*> JOBZ is CHARACTER*1
52*> = 'N': Compute eigenvalues only;
53*> = 'V': Compute eigenvalues and eigenvectors.
54*> \endverbatim
55*>
56*> \param[in] UPLO
57*> \verbatim
58*> UPLO is CHARACTER*1
59*> = 'U': Upper triangle of A is stored;
60*> = 'L': Lower triangle of A is stored.
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrix A. N >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*> A is COMPLEX array, dimension (LDA, N)
72*> On entry, the Hermitian matrix A. If UPLO = 'U', the
73*> leading N-by-N upper triangular part of A contains the
74*> upper triangular part of the matrix A. If UPLO = 'L',
75*> the leading N-by-N lower triangular part of A contains
76*> the lower triangular part of the matrix A.
77*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
78*> orthonormal eigenvectors of the matrix A.
79*> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
80*> or the upper triangle (if UPLO='U') of A, including the
81*> diagonal, is destroyed.
82*> \endverbatim
83*>
84*> \param[in] LDA
85*> \verbatim
86*> LDA is INTEGER
87*> The leading dimension of the array A. LDA >= max(1,N).
88*> \endverbatim
89*>
90*> \param[out] W
91*> \verbatim
92*> W is REAL array, dimension (N)
93*> If INFO = 0, the eigenvalues in ascending order.
94*> \endverbatim
95*>
96*> \param[out] WORK
97*> \verbatim
98*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
99*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
100*> \endverbatim
101*>
102*> \param[in] LWORK
103*> \verbatim
104*> LWORK is INTEGER
105*> The length of the array WORK.
106*> If N <= 1, LWORK must be at least 1.
107*> If JOBZ = 'N' and N > 1, LWORK must be at least N + 1.
108*> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + 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*> \endverbatim
116*>
117*> \param[out] RWORK
118*> \verbatim
119*> RWORK is REAL array,
120*> dimension (LRWORK)
121*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
122*> \endverbatim
123*>
124*> \param[in] LRWORK
125*> \verbatim
126*> LRWORK is INTEGER
127*> The dimension of the array RWORK.
128*> If N <= 1, LRWORK must be at least 1.
129*> If JOBZ = 'N' and N > 1, LRWORK must be at least N.
130*> If JOBZ = 'V' and N > 1, LRWORK must be at least
131*> 1 + 5*N + 2*N**2.
132*>
133*> If LRWORK = -1, then a workspace query is assumed; the
134*> routine only calculates the optimal sizes of the WORK, RWORK
135*> and IWORK arrays, returns these values as the first entries
136*> of the WORK, RWORK and IWORK arrays, and no error message
137*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
138*> \endverbatim
139*>
140*> \param[out] IWORK
141*> \verbatim
142*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
143*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
144*> \endverbatim
145*>
146*> \param[in] LIWORK
147*> \verbatim
148*> LIWORK is INTEGER
149*> The dimension of the array IWORK.
150*> If N <= 1, LIWORK must be at least 1.
151*> If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
152*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
153*>
154*> If LIWORK = -1, then a workspace query is assumed; the
155*> routine only calculates the optimal sizes of the WORK, RWORK
156*> and IWORK arrays, returns these values as the first entries
157*> of the WORK, RWORK and IWORK arrays, and no error message
158*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
159*> \endverbatim
160*>
161*> \param[out] INFO
162*> \verbatim
163*> INFO is INTEGER
164*> = 0: successful exit
165*> < 0: if INFO = -i, the i-th argument had an illegal value
166*> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
167*> to converge; i off-diagonal elements of an intermediate
168*> tridiagonal form did not converge to zero;
169*> if INFO = i and JOBZ = 'V', then the algorithm failed
170*> to compute an eigenvalue while working on the submatrix
171*> lying in rows and columns INFO/(N+1) through
172*> mod(INFO,N+1).
173*> \endverbatim
174*
175* Authors:
176* ========
177*
178*> \author Univ. of Tennessee
179*> \author Univ. of California Berkeley
180*> \author Univ. of Colorado Denver
181*> \author NAG Ltd.
182*
183*> \ingroup heevd
184*
185*> \par Further Details:
186* =====================
187*>
188*> Modified description of INFO. Sven, 16 Feb 05.
189*
190*> \par Contributors:
191* ==================
192*>
193*> Jeff Rutter, Computer Science Division, University of California
194*> at Berkeley, USA
195*>
196* =====================================================================
197 SUBROUTINE cheevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
198 $ LRWORK, IWORK, LIWORK, INFO )
199*
200* -- LAPACK driver routine --
201* -- LAPACK is a software package provided by Univ. of Tennessee, --
202* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203*
204* .. Scalar Arguments ..
205 CHARACTER JOBZ, UPLO
206 INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
207* ..
208* .. Array Arguments ..
209 INTEGER IWORK( * )
210 REAL RWORK( * ), W( * )
211 COMPLEX A( LDA, * ), WORK( * )
212* ..
213*
214* =====================================================================
215*
216* .. Parameters ..
217 REAL ZERO, ONE
218 parameter( zero = 0.0e0, one = 1.0e0 )
219 COMPLEX CONE
220 parameter( cone = ( 1.0e0, 0.0e0 ) )
221* ..
222* .. Local Scalars ..
223 LOGICAL LOWER, LQUERY, WANTZ
224 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
225 $ indwrk, iscale, liopt, liwmin, llrwk, llwork,
226 $ llwrk2, lopt, lropt, lrwmin, lwmin
227 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
228 $ smlnum
229* ..
230* .. External Functions ..
231 LOGICAL LSAME
232 INTEGER ILAENV
233 REAL CLANHE, SLAMCH, SROUNDUP_LWORK
234 EXTERNAL ilaenv, lsame, clanhe, slamch, sroundup_lwork
235* ..
236* .. External Subroutines ..
237 EXTERNAL chetrd, clacpy, clascl, cstedc, cunmtr, sscal,
238 $ ssterf, xerbla
239* ..
240* .. Intrinsic Functions ..
241 INTRINSIC max, sqrt
242* ..
243* .. Executable Statements ..
244*
245* Test the input parameters.
246*
247 wantz = lsame( jobz, 'V' )
248 lower = lsame( uplo, 'L' )
249 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
250*
251 info = 0
252 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
253 info = -1
254 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
255 info = -2
256 ELSE IF( n.LT.0 ) THEN
257 info = -3
258 ELSE IF( lda.LT.max( 1, n ) ) THEN
259 info = -5
260 END IF
261*
262 IF( info.EQ.0 ) THEN
263 IF( n.LE.1 ) THEN
264 lwmin = 1
265 lrwmin = 1
266 liwmin = 1
267 lopt = lwmin
268 lropt = lrwmin
269 liopt = liwmin
270 ELSE
271 IF( wantz ) THEN
272 lwmin = 2*n + n*n
273 lrwmin = 1 + 5*n + 2*n**2
274 liwmin = 3 + 5*n
275 ELSE
276 lwmin = n + 1
277 lrwmin = n
278 liwmin = 1
279 END IF
280 lopt = max( lwmin, n +
281 $ n*ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 ) )
282 lropt = lrwmin
283 liopt = liwmin
284 END IF
285 work( 1 ) = sroundup_lwork(lopt)
286 rwork( 1 ) = lropt
287 iwork( 1 ) = liopt
288*
289 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
290 info = -8
291 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
292 info = -10
293 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
294 info = -12
295 END IF
296 END IF
297*
298 IF( info.NE.0 ) THEN
299 CALL xerbla( 'CHEEVD', -info )
300 RETURN
301 ELSE IF( lquery ) THEN
302 RETURN
303 END IF
304*
305* Quick return if possible
306*
307 IF( n.EQ.0 )
308 $ RETURN
309*
310 IF( n.EQ.1 ) THEN
311 w( 1 ) = real( a( 1, 1 ) )
312 IF( wantz )
313 $ a( 1, 1 ) = cone
314 RETURN
315 END IF
316*
317* Get machine constants.
318*
319 safmin = slamch( 'Safe minimum' )
320 eps = slamch( 'Precision' )
321 smlnum = safmin / eps
322 bignum = one / smlnum
323 rmin = sqrt( smlnum )
324 rmax = sqrt( bignum )
325*
326* Scale matrix to allowable range, if necessary.
327*
328 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
329 iscale = 0
330 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
331 iscale = 1
332 sigma = rmin / anrm
333 ELSE IF( anrm.GT.rmax ) THEN
334 iscale = 1
335 sigma = rmax / anrm
336 END IF
337 IF( iscale.EQ.1 )
338 $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
339*
340* Call CHETRD to reduce Hermitian matrix to tridiagonal form.
341*
342 inde = 1
343 indtau = 1
344 indwrk = indtau + n
345 indrwk = inde + n
346 indwk2 = indwrk + n*n
347 llwork = lwork - indwrk + 1
348 llwrk2 = lwork - indwk2 + 1
349 llrwk = lrwork - indrwk + 1
350 CALL chetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
351 $ work( indwrk ), llwork, iinfo )
352*
353* For eigenvalues only, call SSTERF. For eigenvectors, first call
354* CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
355* tridiagonal matrix, then call CUNMTR to multiply it to the
356* Householder transformations represented as Householder vectors in
357* A.
358*
359 IF( .NOT.wantz ) THEN
360 CALL ssterf( n, w, rwork( inde ), info )
361 ELSE
362 CALL cstedc( 'I', n, w, rwork( inde ), work( indwrk ), n,
363 $ work( indwk2 ), llwrk2, rwork( indrwk ), llrwk,
364 $ iwork, liwork, info )
365 CALL cunmtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
366 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
367 CALL clacpy( 'A', n, n, work( indwrk ), n, a, lda )
368 END IF
369*
370* If matrix was scaled, then rescale eigenvalues appropriately.
371*
372 IF( iscale.EQ.1 ) THEN
373 IF( info.EQ.0 ) THEN
374 imax = n
375 ELSE
376 imax = info - 1
377 END IF
378 CALL sscal( imax, one / sigma, w, 1 )
379 END IF
380*
381 work( 1 ) = sroundup_lwork(lopt)
382 rwork( 1 ) = lropt
383 iwork( 1 ) = liopt
384*
385 RETURN
386*
387* End of CHEEVD
388*
389 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cheevd(jobz, uplo, n, a, lda, w, work, lwork, rwork, lrwork, iwork, liwork, info)
CHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition cheevd.f:199
subroutine chetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
CHETRD
Definition chetrd.f:192
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine cstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
CSTEDC
Definition cstedc.f:206
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine cunmtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
CUNMTR
Definition cunmtr.f:172