LAPACK 3.11.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*> The divide and conquer algorithm makes very mild assumptions about
44*> floating point arithmetic. It will work on machines with a guard
45*> digit in add/subtract, or on those binary machines without guard
46*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
47*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
48*> without guard digits, but we know of none.
49*>
50*> Because of large use of BLAS of level 3, SSYEVD needs N**2 more
51*> workspace than SSYEVX.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] JOBZ
58*> \verbatim
59*> JOBZ is CHARACTER*1
60*> = 'N': Compute eigenvalues only;
61*> = 'V': Compute eigenvalues and eigenvectors.
62*> \endverbatim
63*>
64*> \param[in] UPLO
65*> \verbatim
66*> UPLO is CHARACTER*1
67*> = 'U': Upper triangle of A is stored;
68*> = 'L': Lower triangle of A is stored.
69*> \endverbatim
70*>
71*> \param[in] N
72*> \verbatim
73*> N is INTEGER
74*> The order of the matrix A. N >= 0.
75*> \endverbatim
76*>
77*> \param[in,out] A
78*> \verbatim
79*> A is REAL array, dimension (LDA, N)
80*> On entry, the symmetric matrix A. If UPLO = 'U', the
81*> leading N-by-N upper triangular part of A contains the
82*> upper triangular part of the matrix A. If UPLO = 'L',
83*> the leading N-by-N lower triangular part of A contains
84*> the lower triangular part of the matrix A.
85*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
86*> orthonormal eigenvectors of the matrix A.
87*> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
88*> or the upper triangle (if UPLO='U') of A, including the
89*> diagonal, is destroyed.
90*> \endverbatim
91*>
92*> \param[in] LDA
93*> \verbatim
94*> LDA is INTEGER
95*> The leading dimension of the array A. LDA >= max(1,N).
96*> \endverbatim
97*>
98*> \param[out] W
99*> \verbatim
100*> W is REAL array, dimension (N)
101*> If INFO = 0, the eigenvalues in ascending order.
102*> \endverbatim
103*>
104*> \param[out] WORK
105*> \verbatim
106*> WORK is REAL array,
107*> dimension (LWORK)
108*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
109*> \endverbatim
110*>
111*> \param[in] LWORK
112*> \verbatim
113*> LWORK is INTEGER
114*> The dimension of the array WORK.
115*> If N <= 1, LWORK must be at least 1.
116*> If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
117*> If JOBZ = 'V' and N > 1, LWORK must be at least
118*> 1 + 6*N + 2*N**2.
119*>
120*> If LWORK = -1, then a workspace query is assumed; the routine
121*> only calculates the optimal sizes of the WORK and IWORK
122*> arrays, returns these values as the first entries of the WORK
123*> and IWORK arrays, and no error message related to LWORK or
124*> LIWORK is issued by XERBLA.
125*> \endverbatim
126*>
127*> \param[out] IWORK
128*> \verbatim
129*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
130*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
131*> \endverbatim
132*>
133*> \param[in] LIWORK
134*> \verbatim
135*> LIWORK is INTEGER
136*> The dimension of the array IWORK.
137*> If N <= 1, LIWORK must be at least 1.
138*> If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
139*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
140*>
141*> If LIWORK = -1, then a workspace query is assumed; the
142*> routine only calculates the optimal sizes of the WORK and
143*> IWORK arrays, returns these values as the first entries of
144*> the WORK and IWORK arrays, and no error message related to
145*> LWORK or LIWORK is issued by XERBLA.
146*> \endverbatim
147*>
148*> \param[out] INFO
149*> \verbatim
150*> INFO is INTEGER
151*> = 0: successful exit
152*> < 0: if INFO = -i, the i-th argument had an illegal value
153*> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
154*> to converge; i off-diagonal elements of an intermediate
155*> tridiagonal form did not converge to zero;
156*> if INFO = i and JOBZ = 'V', then the algorithm failed
157*> to compute an eigenvalue while working on the submatrix
158*> lying in rows and columns INFO/(N+1) through
159*> mod(INFO,N+1).
160*> \endverbatim
161*
162* Authors:
163* ========
164*
165*> \author Univ. of Tennessee
166*> \author Univ. of California Berkeley
167*> \author Univ. of Colorado Denver
168*> \author NAG Ltd.
169*
170*> \ingroup realSYeigen
171*
172*> \par Contributors:
173* ==================
174*>
175*> Jeff Rutter, Computer Science Division, University of California
176*> at Berkeley, USA \n
177*> Modified by Francoise Tisseur, University of Tennessee \n
178*> Modified description of INFO. Sven, 16 Feb 05. \n
179*>
180* =====================================================================
181 SUBROUTINE ssyevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
182 $ LIWORK, INFO )
183*
184* -- LAPACK driver routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER JOBZ, UPLO
190 INTEGER INFO, LDA, LIWORK, LWORK, N
191* ..
192* .. Array Arguments ..
193 INTEGER IWORK( * )
194 REAL A( LDA, * ), W( * ), WORK( * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 REAL ZERO, ONE
201 parameter( zero = 0.0e+0, one = 1.0e+0 )
202* ..
203* .. Local Scalars ..
204*
205 LOGICAL LOWER, LQUERY, WANTZ
206 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
207 $ liopt, liwmin, llwork, llwrk2, lopt, lwmin
208 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
209 $ smlnum
210* ..
211* .. External Functions ..
212 LOGICAL LSAME
213 INTEGER ILAENV
214 REAL SLAMCH, SLANSY
215 EXTERNAL ilaenv, lsame, slamch, slansy
216* ..
217* .. External Subroutines ..
218 EXTERNAL slacpy, slascl, sormtr, sscal, sstedc, ssterf,
219 $ ssytrd, xerbla
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC max, sqrt
223* ..
224* .. Executable Statements ..
225*
226* Test the input parameters.
227*
228 wantz = lsame( jobz, 'V' )
229 lower = lsame( uplo, 'L' )
230 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
231*
232 info = 0
233 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
234 info = -1
235 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
236 info = -2
237 ELSE IF( n.LT.0 ) THEN
238 info = -3
239 ELSE IF( lda.LT.max( 1, n ) ) THEN
240 info = -5
241 END IF
242*
243 IF( info.EQ.0 ) THEN
244 IF( n.LE.1 ) THEN
245 liwmin = 1
246 lwmin = 1
247 lopt = lwmin
248 liopt = liwmin
249 ELSE
250 IF( wantz ) THEN
251 liwmin = 3 + 5*n
252 lwmin = 1 + 6*n + 2*n**2
253 ELSE
254 liwmin = 1
255 lwmin = 2*n + 1
256 END IF
257 lopt = max( lwmin, 2*n +
258 $ n*ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 ) )
259 liopt = liwmin
260 END IF
261 work( 1 ) = lopt
262 iwork( 1 ) = liopt
263*
264 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
265 info = -8
266 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
267 info = -10
268 END IF
269 END IF
270*
271 IF( info.NE.0 ) THEN
272 CALL xerbla( 'SSYEVD', -info )
273 RETURN
274 ELSE IF( lquery ) THEN
275 RETURN
276 END IF
277*
278* Quick return if possible
279*
280 IF( n.EQ.0 )
281 $ RETURN
282*
283 IF( n.EQ.1 ) THEN
284 w( 1 ) = a( 1, 1 )
285 IF( wantz )
286 $ a( 1, 1 ) = one
287 RETURN
288 END IF
289*
290* Get machine constants.
291*
292 safmin = slamch( 'Safe minimum' )
293 eps = slamch( 'Precision' )
294 smlnum = safmin / eps
295 bignum = one / smlnum
296 rmin = sqrt( smlnum )
297 rmax = sqrt( bignum )
298*
299* Scale matrix to allowable range, if necessary.
300*
301 anrm = slansy( 'M', uplo, n, a, lda, work )
302 iscale = 0
303 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
304 iscale = 1
305 sigma = rmin / anrm
306 ELSE IF( anrm.GT.rmax ) THEN
307 iscale = 1
308 sigma = rmax / anrm
309 END IF
310 IF( iscale.EQ.1 )
311 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
312*
313* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
314*
315 inde = 1
316 indtau = inde + n
317 indwrk = indtau + n
318 llwork = lwork - indwrk + 1
319 indwk2 = indwrk + n*n
320 llwrk2 = lwork - indwk2 + 1
321*
322 CALL ssytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
323 $ work( indwrk ), llwork, iinfo )
324*
325* For eigenvalues only, call SSTERF. For eigenvectors, first call
326* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
327* tridiagonal matrix, then call SORMTR to multiply it by the
328* Householder transformations stored in A.
329*
330 IF( .NOT.wantz ) THEN
331 CALL ssterf( n, w, work( inde ), info )
332 ELSE
333 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
334 $ work( indwk2 ), llwrk2, iwork, liwork, info )
335 CALL sormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
336 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
337 CALL slacpy( 'A', n, n, work( indwrk ), n, a, lda )
338 END IF
339*
340* If matrix was scaled, then rescale eigenvalues appropriately.
341*
342 IF( iscale.EQ.1 )
343 $ CALL sscal( n, one / sigma, w, 1 )
344*
345 work( 1 ) = lopt
346 iwork( 1 ) = liopt
347*
348 RETURN
349*
350* End of SSYEVD
351*
352 END
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 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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEDC
Definition: sstedc.f:188
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
subroutine ssytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
SSYTRD
Definition: ssytrd.f:192
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:183
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79