1 SUBROUTINE DSTEDC( COMPZ, 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 COMPZ
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 * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
22 * symmetric tridiagonal matrix using the divide and conquer method.
23 * The eigenvectors of a full or band real symmetric matrix can also be
24 * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
25 * matrix to tridiagonal form.
26 *
27 * This code makes very mild assumptions about floating point
28 * arithmetic. It will work on machines with a guard digit in
29 * add/subtract, or on those binary machines without guard digits
30 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
31 * It could conceivably fail on hexadecimal or decimal machines
32 * without guard digits, but we know of none. See DLAED3 for details.
33 *
34 * Arguments
35 * =========
36 *
37 * COMPZ (input) CHARACTER*1
38 * = 'N': Compute eigenvalues only.
39 * = 'I': Compute eigenvectors of tridiagonal matrix also.
40 * = 'V': Compute eigenvectors of original dense symmetric
41 * matrix also. On entry, Z contains the orthogonal
42 * matrix used to reduce the original matrix to
43 * tridiagonal form.
44 *
45 * N (input) INTEGER
46 * The dimension of the symmetric tridiagonal matrix. N >= 0.
47 *
48 * D (input/output) DOUBLE PRECISION array, dimension (N)
49 * On entry, the diagonal elements of the tridiagonal matrix.
50 * On exit, if INFO = 0, the eigenvalues in ascending order.
51 *
52 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
53 * On entry, the subdiagonal elements of the tridiagonal matrix.
54 * On exit, E has been destroyed.
55 *
56 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
57 * On entry, if COMPZ = 'V', then Z contains the orthogonal
58 * matrix used in the reduction to tridiagonal form.
59 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
60 * orthonormal eigenvectors of the original symmetric matrix,
61 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
62 * of the symmetric tridiagonal matrix.
63 * If COMPZ = 'N', then Z is not referenced.
64 *
65 * LDZ (input) INTEGER
66 * The leading dimension of the array Z. LDZ >= 1.
67 * If eigenvectors are desired, then LDZ >= max(1,N).
68 *
69 * WORK (workspace/output) DOUBLE PRECISION array,
70 * dimension (LWORK)
71 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
72 *
73 * LWORK (input) INTEGER
74 * The dimension of the array WORK.
75 * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
76 * If COMPZ = 'V' and N > 1 then LWORK must be at least
77 * ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
78 * where lg( N ) = smallest integer k such
79 * that 2**k >= N.
80 * If COMPZ = 'I' and N > 1 then LWORK must be at least
81 * ( 1 + 4*N + N**2 ).
82 * Note that for COMPZ = 'I' or 'V', then if N is less than or
83 * equal to the minimum divide size, usually 25, then LWORK need
84 * only be max(1,2*(N-1)).
85 *
86 * If LWORK = -1, then a workspace query is assumed; the routine
87 * only calculates the optimal size of the WORK array, returns
88 * this value as the first entry of the WORK array, and no error
89 * message related to LWORK is issued by XERBLA.
90 *
91 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
92 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
93 *
94 * LIWORK (input) INTEGER
95 * The dimension of the array IWORK.
96 * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
97 * If COMPZ = 'V' and N > 1 then LIWORK must be at least
98 * ( 6 + 6*N + 5*N*lg N ).
99 * If COMPZ = 'I' and N > 1 then LIWORK must be at least
100 * ( 3 + 5*N ).
101 * Note that for COMPZ = 'I' or 'V', then if N is less than or
102 * equal to the minimum divide size, usually 25, then LIWORK
103 * need only be 1.
104 *
105 * If LIWORK = -1, then a workspace query is assumed; the
106 * routine only calculates the optimal size of the IWORK array,
107 * returns this value as the first entry of the IWORK array, and
108 * no error message related to LIWORK is issued by XERBLA.
109 *
110 * INFO (output) INTEGER
111 * = 0: successful exit.
112 * < 0: if INFO = -i, the i-th argument had an illegal value.
113 * > 0: The algorithm failed to compute an eigenvalue while
114 * working on the submatrix lying in rows and columns
115 * INFO/(N+1) through mod(INFO,N+1).
116 *
117 * Further Details
118 * ===============
119 *
120 * Based on contributions by
121 * Jeff Rutter, Computer Science Division, University of California
122 * at Berkeley, USA
123 * Modified by Francoise Tisseur, University of Tennessee.
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE, TWO
129 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
130 * ..
131 * .. Local Scalars ..
132 LOGICAL LQUERY
133 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
134 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
135 DOUBLE PRECISION EPS, ORGNRM, P, TINY
136 * ..
137 * .. External Functions ..
138 LOGICAL LSAME
139 INTEGER ILAENV
140 DOUBLE PRECISION DLAMCH, DLANST
141 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
142 * ..
143 * .. External Subroutines ..
144 EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
145 $ DSTEQR, DSTERF, DSWAP, XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
149 * ..
150 * .. Executable Statements ..
151 *
152 * Test the input parameters.
153 *
154 INFO = 0
155 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
156 *
157 IF( LSAME( COMPZ, 'N' ) ) THEN
158 ICOMPZ = 0
159 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
160 ICOMPZ = 1
161 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
162 ICOMPZ = 2
163 ELSE
164 ICOMPZ = -1
165 END IF
166 IF( ICOMPZ.LT.0 ) THEN
167 INFO = -1
168 ELSE IF( N.LT.0 ) THEN
169 INFO = -2
170 ELSE IF( ( LDZ.LT.1 ) .OR.
171 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
172 INFO = -6
173 END IF
174 *
175 IF( INFO.EQ.0 ) THEN
176 *
177 * Compute the workspace requirements
178 *
179 SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
180 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
181 LIWMIN = 1
182 LWMIN = 1
183 ELSE IF( N.LE.SMLSIZ ) THEN
184 LIWMIN = 1
185 LWMIN = 2*( N - 1 )
186 ELSE
187 LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
188 IF( 2**LGN.LT.N )
189 $ LGN = LGN + 1
190 IF( 2**LGN.LT.N )
191 $ LGN = LGN + 1
192 IF( ICOMPZ.EQ.1 ) THEN
193 LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
194 LIWMIN = 6 + 6*N + 5*N*LGN
195 ELSE IF( ICOMPZ.EQ.2 ) THEN
196 LWMIN = 1 + 4*N + N**2
197 LIWMIN = 3 + 5*N
198 END IF
199 END IF
200 WORK( 1 ) = LWMIN
201 IWORK( 1 ) = LIWMIN
202 *
203 IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
204 INFO = -8
205 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
206 INFO = -10
207 END IF
208 END IF
209 *
210 IF( INFO.NE.0 ) THEN
211 CALL XERBLA( 'DSTEDC', -INFO )
212 RETURN
213 ELSE IF (LQUERY) THEN
214 RETURN
215 END IF
216 *
217 * Quick return if possible
218 *
219 IF( N.EQ.0 )
220 $ RETURN
221 IF( N.EQ.1 ) THEN
222 IF( ICOMPZ.NE.0 )
223 $ Z( 1, 1 ) = ONE
224 RETURN
225 END IF
226 *
227 * If the following conditional clause is removed, then the routine
228 * will use the Divide and Conquer routine to compute only the
229 * eigenvalues, which requires (3N + 3N**2) real workspace and
230 * (2 + 5N + 2N lg(N)) integer workspace.
231 * Since on many architectures DSTERF is much faster than any other
232 * algorithm for finding eigenvalues only, it is used here
233 * as the default. If the conditional clause is removed, then
234 * information on the size of workspace needs to be changed.
235 *
236 * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
237 *
238 IF( ICOMPZ.EQ.0 ) THEN
239 CALL DSTERF( N, D, E, INFO )
240 GO TO 50
241 END IF
242 *
243 * If N is smaller than the minimum divide size (SMLSIZ+1), then
244 * solve the problem with another solver.
245 *
246 IF( N.LE.SMLSIZ ) THEN
247 *
248 CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
249 *
250 ELSE
251 *
252 * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
253 * use.
254 *
255 IF( ICOMPZ.EQ.1 ) THEN
256 STOREZ = 1 + N*N
257 ELSE
258 STOREZ = 1
259 END IF
260 *
261 IF( ICOMPZ.EQ.2 ) THEN
262 CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
263 END IF
264 *
265 * Scale.
266 *
267 ORGNRM = DLANST( 'M', N, D, E )
268 IF( ORGNRM.EQ.ZERO )
269 $ GO TO 50
270 *
271 EPS = DLAMCH( 'Epsilon' )
272 *
273 START = 1
274 *
275 * while ( START <= N )
276 *
277 10 CONTINUE
278 IF( START.LE.N ) THEN
279 *
280 * Let FINISH be the position of the next subdiagonal entry
281 * such that E( FINISH ) <= TINY or FINISH = N if no such
282 * subdiagonal exists. The matrix identified by the elements
283 * between START and FINISH constitutes an independent
284 * sub-problem.
285 *
286 FINISH = START
287 20 CONTINUE
288 IF( FINISH.LT.N ) THEN
289 TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
290 $ SQRT( ABS( D( FINISH+1 ) ) )
291 IF( ABS( E( FINISH ) ).GT.TINY ) THEN
292 FINISH = FINISH + 1
293 GO TO 20
294 END IF
295 END IF
296 *
297 * (Sub) Problem determined. Compute its size and solve it.
298 *
299 M = FINISH - START + 1
300 IF( M.EQ.1 ) THEN
301 START = FINISH + 1
302 GO TO 10
303 END IF
304 IF( M.GT.SMLSIZ ) THEN
305 *
306 * Scale.
307 *
308 ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
309 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
310 $ INFO )
311 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
312 $ M-1, INFO )
313 *
314 IF( ICOMPZ.EQ.1 ) THEN
315 STRTRW = 1
316 ELSE
317 STRTRW = START
318 END IF
319 CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
320 $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
321 $ WORK( STOREZ ), IWORK, INFO )
322 IF( INFO.NE.0 ) THEN
323 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
324 $ MOD( INFO, ( M+1 ) ) + START - 1
325 GO TO 50
326 END IF
327 *
328 * Scale back.
329 *
330 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
331 $ INFO )
332 *
333 ELSE
334 IF( ICOMPZ.EQ.1 ) THEN
335 *
336 * Since QR won't update a Z matrix which is larger than
337 * the length of D, we must solve the sub-problem in a
338 * workspace and then multiply back into Z.
339 *
340 CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
341 $ WORK( M*M+1 ), INFO )
342 CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
343 $ WORK( STOREZ ), N )
344 CALL DGEMM( 'N', 'N', N, M, M, ONE,
345 $ WORK( STOREZ ), N, WORK, M, ZERO,
346 $ Z( 1, START ), LDZ )
347 ELSE IF( ICOMPZ.EQ.2 ) THEN
348 CALL DSTEQR( 'I', M, D( START ), E( START ),
349 $ Z( START, START ), LDZ, WORK, INFO )
350 ELSE
351 CALL DSTERF( M, D( START ), E( START ), INFO )
352 END IF
353 IF( INFO.NE.0 ) THEN
354 INFO = START*( N+1 ) + FINISH
355 GO TO 50
356 END IF
357 END IF
358 *
359 START = FINISH + 1
360 GO TO 10
361 END IF
362 *
363 * endwhile
364 *
365 * If the problem split any number of times, then the eigenvalues
366 * will not be properly ordered. Here we permute the eigenvalues
367 * (and the associated eigenvectors) into ascending order.
368 *
369 IF( M.NE.N ) THEN
370 IF( ICOMPZ.EQ.0 ) THEN
371 *
372 * Use Quick Sort
373 *
374 CALL DLASRT( 'I', N, D, INFO )
375 *
376 ELSE
377 *
378 * Use Selection Sort to minimize swaps of eigenvectors
379 *
380 DO 40 II = 2, N
381 I = II - 1
382 K = I
383 P = D( I )
384 DO 30 J = II, N
385 IF( D( J ).LT.P ) THEN
386 K = J
387 P = D( J )
388 END IF
389 30 CONTINUE
390 IF( K.NE.I ) THEN
391 D( K ) = D( I )
392 D( I ) = P
393 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
394 END IF
395 40 CONTINUE
396 END IF
397 END IF
398 END IF
399 *
400 50 CONTINUE
401 WORK( 1 ) = LWMIN
402 IWORK( 1 ) = LIWMIN
403 *
404 RETURN
405 *
406 * End of DSTEDC
407 *
408 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 COMPZ
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 * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
22 * symmetric tridiagonal matrix using the divide and conquer method.
23 * The eigenvectors of a full or band real symmetric matrix can also be
24 * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
25 * matrix to tridiagonal form.
26 *
27 * This code makes very mild assumptions about floating point
28 * arithmetic. It will work on machines with a guard digit in
29 * add/subtract, or on those binary machines without guard digits
30 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
31 * It could conceivably fail on hexadecimal or decimal machines
32 * without guard digits, but we know of none. See DLAED3 for details.
33 *
34 * Arguments
35 * =========
36 *
37 * COMPZ (input) CHARACTER*1
38 * = 'N': Compute eigenvalues only.
39 * = 'I': Compute eigenvectors of tridiagonal matrix also.
40 * = 'V': Compute eigenvectors of original dense symmetric
41 * matrix also. On entry, Z contains the orthogonal
42 * matrix used to reduce the original matrix to
43 * tridiagonal form.
44 *
45 * N (input) INTEGER
46 * The dimension of the symmetric tridiagonal matrix. N >= 0.
47 *
48 * D (input/output) DOUBLE PRECISION array, dimension (N)
49 * On entry, the diagonal elements of the tridiagonal matrix.
50 * On exit, if INFO = 0, the eigenvalues in ascending order.
51 *
52 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
53 * On entry, the subdiagonal elements of the tridiagonal matrix.
54 * On exit, E has been destroyed.
55 *
56 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
57 * On entry, if COMPZ = 'V', then Z contains the orthogonal
58 * matrix used in the reduction to tridiagonal form.
59 * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
60 * orthonormal eigenvectors of the original symmetric matrix,
61 * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
62 * of the symmetric tridiagonal matrix.
63 * If COMPZ = 'N', then Z is not referenced.
64 *
65 * LDZ (input) INTEGER
66 * The leading dimension of the array Z. LDZ >= 1.
67 * If eigenvectors are desired, then LDZ >= max(1,N).
68 *
69 * WORK (workspace/output) DOUBLE PRECISION array,
70 * dimension (LWORK)
71 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
72 *
73 * LWORK (input) INTEGER
74 * The dimension of the array WORK.
75 * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
76 * If COMPZ = 'V' and N > 1 then LWORK must be at least
77 * ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
78 * where lg( N ) = smallest integer k such
79 * that 2**k >= N.
80 * If COMPZ = 'I' and N > 1 then LWORK must be at least
81 * ( 1 + 4*N + N**2 ).
82 * Note that for COMPZ = 'I' or 'V', then if N is less than or
83 * equal to the minimum divide size, usually 25, then LWORK need
84 * only be max(1,2*(N-1)).
85 *
86 * If LWORK = -1, then a workspace query is assumed; the routine
87 * only calculates the optimal size of the WORK array, returns
88 * this value as the first entry of the WORK array, and no error
89 * message related to LWORK is issued by XERBLA.
90 *
91 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
92 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
93 *
94 * LIWORK (input) INTEGER
95 * The dimension of the array IWORK.
96 * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
97 * If COMPZ = 'V' and N > 1 then LIWORK must be at least
98 * ( 6 + 6*N + 5*N*lg N ).
99 * If COMPZ = 'I' and N > 1 then LIWORK must be at least
100 * ( 3 + 5*N ).
101 * Note that for COMPZ = 'I' or 'V', then if N is less than or
102 * equal to the minimum divide size, usually 25, then LIWORK
103 * need only be 1.
104 *
105 * If LIWORK = -1, then a workspace query is assumed; the
106 * routine only calculates the optimal size of the IWORK array,
107 * returns this value as the first entry of the IWORK array, and
108 * no error message related to LIWORK is issued by XERBLA.
109 *
110 * INFO (output) INTEGER
111 * = 0: successful exit.
112 * < 0: if INFO = -i, the i-th argument had an illegal value.
113 * > 0: The algorithm failed to compute an eigenvalue while
114 * working on the submatrix lying in rows and columns
115 * INFO/(N+1) through mod(INFO,N+1).
116 *
117 * Further Details
118 * ===============
119 *
120 * Based on contributions by
121 * Jeff Rutter, Computer Science Division, University of California
122 * at Berkeley, USA
123 * Modified by Francoise Tisseur, University of Tennessee.
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 DOUBLE PRECISION ZERO, ONE, TWO
129 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
130 * ..
131 * .. Local Scalars ..
132 LOGICAL LQUERY
133 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
134 $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
135 DOUBLE PRECISION EPS, ORGNRM, P, TINY
136 * ..
137 * .. External Functions ..
138 LOGICAL LSAME
139 INTEGER ILAENV
140 DOUBLE PRECISION DLAMCH, DLANST
141 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
142 * ..
143 * .. External Subroutines ..
144 EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
145 $ DSTEQR, DSTERF, DSWAP, XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
149 * ..
150 * .. Executable Statements ..
151 *
152 * Test the input parameters.
153 *
154 INFO = 0
155 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
156 *
157 IF( LSAME( COMPZ, 'N' ) ) THEN
158 ICOMPZ = 0
159 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
160 ICOMPZ = 1
161 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
162 ICOMPZ = 2
163 ELSE
164 ICOMPZ = -1
165 END IF
166 IF( ICOMPZ.LT.0 ) THEN
167 INFO = -1
168 ELSE IF( N.LT.0 ) THEN
169 INFO = -2
170 ELSE IF( ( LDZ.LT.1 ) .OR.
171 $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
172 INFO = -6
173 END IF
174 *
175 IF( INFO.EQ.0 ) THEN
176 *
177 * Compute the workspace requirements
178 *
179 SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
180 IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
181 LIWMIN = 1
182 LWMIN = 1
183 ELSE IF( N.LE.SMLSIZ ) THEN
184 LIWMIN = 1
185 LWMIN = 2*( N - 1 )
186 ELSE
187 LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
188 IF( 2**LGN.LT.N )
189 $ LGN = LGN + 1
190 IF( 2**LGN.LT.N )
191 $ LGN = LGN + 1
192 IF( ICOMPZ.EQ.1 ) THEN
193 LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
194 LIWMIN = 6 + 6*N + 5*N*LGN
195 ELSE IF( ICOMPZ.EQ.2 ) THEN
196 LWMIN = 1 + 4*N + N**2
197 LIWMIN = 3 + 5*N
198 END IF
199 END IF
200 WORK( 1 ) = LWMIN
201 IWORK( 1 ) = LIWMIN
202 *
203 IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
204 INFO = -8
205 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
206 INFO = -10
207 END IF
208 END IF
209 *
210 IF( INFO.NE.0 ) THEN
211 CALL XERBLA( 'DSTEDC', -INFO )
212 RETURN
213 ELSE IF (LQUERY) THEN
214 RETURN
215 END IF
216 *
217 * Quick return if possible
218 *
219 IF( N.EQ.0 )
220 $ RETURN
221 IF( N.EQ.1 ) THEN
222 IF( ICOMPZ.NE.0 )
223 $ Z( 1, 1 ) = ONE
224 RETURN
225 END IF
226 *
227 * If the following conditional clause is removed, then the routine
228 * will use the Divide and Conquer routine to compute only the
229 * eigenvalues, which requires (3N + 3N**2) real workspace and
230 * (2 + 5N + 2N lg(N)) integer workspace.
231 * Since on many architectures DSTERF is much faster than any other
232 * algorithm for finding eigenvalues only, it is used here
233 * as the default. If the conditional clause is removed, then
234 * information on the size of workspace needs to be changed.
235 *
236 * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
237 *
238 IF( ICOMPZ.EQ.0 ) THEN
239 CALL DSTERF( N, D, E, INFO )
240 GO TO 50
241 END IF
242 *
243 * If N is smaller than the minimum divide size (SMLSIZ+1), then
244 * solve the problem with another solver.
245 *
246 IF( N.LE.SMLSIZ ) THEN
247 *
248 CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
249 *
250 ELSE
251 *
252 * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
253 * use.
254 *
255 IF( ICOMPZ.EQ.1 ) THEN
256 STOREZ = 1 + N*N
257 ELSE
258 STOREZ = 1
259 END IF
260 *
261 IF( ICOMPZ.EQ.2 ) THEN
262 CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
263 END IF
264 *
265 * Scale.
266 *
267 ORGNRM = DLANST( 'M', N, D, E )
268 IF( ORGNRM.EQ.ZERO )
269 $ GO TO 50
270 *
271 EPS = DLAMCH( 'Epsilon' )
272 *
273 START = 1
274 *
275 * while ( START <= N )
276 *
277 10 CONTINUE
278 IF( START.LE.N ) THEN
279 *
280 * Let FINISH be the position of the next subdiagonal entry
281 * such that E( FINISH ) <= TINY or FINISH = N if no such
282 * subdiagonal exists. The matrix identified by the elements
283 * between START and FINISH constitutes an independent
284 * sub-problem.
285 *
286 FINISH = START
287 20 CONTINUE
288 IF( FINISH.LT.N ) THEN
289 TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
290 $ SQRT( ABS( D( FINISH+1 ) ) )
291 IF( ABS( E( FINISH ) ).GT.TINY ) THEN
292 FINISH = FINISH + 1
293 GO TO 20
294 END IF
295 END IF
296 *
297 * (Sub) Problem determined. Compute its size and solve it.
298 *
299 M = FINISH - START + 1
300 IF( M.EQ.1 ) THEN
301 START = FINISH + 1
302 GO TO 10
303 END IF
304 IF( M.GT.SMLSIZ ) THEN
305 *
306 * Scale.
307 *
308 ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
309 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
310 $ INFO )
311 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
312 $ M-1, INFO )
313 *
314 IF( ICOMPZ.EQ.1 ) THEN
315 STRTRW = 1
316 ELSE
317 STRTRW = START
318 END IF
319 CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
320 $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
321 $ WORK( STOREZ ), IWORK, INFO )
322 IF( INFO.NE.0 ) THEN
323 INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
324 $ MOD( INFO, ( M+1 ) ) + START - 1
325 GO TO 50
326 END IF
327 *
328 * Scale back.
329 *
330 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
331 $ INFO )
332 *
333 ELSE
334 IF( ICOMPZ.EQ.1 ) THEN
335 *
336 * Since QR won't update a Z matrix which is larger than
337 * the length of D, we must solve the sub-problem in a
338 * workspace and then multiply back into Z.
339 *
340 CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
341 $ WORK( M*M+1 ), INFO )
342 CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
343 $ WORK( STOREZ ), N )
344 CALL DGEMM( 'N', 'N', N, M, M, ONE,
345 $ WORK( STOREZ ), N, WORK, M, ZERO,
346 $ Z( 1, START ), LDZ )
347 ELSE IF( ICOMPZ.EQ.2 ) THEN
348 CALL DSTEQR( 'I', M, D( START ), E( START ),
349 $ Z( START, START ), LDZ, WORK, INFO )
350 ELSE
351 CALL DSTERF( M, D( START ), E( START ), INFO )
352 END IF
353 IF( INFO.NE.0 ) THEN
354 INFO = START*( N+1 ) + FINISH
355 GO TO 50
356 END IF
357 END IF
358 *
359 START = FINISH + 1
360 GO TO 10
361 END IF
362 *
363 * endwhile
364 *
365 * If the problem split any number of times, then the eigenvalues
366 * will not be properly ordered. Here we permute the eigenvalues
367 * (and the associated eigenvectors) into ascending order.
368 *
369 IF( M.NE.N ) THEN
370 IF( ICOMPZ.EQ.0 ) THEN
371 *
372 * Use Quick Sort
373 *
374 CALL DLASRT( 'I', N, D, INFO )
375 *
376 ELSE
377 *
378 * Use Selection Sort to minimize swaps of eigenvectors
379 *
380 DO 40 II = 2, N
381 I = II - 1
382 K = I
383 P = D( I )
384 DO 30 J = II, N
385 IF( D( J ).LT.P ) THEN
386 K = J
387 P = D( J )
388 END IF
389 30 CONTINUE
390 IF( K.NE.I ) THEN
391 D( K ) = D( I )
392 D( I ) = P
393 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
394 END IF
395 40 CONTINUE
396 END IF
397 END IF
398 END IF
399 *
400 50 CONTINUE
401 WORK( 1 ) = LWMIN
402 IWORK( 1 ) = LIWMIN
403 *
404 RETURN
405 *
406 * End of DSTEDC
407 *
408 END