LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssyevd.f
Go to the documentation of this file.
1*> \brief <b> SSYEVD 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
9*> Download SSYEVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyevd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyevd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyevd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSYEVD( 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* REAL A( LDA, * ), W( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SSYEVD 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, SSYEVD needs N**2 more
44*> workspace than SSYEVX.
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 REAL 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 REAL array, dimension (N)
94*> If INFO = 0, the eigenvalues in ascending order.
95*> \endverbatim
96*>
97*> \param[out] WORK
98*> \verbatim
99*> WORK is REAL 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 SUBROUTINE ssyevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
175 $ LIWORK, INFO )
176*
177* -- LAPACK driver routine --
178* -- LAPACK is a software package provided by Univ. of Tennessee, --
179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181* .. Scalar Arguments ..
182 CHARACTER JOBZ, UPLO
183 INTEGER INFO, LDA, LIWORK, LWORK, N
184* ..
185* .. Array Arguments ..
186 INTEGER IWORK( * )
187 REAL A( LDA, * ), W( * ), WORK( * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 REAL ZERO, ONE
194 parameter( zero = 0.0e+0, one = 1.0e+0 )
195* ..
196* .. Local Scalars ..
197*
198 LOGICAL LOWER, LQUERY, WANTZ
199 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
200 $ liopt, liwmin, llwork, llwrk2, lopt, lwmin
201 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
202 $ smlnum
203* ..
204* .. External Functions ..
205 LOGICAL LSAME
206 INTEGER ILAENV
207 REAL SLAMCH, SLANSY, SROUNDUP_LWORK
208 EXTERNAL ilaenv, lsame, slamch, slansy, sroundup_lwork
209* ..
210* .. External Subroutines ..
211 EXTERNAL slacpy, slascl, sormtr, sscal, sstedc, ssterf,
212 $ ssytrd, xerbla
213* ..
214* .. Intrinsic Functions ..
215 INTRINSIC max, sqrt
216* ..
217* .. Executable Statements ..
218*
219* Test the input parameters.
220*
221 wantz = lsame( jobz, 'V' )
222 lower = lsame( uplo, 'L' )
223 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
224*
225 info = 0
226 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
227 info = -1
228 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
229 info = -2
230 ELSE IF( n.LT.0 ) THEN
231 info = -3
232 ELSE IF( lda.LT.max( 1, n ) ) THEN
233 info = -5
234 END IF
235*
236 IF( info.EQ.0 ) THEN
237 IF( n.LE.1 ) THEN
238 liwmin = 1
239 lwmin = 1
240 lopt = lwmin
241 liopt = liwmin
242 ELSE
243 IF( wantz ) THEN
244 liwmin = 3 + 5*n
245 lwmin = 1 + 6*n + 2*n**2
246 ELSE
247 liwmin = 1
248 lwmin = 2*n + 1
249 END IF
250 lopt = max( lwmin, 2*n +
251 $ n*ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 ) )
252 liopt = liwmin
253 END IF
254 work( 1 ) = sroundup_lwork(lopt)
255 iwork( 1 ) = liopt
256*
257 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
258 info = -8
259 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
260 info = -10
261 END IF
262 END IF
263*
264 IF( info.NE.0 ) THEN
265 CALL xerbla( 'SSYEVD', -info )
266 RETURN
267 ELSE IF( lquery ) THEN
268 RETURN
269 END IF
270*
271* Quick return if possible
272*
273 IF( n.EQ.0 )
274 $ RETURN
275*
276 IF( n.EQ.1 ) THEN
277 w( 1 ) = a( 1, 1 )
278 IF( wantz )
279 $ a( 1, 1 ) = one
280 RETURN
281 END IF
282*
283* Get machine constants.
284*
285 safmin = slamch( 'Safe minimum' )
286 eps = slamch( 'Precision' )
287 smlnum = safmin / eps
288 bignum = one / smlnum
289 rmin = sqrt( smlnum )
290 rmax = sqrt( bignum )
291*
292* Scale matrix to allowable range, if necessary.
293*
294 anrm = slansy( 'M', uplo, n, a, lda, work )
295 iscale = 0
296 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
297 iscale = 1
298 sigma = rmin / anrm
299 ELSE IF( anrm.GT.rmax ) THEN
300 iscale = 1
301 sigma = rmax / anrm
302 END IF
303 IF( iscale.EQ.1 )
304 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
305*
306* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
307*
308 inde = 1
309 indtau = inde + n
310 indwrk = indtau + n
311 llwork = lwork - indwrk + 1
312 indwk2 = indwrk + n*n
313 llwrk2 = lwork - indwk2 + 1
314*
315 CALL ssytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
316 $ work( indwrk ), llwork, iinfo )
317*
318* For eigenvalues only, call SSTERF. For eigenvectors, first call
319* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
320* tridiagonal matrix, then call SORMTR to multiply it by the
321* Householder transformations stored in A.
322*
323 IF( .NOT.wantz ) THEN
324 CALL ssterf( n, w, work( inde ), info )
325 ELSE
326 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
327 $ work( indwk2 ), llwrk2, iwork, liwork, info )
328 CALL sormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
329 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
330 CALL slacpy( 'A', n, n, work( indwrk ), n, a, lda )
331 END IF
332*
333* If matrix was scaled, then rescale eigenvalues appropriately.
334*
335 IF( iscale.EQ.1 )
336 $ CALL sscal( n, one / sigma, w, 1 )
337*
338 work( 1 ) = sroundup_lwork(lopt)
339 iwork( 1 ) = liopt
340*
341 RETURN
342*
343* End of SSYEVD
344*
345 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ssyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyevd.f:176
subroutine ssytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
SSYTRD
Definition ssytrd.f:192
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:143
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:182
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine sormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
SORMTR
Definition sormtr.f:172