Annotation of rpl/lapack/lapack/zstedc.f, revision 1.7
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zstedc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack