1 SUBROUTINE DSTEVD( JOBZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
2 $ 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
11 INTEGER INFO, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSTEVD computes all eigenvalues and, optionally, eigenvectors of a
22 * real symmetric tridiagonal matrix. If eigenvectors are desired, it
23 * 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 * N (input) INTEGER
40 * The order of the matrix. N >= 0.
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry, the n diagonal elements of the tridiagonal matrix
44 * A.
45 * On exit, if INFO = 0, the eigenvalues in ascending order.
46 *
47 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
48 * On entry, the (n-1) subdiagonal elements of the tridiagonal
49 * matrix A, stored in elements 1 to N-1 of E.
50 * On exit, the contents of E are destroyed.
51 *
52 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
53 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
54 * eigenvectors of the matrix A, with the i-th column of Z
55 * holding the eigenvector associated with D(i).
56 * If JOBZ = 'N', then Z is not referenced.
57 *
58 * LDZ (input) INTEGER
59 * The leading dimension of the array Z. LDZ >= 1, and if
60 * JOBZ = 'V', LDZ >= max(1,N).
61 *
62 * WORK (workspace/output) DOUBLE PRECISION array,
63 * dimension (LWORK)
64 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
65 *
66 * LWORK (input) INTEGER
67 * The dimension of the array WORK.
68 * If JOBZ = 'N' or N <= 1 then LWORK must be at least 1.
69 * If JOBZ = 'V' and N > 1 then LWORK must be at least
70 * ( 1 + 4*N + N**2 ).
71 *
72 * If LWORK = -1, then a workspace query is assumed; the routine
73 * only calculates the optimal sizes of the WORK and IWORK
74 * arrays, returns these values as the first entries of the WORK
75 * and IWORK arrays, and no error message related to LWORK or
76 * LIWORK is issued by XERBLA.
77 *
78 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
79 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
80 *
81 * LIWORK (input) INTEGER
82 * The dimension of the array IWORK.
83 * If JOBZ = 'N' or N <= 1 then LIWORK must be at least 1.
84 * If JOBZ = 'V' and N > 1 then LIWORK must be at least 3+5*N.
85 *
86 * If LIWORK = -1, then a workspace query is assumed; the
87 * routine only calculates the optimal sizes of the WORK and
88 * IWORK arrays, returns these values as the first entries of
89 * the WORK and IWORK arrays, and no error message related to
90 * LWORK or LIWORK is issued by XERBLA.
91 *
92 * INFO (output) INTEGER
93 * = 0: successful exit
94 * < 0: if INFO = -i, the i-th argument had an illegal value
95 * > 0: if INFO = i, the algorithm failed to converge; i
96 * off-diagonal elements of E did not converge to zero.
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 DOUBLE PRECISION ZERO, ONE
102 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
103 * ..
104 * .. Local Scalars ..
105 LOGICAL LQUERY, WANTZ
106 INTEGER ISCALE, LIWMIN, LWMIN
107 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
108 $ TNRM
109 * ..
110 * .. External Functions ..
111 LOGICAL LSAME
112 DOUBLE PRECISION DLAMCH, DLANST
113 EXTERNAL LSAME, DLAMCH, DLANST
114 * ..
115 * .. External Subroutines ..
116 EXTERNAL DSCAL, DSTEDC, DSTERF, XERBLA
117 * ..
118 * .. Intrinsic Functions ..
119 INTRINSIC SQRT
120 * ..
121 * .. Executable Statements ..
122 *
123 * Test the input parameters.
124 *
125 WANTZ = LSAME( JOBZ, 'V' )
126 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
127 *
128 INFO = 0
129 LIWMIN = 1
130 LWMIN = 1
131 IF( N.GT.1 .AND. WANTZ ) THEN
132 LWMIN = 1 + 4*N + N**2
133 LIWMIN = 3 + 5*N
134 END IF
135 *
136 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
137 INFO = -1
138 ELSE IF( N.LT.0 ) THEN
139 INFO = -2
140 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
141 INFO = -6
142 END IF
143 *
144 IF( INFO.EQ.0 ) THEN
145 WORK( 1 ) = LWMIN
146 IWORK( 1 ) = LIWMIN
147 *
148 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
149 INFO = -8
150 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
151 INFO = -10
152 END IF
153 END IF
154 *
155 IF( INFO.NE.0 ) THEN
156 CALL XERBLA( 'DSTEVD', -INFO )
157 RETURN
158 ELSE IF( LQUERY ) THEN
159 RETURN
160 END IF
161 *
162 * Quick return if possible
163 *
164 IF( N.EQ.0 )
165 $ RETURN
166 *
167 IF( N.EQ.1 ) THEN
168 IF( WANTZ )
169 $ Z( 1, 1 ) = ONE
170 RETURN
171 END IF
172 *
173 * Get machine constants.
174 *
175 SAFMIN = DLAMCH( 'Safe minimum' )
176 EPS = DLAMCH( 'Precision' )
177 SMLNUM = SAFMIN / EPS
178 BIGNUM = ONE / SMLNUM
179 RMIN = SQRT( SMLNUM )
180 RMAX = SQRT( BIGNUM )
181 *
182 * Scale matrix to allowable range, if necessary.
183 *
184 ISCALE = 0
185 TNRM = DLANST( 'M', N, D, E )
186 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
187 ISCALE = 1
188 SIGMA = RMIN / TNRM
189 ELSE IF( TNRM.GT.RMAX ) THEN
190 ISCALE = 1
191 SIGMA = RMAX / TNRM
192 END IF
193 IF( ISCALE.EQ.1 ) THEN
194 CALL DSCAL( N, SIGMA, D, 1 )
195 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
196 END IF
197 *
198 * For eigenvalues only, call DSTERF. For eigenvalues and
199 * eigenvectors, call DSTEDC.
200 *
201 IF( .NOT.WANTZ ) THEN
202 CALL DSTERF( N, D, E, INFO )
203 ELSE
204 CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK,
205 $ INFO )
206 END IF
207 *
208 * If matrix was scaled, then rescale eigenvalues appropriately.
209 *
210 IF( ISCALE.EQ.1 )
211 $ CALL DSCAL( N, ONE / SIGMA, D, 1 )
212 *
213 WORK( 1 ) = LWMIN
214 IWORK( 1 ) = LIWMIN
215 *
216 RETURN
217 *
218 * End of DSTEVD
219 *
220 END
2 $ 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
11 INTEGER INFO, LDZ, LIWORK, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSTEVD computes all eigenvalues and, optionally, eigenvectors of a
22 * real symmetric tridiagonal matrix. If eigenvectors are desired, it
23 * 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 * N (input) INTEGER
40 * The order of the matrix. N >= 0.
41 *
42 * D (input/output) DOUBLE PRECISION array, dimension (N)
43 * On entry, the n diagonal elements of the tridiagonal matrix
44 * A.
45 * On exit, if INFO = 0, the eigenvalues in ascending order.
46 *
47 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
48 * On entry, the (n-1) subdiagonal elements of the tridiagonal
49 * matrix A, stored in elements 1 to N-1 of E.
50 * On exit, the contents of E are destroyed.
51 *
52 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
53 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
54 * eigenvectors of the matrix A, with the i-th column of Z
55 * holding the eigenvector associated with D(i).
56 * If JOBZ = 'N', then Z is not referenced.
57 *
58 * LDZ (input) INTEGER
59 * The leading dimension of the array Z. LDZ >= 1, and if
60 * JOBZ = 'V', LDZ >= max(1,N).
61 *
62 * WORK (workspace/output) DOUBLE PRECISION array,
63 * dimension (LWORK)
64 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
65 *
66 * LWORK (input) INTEGER
67 * The dimension of the array WORK.
68 * If JOBZ = 'N' or N <= 1 then LWORK must be at least 1.
69 * If JOBZ = 'V' and N > 1 then LWORK must be at least
70 * ( 1 + 4*N + N**2 ).
71 *
72 * If LWORK = -1, then a workspace query is assumed; the routine
73 * only calculates the optimal sizes of the WORK and IWORK
74 * arrays, returns these values as the first entries of the WORK
75 * and IWORK arrays, and no error message related to LWORK or
76 * LIWORK is issued by XERBLA.
77 *
78 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
79 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
80 *
81 * LIWORK (input) INTEGER
82 * The dimension of the array IWORK.
83 * If JOBZ = 'N' or N <= 1 then LIWORK must be at least 1.
84 * If JOBZ = 'V' and N > 1 then LIWORK must be at least 3+5*N.
85 *
86 * If LIWORK = -1, then a workspace query is assumed; the
87 * routine only calculates the optimal sizes of the WORK and
88 * IWORK arrays, returns these values as the first entries of
89 * the WORK and IWORK arrays, and no error message related to
90 * LWORK or LIWORK is issued by XERBLA.
91 *
92 * INFO (output) INTEGER
93 * = 0: successful exit
94 * < 0: if INFO = -i, the i-th argument had an illegal value
95 * > 0: if INFO = i, the algorithm failed to converge; i
96 * off-diagonal elements of E did not converge to zero.
97 *
98 * =====================================================================
99 *
100 * .. Parameters ..
101 DOUBLE PRECISION ZERO, ONE
102 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
103 * ..
104 * .. Local Scalars ..
105 LOGICAL LQUERY, WANTZ
106 INTEGER ISCALE, LIWMIN, LWMIN
107 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
108 $ TNRM
109 * ..
110 * .. External Functions ..
111 LOGICAL LSAME
112 DOUBLE PRECISION DLAMCH, DLANST
113 EXTERNAL LSAME, DLAMCH, DLANST
114 * ..
115 * .. External Subroutines ..
116 EXTERNAL DSCAL, DSTEDC, DSTERF, XERBLA
117 * ..
118 * .. Intrinsic Functions ..
119 INTRINSIC SQRT
120 * ..
121 * .. Executable Statements ..
122 *
123 * Test the input parameters.
124 *
125 WANTZ = LSAME( JOBZ, 'V' )
126 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
127 *
128 INFO = 0
129 LIWMIN = 1
130 LWMIN = 1
131 IF( N.GT.1 .AND. WANTZ ) THEN
132 LWMIN = 1 + 4*N + N**2
133 LIWMIN = 3 + 5*N
134 END IF
135 *
136 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
137 INFO = -1
138 ELSE IF( N.LT.0 ) THEN
139 INFO = -2
140 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
141 INFO = -6
142 END IF
143 *
144 IF( INFO.EQ.0 ) THEN
145 WORK( 1 ) = LWMIN
146 IWORK( 1 ) = LIWMIN
147 *
148 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
149 INFO = -8
150 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
151 INFO = -10
152 END IF
153 END IF
154 *
155 IF( INFO.NE.0 ) THEN
156 CALL XERBLA( 'DSTEVD', -INFO )
157 RETURN
158 ELSE IF( LQUERY ) THEN
159 RETURN
160 END IF
161 *
162 * Quick return if possible
163 *
164 IF( N.EQ.0 )
165 $ RETURN
166 *
167 IF( N.EQ.1 ) THEN
168 IF( WANTZ )
169 $ Z( 1, 1 ) = ONE
170 RETURN
171 END IF
172 *
173 * Get machine constants.
174 *
175 SAFMIN = DLAMCH( 'Safe minimum' )
176 EPS = DLAMCH( 'Precision' )
177 SMLNUM = SAFMIN / EPS
178 BIGNUM = ONE / SMLNUM
179 RMIN = SQRT( SMLNUM )
180 RMAX = SQRT( BIGNUM )
181 *
182 * Scale matrix to allowable range, if necessary.
183 *
184 ISCALE = 0
185 TNRM = DLANST( 'M', N, D, E )
186 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
187 ISCALE = 1
188 SIGMA = RMIN / TNRM
189 ELSE IF( TNRM.GT.RMAX ) THEN
190 ISCALE = 1
191 SIGMA = RMAX / TNRM
192 END IF
193 IF( ISCALE.EQ.1 ) THEN
194 CALL DSCAL( N, SIGMA, D, 1 )
195 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
196 END IF
197 *
198 * For eigenvalues only, call DSTERF. For eigenvalues and
199 * eigenvectors, call DSTEDC.
200 *
201 IF( .NOT.WANTZ ) THEN
202 CALL DSTERF( N, D, E, INFO )
203 ELSE
204 CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK,
205 $ INFO )
206 END IF
207 *
208 * If matrix was scaled, then rescale eigenvalues appropriately.
209 *
210 IF( ISCALE.EQ.1 )
211 $ CALL DSCAL( N, ONE / SIGMA, D, 1 )
212 *
213 WORK( 1 ) = LWMIN
214 IWORK( 1 ) = LIWMIN
215 *
216 RETURN
217 *
218 * End of DSTEVD
219 *
220 END