LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dsyev.f
Go to the documentation of this file.
1 *> \brief <b> DSYEV 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 DSYEV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyev.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyev.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyev.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER JOBZ, UPLO
25 * INTEGER INFO, LDA, LWORK, N
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DSYEV computes all eigenvalues and, optionally, eigenvectors of a
38 *> real symmetric matrix A.
39 *> \endverbatim
40 *
41 * Arguments:
42 * ==========
43 *
44 *> \param[in] JOBZ
45 *> \verbatim
46 *> JOBZ is CHARACTER*1
47 *> = 'N': Compute eigenvalues only;
48 *> = 'V': Compute eigenvalues and eigenvectors.
49 *> \endverbatim
50 *>
51 *> \param[in] UPLO
52 *> \verbatim
53 *> UPLO is CHARACTER*1
54 *> = 'U': Upper triangle of A is stored;
55 *> = 'L': Lower triangle of A is stored.
56 *> \endverbatim
57 *>
58 *> \param[in] N
59 *> \verbatim
60 *> N is INTEGER
61 *> The order of the matrix A. N >= 0.
62 *> \endverbatim
63 *>
64 *> \param[in,out] A
65 *> \verbatim
66 *> A is DOUBLE PRECISION array, dimension (LDA, N)
67 *> On entry, the symmetric matrix A. If UPLO = 'U', the
68 *> leading N-by-N upper triangular part of A contains the
69 *> upper triangular part of the matrix A. If UPLO = 'L',
70 *> the leading N-by-N lower triangular part of A contains
71 *> the lower triangular part of the matrix A.
72 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
73 *> orthonormal eigenvectors of the matrix A.
74 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
75 *> or the upper triangle (if UPLO='U') of A, including the
76 *> diagonal, is destroyed.
77 *> \endverbatim
78 *>
79 *> \param[in] LDA
80 *> \verbatim
81 *> LDA is INTEGER
82 *> The leading dimension of the array A. LDA >= max(1,N).
83 *> \endverbatim
84 *>
85 *> \param[out] W
86 *> \verbatim
87 *> W is DOUBLE PRECISION array, dimension (N)
88 *> If INFO = 0, the eigenvalues in ascending order.
89 *> \endverbatim
90 *>
91 *> \param[out] WORK
92 *> \verbatim
93 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
94 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
95 *> \endverbatim
96 *>
97 *> \param[in] LWORK
98 *> \verbatim
99 *> LWORK is INTEGER
100 *> The length of the array WORK. LWORK >= max(1,3*N-1).
101 *> For optimal efficiency, LWORK >= (NB+2)*N,
102 *> where NB is the blocksize for DSYTRD returned by ILAENV.
103 *>
104 *> If LWORK = -1, then a workspace query is assumed; the routine
105 *> only calculates the optimal size of the WORK array, returns
106 *> this value as the first entry of the WORK array, and no error
107 *> message related to LWORK is issued by XERBLA.
108 *> \endverbatim
109 *>
110 *> \param[out] INFO
111 *> \verbatim
112 *> INFO is INTEGER
113 *> = 0: successful exit
114 *> < 0: if INFO = -i, the i-th argument had an illegal value
115 *> > 0: if INFO = i, the algorithm failed to converge; i
116 *> off-diagonal elements of an intermediate tridiagonal
117 *> form did not converge to zero.
118 *> \endverbatim
119 *
120 * Authors:
121 * ========
122 *
123 *> \author Univ. of Tennessee
124 *> \author Univ. of California Berkeley
125 *> \author Univ. of Colorado Denver
126 *> \author NAG Ltd.
127 *
128 *> \date November 2011
129 *
130 *> \ingroup doubleSYeigen
131 *
132 * =====================================================================
133  SUBROUTINE dsyev( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
134 *
135 * -- LAPACK driver routine (version 3.4.0) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * November 2011
139 *
140 * .. Scalar Arguments ..
141  CHARACTER jobz, uplo
142  INTEGER info, lda, lwork, n
143 * ..
144 * .. Array Arguments ..
145  DOUBLE PRECISION a( lda, * ), w( * ), work( * )
146 * ..
147 *
148 * =====================================================================
149 *
150 * .. Parameters ..
151  DOUBLE PRECISION zero, one
152  parameter( zero = 0.0d0, one = 1.0d0 )
153 * ..
154 * .. Local Scalars ..
155  LOGICAL lower, lquery, wantz
156  INTEGER iinfo, imax, inde, indtau, indwrk, iscale,
157  $ llwork, lwkopt, nb
158  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
159  $ smlnum
160 * ..
161 * .. External Functions ..
162  LOGICAL lsame
163  INTEGER ilaenv
164  DOUBLE PRECISION dlamch, dlansy
165  EXTERNAL lsame, ilaenv, dlamch, dlansy
166 * ..
167 * .. External Subroutines ..
168  EXTERNAL dlascl, dorgtr, dscal, dsteqr, dsterf, dsytrd,
169  $ xerbla
170 * ..
171 * .. Intrinsic Functions ..
172  INTRINSIC max, sqrt
173 * ..
174 * .. Executable Statements ..
175 *
176 * Test the input parameters.
177 *
178  wantz = lsame( jobz, 'V' )
179  lower = lsame( uplo, 'L' )
180  lquery = ( lwork.EQ.-1 )
181 *
182  info = 0
183  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
184  info = -1
185  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
186  info = -2
187  ELSE IF( n.LT.0 ) THEN
188  info = -3
189  ELSE IF( lda.LT.max( 1, n ) ) THEN
190  info = -5
191  END IF
192 *
193  IF( info.EQ.0 ) THEN
194  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
195  lwkopt = max( 1, ( nb+2 )*n )
196  work( 1 ) = lwkopt
197 *
198  IF( lwork.LT.max( 1, 3*n-1 ) .AND. .NOT.lquery )
199  $ info = -8
200  END IF
201 *
202  IF( info.NE.0 ) THEN
203  CALL xerbla( 'DSYEV ', -info )
204  RETURN
205  ELSE IF( lquery ) THEN
206  RETURN
207  END IF
208 *
209 * Quick return if possible
210 *
211  IF( n.EQ.0 ) THEN
212  RETURN
213  END IF
214 *
215  IF( n.EQ.1 ) THEN
216  w( 1 ) = a( 1, 1 )
217  work( 1 ) = 2
218  IF( wantz )
219  $ a( 1, 1 ) = one
220  RETURN
221  END IF
222 *
223 * Get machine constants.
224 *
225  safmin = dlamch( 'Safe minimum' )
226  eps = dlamch( 'Precision' )
227  smlnum = safmin / eps
228  bignum = one / smlnum
229  rmin = sqrt( smlnum )
230  rmax = sqrt( bignum )
231 *
232 * Scale matrix to allowable range, if necessary.
233 *
234  anrm = dlansy( 'M', uplo, n, a, lda, work )
235  iscale = 0
236  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
237  iscale = 1
238  sigma = rmin / anrm
239  ELSE IF( anrm.GT.rmax ) THEN
240  iscale = 1
241  sigma = rmax / anrm
242  END IF
243  IF( iscale.EQ.1 )
244  $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
245 *
246 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
247 *
248  inde = 1
249  indtau = inde + n
250  indwrk = indtau + n
251  llwork = lwork - indwrk + 1
252  CALL dsytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
253  $ work( indwrk ), llwork, iinfo )
254 *
255 * For eigenvalues only, call DSTERF. For eigenvectors, first call
256 * DORGTR to generate the orthogonal matrix, then call DSTEQR.
257 *
258  IF( .NOT.wantz ) THEN
259  CALL dsterf( n, w, work( inde ), info )
260  ELSE
261  CALL dorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
262  $ llwork, iinfo )
263  CALL dsteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
264  $ info )
265  END IF
266 *
267 * If matrix was scaled, then rescale eigenvalues appropriately.
268 *
269  IF( iscale.EQ.1 ) THEN
270  IF( info.EQ.0 ) THEN
271  imax = n
272  ELSE
273  imax = info - 1
274  END IF
275  CALL dscal( imax, one / sigma, w, 1 )
276  END IF
277 *
278 * Set WORK(1) to optimal workspace size.
279 *
280  work( 1 ) = lwkopt
281 *
282  RETURN
283 *
284 * End of DSYEV
285 *
286  END