LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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
9 *> Download DSYEVD + dependencies
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 *> 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, DSYEVD needs N**2 more
51 *> workspace than DSYEVX.
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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
101 *> If INFO = 0, the eigenvalues in ascending order.
102 *> \endverbatim
103 *>
104 *> \param[out] WORK
105 *> \verbatim
106 *> WORK is DOUBLE PRECISION 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 doubleSYeigen
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 *>
182 * =====================================================================
183  SUBROUTINE dsyevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
184  $ LIWORK, INFO )
185 *
186 * -- LAPACK driver routine --
187 * -- LAPACK is a software package provided by Univ. of Tennessee, --
188 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189 *
190 * .. Scalar Arguments ..
191  CHARACTER JOBZ, UPLO
192  INTEGER INFO, LDA, LIWORK, LWORK, N
193 * ..
194 * .. Array Arguments ..
195  INTEGER IWORK( * )
196  DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
197 * ..
198 *
199 * =====================================================================
200 *
201 * .. Parameters ..
202  DOUBLE PRECISION ZERO, ONE
203  parameter( zero = 0.0d+0, one = 1.0d+0 )
204 * ..
205 * .. Local Scalars ..
206 *
207  LOGICAL LOWER, LQUERY, WANTZ
208  INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
209  $ liopt, liwmin, llwork, llwrk2, lopt, lwmin
210  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
211  $ smlnum
212 * ..
213 * .. External Functions ..
214  LOGICAL LSAME
215  INTEGER ILAENV
216  DOUBLE PRECISION DLAMCH, DLANSY
217  EXTERNAL lsame, dlamch, dlansy, ilaenv
218 * ..
219 * .. External Subroutines ..
220  EXTERNAL dlacpy, dlascl, dormtr, dscal, dstedc, dsterf,
221  $ dsytrd, xerbla
222 * ..
223 * .. Intrinsic Functions ..
224  INTRINSIC max, sqrt
225 * ..
226 * .. Executable Statements ..
227 *
228 * Test the input parameters.
229 *
230  wantz = lsame( jobz, 'V' )
231  lower = lsame( uplo, 'L' )
232  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
233 *
234  info = 0
235  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
236  info = -1
237  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
238  info = -2
239  ELSE IF( n.LT.0 ) THEN
240  info = -3
241  ELSE IF( lda.LT.max( 1, n ) ) THEN
242  info = -5
243  END IF
244 *
245  IF( info.EQ.0 ) THEN
246  IF( n.LE.1 ) THEN
247  liwmin = 1
248  lwmin = 1
249  lopt = lwmin
250  liopt = liwmin
251  ELSE
252  IF( wantz ) THEN
253  liwmin = 3 + 5*n
254  lwmin = 1 + 6*n + 2*n**2
255  ELSE
256  liwmin = 1
257  lwmin = 2*n + 1
258  END IF
259  lopt = max( lwmin, 2*n +
260  $ ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 ) )
261  liopt = liwmin
262  END IF
263  work( 1 ) = lopt
264  iwork( 1 ) = liopt
265 *
266  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
267  info = -8
268  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
269  info = -10
270  END IF
271  END IF
272 *
273  IF( info.NE.0 ) THEN
274  CALL xerbla( 'DSYEVD', -info )
275  RETURN
276  ELSE IF( lquery ) THEN
277  RETURN
278  END IF
279 *
280 * Quick return if possible
281 *
282  IF( n.EQ.0 )
283  $ RETURN
284 *
285  IF( n.EQ.1 ) THEN
286  w( 1 ) = a( 1, 1 )
287  IF( wantz )
288  $ a( 1, 1 ) = one
289  RETURN
290  END IF
291 *
292 * Get machine constants.
293 *
294  safmin = dlamch( 'Safe minimum' )
295  eps = dlamch( 'Precision' )
296  smlnum = safmin / eps
297  bignum = one / smlnum
298  rmin = sqrt( smlnum )
299  rmax = sqrt( bignum )
300 *
301 * Scale matrix to allowable range, if necessary.
302 *
303  anrm = dlansy( 'M', uplo, n, a, lda, work )
304  iscale = 0
305  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
306  iscale = 1
307  sigma = rmin / anrm
308  ELSE IF( anrm.GT.rmax ) THEN
309  iscale = 1
310  sigma = rmax / anrm
311  END IF
312  IF( iscale.EQ.1 )
313  $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
314 *
315 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
316 *
317  inde = 1
318  indtau = inde + n
319  indwrk = indtau + n
320  llwork = lwork - indwrk + 1
321  indwk2 = indwrk + n*n
322  llwrk2 = lwork - indwk2 + 1
323 *
324  CALL dsytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
325  $ work( indwrk ), llwork, iinfo )
326 *
327 * For eigenvalues only, call DSTERF. For eigenvectors, first call
328 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
329 * tridiagonal matrix, then call DORMTR to multiply it by the
330 * Householder transformations stored in A.
331 *
332  IF( .NOT.wantz ) THEN
333  CALL dsterf( n, w, work( inde ), info )
334  ELSE
335  CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
336  $ work( indwk2 ), llwrk2, iwork, liwork, info )
337  CALL dormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
338  $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
339  CALL dlacpy( 'A', n, n, work( indwrk ), n, a, lda )
340  END IF
341 *
342 * If matrix was scaled, then rescale eigenvalues appropriately.
343 *
344  IF( iscale.EQ.1 )
345  $ CALL dscal( n, one / sigma, w, 1 )
346 *
347  work( 1 ) = lopt
348  iwork( 1 ) = liopt
349 *
350  RETURN
351 *
352 * End of DSYEVD
353 *
354  END
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 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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO)
DSTEDC
Definition: dstedc.f:188
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMTR
Definition: dormtr.f:171
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
DSYTRD
Definition: dsytrd.f:192
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:185