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