LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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  $ 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