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