LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> \date November 2011
171 *
172 *> \ingroup realSYeigen
173 *
174 *> \par Contributors:
175 * ==================
176 *>
177 *> Jeff Rutter, Computer Science Division, University of California
178 *> at Berkeley, USA \n
179 *> Modified by Francoise Tisseur, University of Tennessee \n
180 *> Modified description of INFO. Sven, 16 Feb 05. \n
181 *>
182 * =====================================================================
183  SUBROUTINE ssyevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
184  $ liwork, info )
185 *
186 * -- LAPACK driver routine (version 3.4.0) --
187 * -- LAPACK is a software package provided by Univ. of Tennessee, --
188 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189 * November 2011
190 *
191 * .. Scalar Arguments ..
192  CHARACTER jobz, uplo
193  INTEGER info, lda, liwork, lwork, n
194 * ..
195 * .. Array Arguments ..
196  INTEGER iwork( * )
197  REAL a( lda, * ), w( * ), work( * )
198 * ..
199 *
200 * =====================================================================
201 *
202 * .. Parameters ..
203  REAL zero, one
204  parameter( zero = 0.0e+0, one = 1.0e+0 )
205 * ..
206 * .. Local Scalars ..
207 *
208  LOGICAL lower, lquery, wantz
209  INTEGER iinfo, inde, indtau, indwk2, indwrk, iscale,
210  $ liopt, liwmin, llwork, llwrk2, lopt, lwmin
211  REAL anrm, bignum, eps, rmax, rmin, safmin, sigma,
212  $ smlnum
213 * ..
214 * .. External Functions ..
215  LOGICAL lsame
216  INTEGER ilaenv
217  REAL slamch, slansy
218  EXTERNAL ilaenv, lsame, slamch, slansy
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL slacpy, slascl, sormtr, sscal, sstedc, ssterf,
222  $ ssytrd, xerbla
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC max, sqrt
226 * ..
227 * .. Executable Statements ..
228 *
229 * Test the input parameters.
230 *
231  wantz = lsame( jobz, 'V' )
232  lower = lsame( uplo, 'L' )
233  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
234 *
235  info = 0
236  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
237  info = -1
238  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
239  info = -2
240  ELSE IF( n.LT.0 ) THEN
241  info = -3
242  ELSE IF( lda.LT.max( 1, n ) ) THEN
243  info = -5
244  END IF
245 *
246  IF( info.EQ.0 ) THEN
247  IF( n.LE.1 ) THEN
248  liwmin = 1
249  lwmin = 1
250  lopt = lwmin
251  liopt = liwmin
252  ELSE
253  IF( wantz ) THEN
254  liwmin = 3 + 5*n
255  lwmin = 1 + 6*n + 2*n**2
256  ELSE
257  liwmin = 1
258  lwmin = 2*n + 1
259  END IF
260  lopt = max( lwmin, 2*n +
261  $ ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 ) )
262  liopt = liwmin
263  END IF
264  work( 1 ) = lopt
265  iwork( 1 ) = liopt
266 *
267  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
268  info = -8
269  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
270  info = -10
271  END IF
272  END IF
273 *
274  IF( info.NE.0 ) THEN
275  CALL xerbla( 'SSYEVD', -info )
276  RETURN
277  ELSE IF( lquery ) THEN
278  RETURN
279  END IF
280 *
281 * Quick return if possible
282 *
283  IF( n.EQ.0 )
284  $ RETURN
285 *
286  IF( n.EQ.1 ) THEN
287  w( 1 ) = a( 1, 1 )
288  IF( wantz )
289  $ a( 1, 1 ) = one
290  RETURN
291  END IF
292 *
293 * Get machine constants.
294 *
295  safmin = slamch( 'Safe minimum' )
296  eps = slamch( 'Precision' )
297  smlnum = safmin / eps
298  bignum = one / smlnum
299  rmin = sqrt( smlnum )
300  rmax = sqrt( bignum )
301 *
302 * Scale matrix to allowable range, if necessary.
303 *
304  anrm = slansy( 'M', uplo, n, a, lda, work )
305  iscale = 0
306  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
307  iscale = 1
308  sigma = rmin / anrm
309  ELSE IF( anrm.GT.rmax ) THEN
310  iscale = 1
311  sigma = rmax / anrm
312  END IF
313  IF( iscale.EQ.1 )
314  $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
315 *
316 * Call SSYTRD to reduce symmetric matrix to tridiagonal form.
317 *
318  inde = 1
319  indtau = inde + n
320  indwrk = indtau + n
321  llwork = lwork - indwrk + 1
322  indwk2 = indwrk + n*n
323  llwrk2 = lwork - indwk2 + 1
324 *
325  CALL ssytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
326  $ work( indwrk ), llwork, iinfo )
327 *
328 * For eigenvalues only, call SSTERF. For eigenvectors, first call
329 * SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
330 * tridiagonal matrix, then call SORMTR to multiply it by the
331 * Householder transformations stored in A.
332 *
333  IF( .NOT.wantz ) THEN
334  CALL ssterf( n, w, work( inde ), info )
335  ELSE
336  CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
337  $ work( indwk2 ), llwrk2, iwork, liwork, info )
338  CALL sormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
339  $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
340  CALL slacpy( 'A', n, n, work( indwrk ), n, a, lda )
341  END IF
342 *
343 * If matrix was scaled, then rescale eigenvalues appropriately.
344 *
345  IF( iscale.EQ.1 )
346  $ CALL sscal( n, one / sigma, w, 1 )
347 *
348  work( 1 ) = lopt
349  iwork( 1 ) = liopt
350 *
351  RETURN
352 *
353 * End of SSYEVD
354 *
355  END