KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsyev.f
Go to the documentation of this file.
1  SUBROUTINE dsyev( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
2 *
3 * -- LAPACK driver routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * June 30, 1999
7 *
8 * .. Scalar Arguments ..
9  CHARACTER JOBZ, UPLO
10  INTEGER INFO, LDA, LWORK, N
11 * ..
12 * .. Array Arguments ..
13  DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSYEV computes all eigenvalues and, optionally, eigenvectors of a
20 * real symmetric matrix A.
21 *
22 * Arguments
23 * =========
24 *
25 * JOBZ (input) CHARACTER*1
26 * = 'N': Compute eigenvalues only;
27 * = 'V': Compute eigenvalues and eigenvectors.
28 *
29 * UPLO (input) CHARACTER*1
30 * = 'U': Upper triangle of A is stored;
31 * = 'L': Lower triangle of A is stored.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
37 * On entry, the symmetric matrix A. If UPLO = 'U', the
38 * leading N-by-N upper triangular part of A contains the
39 * upper triangular part of the matrix A. If UPLO = 'L',
40 * the leading N-by-N lower triangular part of A contains
41 * the lower triangular part of the matrix A.
42 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
43 * orthonormal eigenvectors of the matrix A.
44 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
45 * or the upper triangle (if UPLO='U') of A, including the
46 * diagonal, is destroyed.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,N).
50 *
51 * W (output) DOUBLE PRECISION array, dimension (N)
52 * If INFO = 0, the eigenvalues in ascending order.
53 *
54 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
55 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
56 *
57 * LWORK (input) INTEGER
58 * The length of the array WORK. LWORK >= max(1,3*N-1).
59 * For optimal efficiency, LWORK >= (NB+2)*N,
60 * where NB is the blocksize for DSYTRD returned by ILAENV.
61 *
62 * If LWORK = -1, then a workspace query is assumed; the routine
63 * only calculates the optimal size of the WORK array, returns
64 * this value as the first entry of the WORK array, and no error
65 * message related to LWORK is issued by XERBLA.
66 *
67 * INFO (output) INTEGER
68 * = 0: successful exit
69 * < 0: if INFO = -i, the i-th argument had an illegal value
70 * > 0: if INFO = i, the algorithm failed to converge; i
71 * off-diagonal elements of an intermediate tridiagonal
72 * form did not converge to zero.
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77  DOUBLE PRECISION ZERO, ONE
78  parameter( zero = 0.0d0, one = 1.0d0 )
79 * ..
80 * .. Local Scalars ..
81  LOGICAL LOWER, LQUERY, WANTZ
82  INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
83  $ LLWORK, LOPT, LWKOPT, NB
84  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
85  $ SMLNUM
86 * ..
87 * .. External Functions ..
88  LOGICAL LSAME
89  INTEGER ILAENV
90  DOUBLE PRECISION DLAMCH, DLANSY
91  EXTERNAL lsame, ilaenv, dlamch, dlansy
92 * ..
93 * .. External Subroutines ..
94  EXTERNAL dlascl, dorgtr, dscal, dsteqr, dsterf, dsytrd,
95  $ xerbla
96 * ..
97 * .. Intrinsic Functions ..
98  INTRINSIC max, sqrt
99 * ..
100 * .. Executable Statements ..
101 *
102 * Test the input parameters.
103 *
104  wantz = lsame( jobz, 'V' )
105  lower = lsame( uplo, 'L' )
106  lquery = ( lwork.EQ.-1 )
107 *
108  info = 0
109  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
110  info = -1
111  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
112  info = -2
113  ELSE IF( n.LT.0 ) THEN
114  info = -3
115  ELSE IF( lda.LT.max( 1, n ) ) THEN
116  info = -5
117  ELSE IF( lwork.LT.max( 1, 3*n-1 ) .AND. .NOT.lquery ) THEN
118  info = -8
119  END IF
120 *
121  IF( info.EQ.0 ) THEN
122  nb = ilaenv( 1, 'DSYTRD', uplo, n, -1, -1, -1 )
123  lwkopt = max( 1, ( nb+2 )*n )
124  work( 1 ) = lwkopt
125  END IF
126 *
127  IF( info.NE.0 ) THEN
128  CALL xerbla( 'DSYEV ', -info )
129  RETURN
130  ELSE IF( lquery ) THEN
131  RETURN
132  END IF
133 *
134 * Quick return if possible
135 *
136  IF( n.EQ.0 ) THEN
137  work( 1 ) = 1
138  RETURN
139  END IF
140 *
141  IF( n.EQ.1 ) THEN
142  w( 1 ) = a( 1, 1 )
143  work( 1 ) = 3
144  IF( wantz )
145  $ a( 1, 1 ) = one
146  RETURN
147  END IF
148 *
149 * Get machine constants.
150 *
151  safmin = dlamch( 'Safe minimum' )
152  eps = dlamch( 'Precision' )
153  smlnum = safmin / eps
154  bignum = one / smlnum
155  rmin = sqrt( smlnum )
156  rmax = sqrt( bignum )
157 *
158 * Scale matrix to allowable range, if necessary.
159 *
160  anrm = dlansy( 'M', uplo, n, a, lda, work )
161  iscale = 0
162  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
163  iscale = 1
164  sigma = rmin / anrm
165  ELSE IF( anrm.GT.rmax ) THEN
166  iscale = 1
167  sigma = rmax / anrm
168  END IF
169  IF( iscale.EQ.1 )
170  $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
171 *
172 * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
173 *
174  inde = 1
175  indtau = inde + n
176  indwrk = indtau + n
177  llwork = lwork - indwrk + 1
178  CALL dsytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
179  $ work( indwrk ), llwork, iinfo )
180  lopt = 2*n + work( indwrk )
181 *
182 * For eigenvalues only, call DSTERF. For eigenvectors, first call
183 * DORGTR to generate the orthogonal matrix, then call DSTEQR.
184 *
185  IF( .NOT.wantz ) THEN
186  CALL dsterf( n, w, work( inde ), info )
187  ELSE
188  CALL dorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
189  $ llwork, iinfo )
190  CALL dsteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
191  $ info )
192  END IF
193 *
194 * If matrix was scaled, then rescale eigenvalues appropriately.
195 *
196  IF( iscale.EQ.1 ) THEN
197  IF( info.EQ.0 ) THEN
198  imax = n
199  ELSE
200  imax = info - 1
201  END IF
202  CALL dscal( imax, one / sigma, w, 1 )
203  END IF
204 *
205 * Set WORK(1) to optimal workspace size.
206 *
207  work( 1 ) = lwkopt
208 *
209  RETURN
210 *
211 * End of DSYEV
212 *
213  END
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
Definition: dlascl.f:2
subroutine dorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
Definition: dorgtr.f:2
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2
subroutine dsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
Definition: dsteqr.f:2
subroutine dsterf(N, D, E, INFO)
Definition: dsterf.f:2
subroutine dsyev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)
Definition: dsyev.f:2
subroutine dsytrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
Definition: dsytrd.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2