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