1 SUBROUTINE DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
2 $ IWORK, LIWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSPEVD computes all the eigenvalues and, optionally, eigenvectors
22 * of a real symmetric matrix A in packed storage. If eigenvectors are
23 * desired, it uses a divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
31 *
32 * Arguments
33 * =========
34 *
35 * JOBZ (input) CHARACTER*1
36 * = 'N': Compute eigenvalues only;
37 * = 'V': Compute eigenvalues and eigenvectors.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of A is stored;
41 * = 'L': Lower triangle of A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
47 * On entry, the upper or lower triangle of the symmetric matrix
48 * A, packed columnwise in a linear array. The j-th column of A
49 * is stored in the array AP as follows:
50 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
51 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
52 *
53 * On exit, AP is overwritten by values generated during the
54 * reduction to tridiagonal form. If UPLO = 'U', the diagonal
55 * and first superdiagonal of the tridiagonal matrix T overwrite
56 * the corresponding elements of A, and if UPLO = 'L', the
57 * diagonal and first subdiagonal of T overwrite the
58 * corresponding elements of A.
59 *
60 * W (output) DOUBLE PRECISION array, dimension (N)
61 * If INFO = 0, the eigenvalues in ascending order.
62 *
63 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
64 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
65 * eigenvectors of the matrix A, with the i-th column of Z
66 * holding the eigenvector associated with W(i).
67 * If JOBZ = 'N', then Z is not referenced.
68 *
69 * LDZ (input) INTEGER
70 * The leading dimension of the array Z. LDZ >= 1, and if
71 * JOBZ = 'V', LDZ >= max(1,N).
72 *
73 * WORK (workspace/output) DOUBLE PRECISION array,
74 * dimension (LWORK)
75 * On exit, if INFO = 0, WORK(1) returns the required LWORK.
76 *
77 * LWORK (input) INTEGER
78 * The dimension of the array WORK.
79 * If N <= 1, LWORK must be at least 1.
80 * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
81 * If JOBZ = 'V' and N > 1, LWORK must be at least
82 * 1 + 6*N + N**2.
83 *
84 * If LWORK = -1, then a workspace query is assumed; the routine
85 * only calculates the required sizes of the WORK and IWORK
86 * arrays, returns these values as the first entries of the WORK
87 * and IWORK arrays, and no error message related to LWORK or
88 * LIWORK is issued by XERBLA.
89 *
90 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
91 * On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
92 *
93 * LIWORK (input) INTEGER
94 * The dimension of the array IWORK.
95 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
96 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
97 *
98 * If LIWORK = -1, then a workspace query is assumed; the
99 * routine only calculates the required sizes of the WORK and
100 * IWORK arrays, returns these values as the first entries of
101 * the WORK and IWORK arrays, and no error message related to
102 * LWORK or LIWORK is issued by XERBLA.
103 *
104 * INFO (output) INTEGER
105 * = 0: successful exit
106 * < 0: if INFO = -i, the i-th argument had an illegal value.
107 * > 0: if INFO = i, the algorithm failed to converge; i
108 * off-diagonal elements of an intermediate tridiagonal
109 * form did not converge to zero.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 DOUBLE PRECISION ZERO, ONE
115 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
116 * ..
117 * .. Local Scalars ..
118 LOGICAL LQUERY, WANTZ
119 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
120 $ LLWORK, LWMIN
121 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
122 $ SMLNUM
123 * ..
124 * .. External Functions ..
125 LOGICAL LSAME
126 DOUBLE PRECISION DLAMCH, DLANSP
127 EXTERNAL LSAME, DLAMCH, DLANSP
128 * ..
129 * .. External Subroutines ..
130 EXTERNAL DOPMTR, DSCAL, DSPTRD, DSTEDC, DSTERF, XERBLA
131 * ..
132 * .. Intrinsic Functions ..
133 INTRINSIC SQRT
134 * ..
135 * .. Executable Statements ..
136 *
137 * Test the input parameters.
138 *
139 WANTZ = LSAME( JOBZ, 'V' )
140 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
141 *
142 INFO = 0
143 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
144 INFO = -1
145 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
146 $ THEN
147 INFO = -2
148 ELSE IF( N.LT.0 ) THEN
149 INFO = -3
150 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
151 INFO = -7
152 END IF
153 *
154 IF( INFO.EQ.0 ) THEN
155 IF( N.LE.1 ) THEN
156 LIWMIN = 1
157 LWMIN = 1
158 ELSE
159 IF( WANTZ ) THEN
160 LIWMIN = 3 + 5*N
161 LWMIN = 1 + 6*N + N**2
162 ELSE
163 LIWMIN = 1
164 LWMIN = 2*N
165 END IF
166 END IF
167 IWORK( 1 ) = LIWMIN
168 WORK( 1 ) = LWMIN
169 *
170 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
171 INFO = -9
172 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
173 INFO = -11
174 END IF
175 END IF
176 *
177 IF( INFO.NE.0 ) THEN
178 CALL XERBLA( 'DSPEVD', -INFO )
179 RETURN
180 ELSE IF( LQUERY ) THEN
181 RETURN
182 END IF
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 IF( N.EQ.1 ) THEN
190 W( 1 ) = AP( 1 )
191 IF( WANTZ )
192 $ Z( 1, 1 ) = ONE
193 RETURN
194 END IF
195 *
196 * Get machine constants.
197 *
198 SAFMIN = DLAMCH( 'Safe minimum' )
199 EPS = DLAMCH( 'Precision' )
200 SMLNUM = SAFMIN / EPS
201 BIGNUM = ONE / SMLNUM
202 RMIN = SQRT( SMLNUM )
203 RMAX = SQRT( BIGNUM )
204 *
205 * Scale matrix to allowable range, if necessary.
206 *
207 ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
208 ISCALE = 0
209 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
210 ISCALE = 1
211 SIGMA = RMIN / ANRM
212 ELSE IF( ANRM.GT.RMAX ) THEN
213 ISCALE = 1
214 SIGMA = RMAX / ANRM
215 END IF
216 IF( ISCALE.EQ.1 ) THEN
217 CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
218 END IF
219 *
220 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
221 *
222 INDE = 1
223 INDTAU = INDE + N
224 CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
225 *
226 * For eigenvalues only, call DSTERF. For eigenvectors, first call
227 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
228 * tridiagonal matrix, then call DOPMTR to multiply it by the
229 * Householder transformations represented in AP.
230 *
231 IF( .NOT.WANTZ ) THEN
232 CALL DSTERF( N, W, WORK( INDE ), INFO )
233 ELSE
234 INDWRK = INDTAU + N
235 LLWORK = LWORK - INDWRK + 1
236 CALL DSTEDC( 'I', N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
237 $ LLWORK, IWORK, LIWORK, INFO )
238 CALL DOPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ,
239 $ WORK( INDWRK ), IINFO )
240 END IF
241 *
242 * If matrix was scaled, then rescale eigenvalues appropriately.
243 *
244 IF( ISCALE.EQ.1 )
245 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
246 *
247 WORK( 1 ) = LWMIN
248 IWORK( 1 ) = LIWMIN
249 RETURN
250 *
251 * End of DSPEVD
252 *
253 END
2 $ IWORK, LIWORK, INFO )
3 *
4 * -- LAPACK driver routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER JOBZ, UPLO
11 INTEGER INFO, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSPEVD computes all the eigenvalues and, optionally, eigenvectors
22 * of a real symmetric matrix A in packed storage. If eigenvectors are
23 * desired, it uses a divide and conquer algorithm.
24 *
25 * The divide and conquer algorithm makes very mild assumptions about
26 * floating point arithmetic. It will work on machines with a guard
27 * digit in add/subtract, or on those binary machines without guard
28 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
29 * Cray-2. It could conceivably fail on hexadecimal or decimal machines
30 * without guard digits, but we know of none.
31 *
32 * Arguments
33 * =========
34 *
35 * JOBZ (input) CHARACTER*1
36 * = 'N': Compute eigenvalues only;
37 * = 'V': Compute eigenvalues and eigenvectors.
38 *
39 * UPLO (input) CHARACTER*1
40 * = 'U': Upper triangle of A is stored;
41 * = 'L': Lower triangle of A is stored.
42 *
43 * N (input) INTEGER
44 * The order of the matrix A. N >= 0.
45 *
46 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
47 * On entry, the upper or lower triangle of the symmetric matrix
48 * A, packed columnwise in a linear array. The j-th column of A
49 * is stored in the array AP as follows:
50 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
51 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
52 *
53 * On exit, AP is overwritten by values generated during the
54 * reduction to tridiagonal form. If UPLO = 'U', the diagonal
55 * and first superdiagonal of the tridiagonal matrix T overwrite
56 * the corresponding elements of A, and if UPLO = 'L', the
57 * diagonal and first subdiagonal of T overwrite the
58 * corresponding elements of A.
59 *
60 * W (output) DOUBLE PRECISION array, dimension (N)
61 * If INFO = 0, the eigenvalues in ascending order.
62 *
63 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
64 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
65 * eigenvectors of the matrix A, with the i-th column of Z
66 * holding the eigenvector associated with W(i).
67 * If JOBZ = 'N', then Z is not referenced.
68 *
69 * LDZ (input) INTEGER
70 * The leading dimension of the array Z. LDZ >= 1, and if
71 * JOBZ = 'V', LDZ >= max(1,N).
72 *
73 * WORK (workspace/output) DOUBLE PRECISION array,
74 * dimension (LWORK)
75 * On exit, if INFO = 0, WORK(1) returns the required LWORK.
76 *
77 * LWORK (input) INTEGER
78 * The dimension of the array WORK.
79 * If N <= 1, LWORK must be at least 1.
80 * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N.
81 * If JOBZ = 'V' and N > 1, LWORK must be at least
82 * 1 + 6*N + N**2.
83 *
84 * If LWORK = -1, then a workspace query is assumed; the routine
85 * only calculates the required sizes of the WORK and IWORK
86 * arrays, returns these values as the first entries of the WORK
87 * and IWORK arrays, and no error message related to LWORK or
88 * LIWORK is issued by XERBLA.
89 *
90 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
91 * On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
92 *
93 * LIWORK (input) INTEGER
94 * The dimension of the array IWORK.
95 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
96 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
97 *
98 * If LIWORK = -1, then a workspace query is assumed; the
99 * routine only calculates the required sizes of the WORK and
100 * IWORK arrays, returns these values as the first entries of
101 * the WORK and IWORK arrays, and no error message related to
102 * LWORK or LIWORK is issued by XERBLA.
103 *
104 * INFO (output) INTEGER
105 * = 0: successful exit
106 * < 0: if INFO = -i, the i-th argument had an illegal value.
107 * > 0: if INFO = i, the algorithm failed to converge; i
108 * off-diagonal elements of an intermediate tridiagonal
109 * form did not converge to zero.
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114 DOUBLE PRECISION ZERO, ONE
115 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
116 * ..
117 * .. Local Scalars ..
118 LOGICAL LQUERY, WANTZ
119 INTEGER IINFO, INDE, INDTAU, INDWRK, ISCALE, LIWMIN,
120 $ LLWORK, LWMIN
121 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
122 $ SMLNUM
123 * ..
124 * .. External Functions ..
125 LOGICAL LSAME
126 DOUBLE PRECISION DLAMCH, DLANSP
127 EXTERNAL LSAME, DLAMCH, DLANSP
128 * ..
129 * .. External Subroutines ..
130 EXTERNAL DOPMTR, DSCAL, DSPTRD, DSTEDC, DSTERF, XERBLA
131 * ..
132 * .. Intrinsic Functions ..
133 INTRINSIC SQRT
134 * ..
135 * .. Executable Statements ..
136 *
137 * Test the input parameters.
138 *
139 WANTZ = LSAME( JOBZ, 'V' )
140 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
141 *
142 INFO = 0
143 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
144 INFO = -1
145 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
146 $ THEN
147 INFO = -2
148 ELSE IF( N.LT.0 ) THEN
149 INFO = -3
150 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
151 INFO = -7
152 END IF
153 *
154 IF( INFO.EQ.0 ) THEN
155 IF( N.LE.1 ) THEN
156 LIWMIN = 1
157 LWMIN = 1
158 ELSE
159 IF( WANTZ ) THEN
160 LIWMIN = 3 + 5*N
161 LWMIN = 1 + 6*N + N**2
162 ELSE
163 LIWMIN = 1
164 LWMIN = 2*N
165 END IF
166 END IF
167 IWORK( 1 ) = LIWMIN
168 WORK( 1 ) = LWMIN
169 *
170 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
171 INFO = -9
172 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
173 INFO = -11
174 END IF
175 END IF
176 *
177 IF( INFO.NE.0 ) THEN
178 CALL XERBLA( 'DSPEVD', -INFO )
179 RETURN
180 ELSE IF( LQUERY ) THEN
181 RETURN
182 END IF
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 IF( N.EQ.1 ) THEN
190 W( 1 ) = AP( 1 )
191 IF( WANTZ )
192 $ Z( 1, 1 ) = ONE
193 RETURN
194 END IF
195 *
196 * Get machine constants.
197 *
198 SAFMIN = DLAMCH( 'Safe minimum' )
199 EPS = DLAMCH( 'Precision' )
200 SMLNUM = SAFMIN / EPS
201 BIGNUM = ONE / SMLNUM
202 RMIN = SQRT( SMLNUM )
203 RMAX = SQRT( BIGNUM )
204 *
205 * Scale matrix to allowable range, if necessary.
206 *
207 ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
208 ISCALE = 0
209 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
210 ISCALE = 1
211 SIGMA = RMIN / ANRM
212 ELSE IF( ANRM.GT.RMAX ) THEN
213 ISCALE = 1
214 SIGMA = RMAX / ANRM
215 END IF
216 IF( ISCALE.EQ.1 ) THEN
217 CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
218 END IF
219 *
220 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
221 *
222 INDE = 1
223 INDTAU = INDE + N
224 CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
225 *
226 * For eigenvalues only, call DSTERF. For eigenvectors, first call
227 * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
228 * tridiagonal matrix, then call DOPMTR to multiply it by the
229 * Householder transformations represented in AP.
230 *
231 IF( .NOT.WANTZ ) THEN
232 CALL DSTERF( N, W, WORK( INDE ), INFO )
233 ELSE
234 INDWRK = INDTAU + N
235 LLWORK = LWORK - INDWRK + 1
236 CALL DSTEDC( 'I', N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
237 $ LLWORK, IWORK, LIWORK, INFO )
238 CALL DOPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ,
239 $ WORK( INDWRK ), IINFO )
240 END IF
241 *
242 * If matrix was scaled, then rescale eigenvalues appropriately.
243 *
244 IF( ISCALE.EQ.1 )
245 $ CALL DSCAL( N, ONE / SIGMA, W, 1 )
246 *
247 WORK( 1 ) = LWMIN
248 IWORK( 1 ) = LIWMIN
249 RETURN
250 *
251 * End of DSPEVD
252 *
253 END