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