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