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