Annotation of rpl/lapack/lapack/dstedc.f, revision 1.10
1.8 bertrand 1: *> \brief \b DSTEBZ
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,
109: *> dimension (LWORK)
110: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
111: *> \endverbatim
112: *>
113: *> \param[in] LWORK
114: *> \verbatim
115: *> LWORK is INTEGER
116: *> The dimension of the array WORK.
117: *> If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
118: *> If COMPZ = 'V' and N > 1 then LWORK must be at least
119: *> ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
120: *> where lg( N ) = smallest integer k such
121: *> that 2**k >= N.
122: *> If COMPZ = 'I' and N > 1 then LWORK must be at least
123: *> ( 1 + 4*N + N**2 ).
124: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
125: *> equal to the minimum divide size, usually 25, then LWORK need
126: *> only be max(1,2*(N-1)).
127: *>
128: *> If LWORK = -1, then a workspace query is assumed; the routine
129: *> only calculates the optimal size of the WORK array, returns
130: *> this value as the first entry of the WORK array, and no error
131: *> message related to LWORK is issued by XERBLA.
132: *> \endverbatim
133: *>
134: *> \param[out] IWORK
135: *> \verbatim
136: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
137: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
138: *> \endverbatim
139: *>
140: *> \param[in] LIWORK
141: *> \verbatim
142: *> LIWORK is INTEGER
143: *> The dimension of the array IWORK.
144: *> If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
145: *> If COMPZ = 'V' and N > 1 then LIWORK must be at least
146: *> ( 6 + 6*N + 5*N*lg N ).
147: *> If COMPZ = 'I' and N > 1 then LIWORK must be at least
148: *> ( 3 + 5*N ).
149: *> Note that for COMPZ = 'I' or 'V', then if N is less than or
150: *> equal to the minimum divide size, usually 25, then LIWORK
151: *> need only be 1.
152: *>
153: *> If LIWORK = -1, then a workspace query is assumed; the
154: *> routine only calculates the optimal size of the IWORK array,
155: *> returns this value as the first entry of the IWORK array, and
156: *> no error message related to LIWORK is issued by XERBLA.
157: *> \endverbatim
158: *>
159: *> \param[out] INFO
160: *> \verbatim
161: *> INFO is INTEGER
162: *> = 0: successful exit.
163: *> < 0: if INFO = -i, the i-th argument had an illegal value.
164: *> > 0: The algorithm failed to compute an eigenvalue while
165: *> working on the submatrix lying in rows and columns
166: *> INFO/(N+1) through mod(INFO,N+1).
167: *> \endverbatim
168: *
169: * Authors:
170: * ========
171: *
172: *> \author Univ. of Tennessee
173: *> \author Univ. of California Berkeley
174: *> \author Univ. of Colorado Denver
175: *> \author NAG Ltd.
176: *
177: *> \date November 2011
178: *
179: *> \ingroup auxOTHERcomputational
180: *
181: *> \par Contributors:
182: * ==================
183: *>
184: *> Jeff Rutter, Computer Science Division, University of California
185: *> at Berkeley, USA \n
186: *> Modified by Francoise Tisseur, University of Tennessee
187: *>
188: * =====================================================================
1.1 bertrand 189: SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
190: $ LIWORK, INFO )
191: *
1.8 bertrand 192: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 193: * -- LAPACK is a software package provided by Univ. of Tennessee, --
194: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 bertrand 195: * November 2011
1.1 bertrand 196: *
197: * .. Scalar Arguments ..
198: CHARACTER COMPZ
199: INTEGER INFO, LDZ, LIWORK, LWORK, N
200: * ..
201: * .. Array Arguments ..
202: INTEGER IWORK( * )
203: DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
204: * ..
205: *
206: * =====================================================================
207: *
208: * .. Parameters ..
209: DOUBLE PRECISION ZERO, ONE, TWO
210: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
211: * ..
212: * .. Local Scalars ..
213: LOGICAL LQUERY
214: INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
215: $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
216: DOUBLE PRECISION EPS, ORGNRM, P, TINY
217: * ..
218: * .. External Functions ..
219: LOGICAL LSAME
220: INTEGER ILAENV
221: DOUBLE PRECISION DLAMCH, DLANST
222: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
223: * ..
224: * .. External Subroutines ..
225: EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
226: $ DSTEQR, DSTERF, DSWAP, XERBLA
227: * ..
228: * .. Intrinsic Functions ..
229: INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
230: * ..
231: * .. Executable Statements ..
232: *
233: * Test the input parameters.
234: *
235: INFO = 0
236: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
237: *
238: IF( LSAME( COMPZ, 'N' ) ) THEN
239: ICOMPZ = 0
240: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
241: ICOMPZ = 1
242: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
243: ICOMPZ = 2
244: ELSE
245: ICOMPZ = -1
246: END IF
247: IF( ICOMPZ.LT.0 ) THEN
248: INFO = -1
249: ELSE IF( N.LT.0 ) THEN
250: INFO = -2
251: ELSE IF( ( LDZ.LT.1 ) .OR.
252: $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
253: INFO = -6
254: END IF
255: *
256: IF( INFO.EQ.0 ) THEN
257: *
258: * Compute the workspace requirements
259: *
260: SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
261: IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
262: LIWMIN = 1
263: LWMIN = 1
264: ELSE IF( N.LE.SMLSIZ ) THEN
265: LIWMIN = 1
266: LWMIN = 2*( N - 1 )
267: ELSE
268: LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
269: IF( 2**LGN.LT.N )
270: $ LGN = LGN + 1
271: IF( 2**LGN.LT.N )
272: $ LGN = LGN + 1
273: IF( ICOMPZ.EQ.1 ) THEN
1.8 bertrand 274: LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
1.1 bertrand 275: LIWMIN = 6 + 6*N + 5*N*LGN
276: ELSE IF( ICOMPZ.EQ.2 ) THEN
277: LWMIN = 1 + 4*N + N**2
278: LIWMIN = 3 + 5*N
279: END IF
280: END IF
281: WORK( 1 ) = LWMIN
282: IWORK( 1 ) = LIWMIN
283: *
284: IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
285: INFO = -8
286: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
287: INFO = -10
288: END IF
289: END IF
290: *
291: IF( INFO.NE.0 ) THEN
292: CALL XERBLA( 'DSTEDC', -INFO )
293: RETURN
294: ELSE IF (LQUERY) THEN
295: RETURN
296: END IF
297: *
298: * Quick return if possible
299: *
300: IF( N.EQ.0 )
301: $ RETURN
302: IF( N.EQ.1 ) THEN
303: IF( ICOMPZ.NE.0 )
304: $ Z( 1, 1 ) = ONE
305: RETURN
306: END IF
307: *
308: * If the following conditional clause is removed, then the routine
309: * will use the Divide and Conquer routine to compute only the
310: * eigenvalues, which requires (3N + 3N**2) real workspace and
311: * (2 + 5N + 2N lg(N)) integer workspace.
312: * Since on many architectures DSTERF is much faster than any other
313: * algorithm for finding eigenvalues only, it is used here
314: * as the default. If the conditional clause is removed, then
315: * information on the size of workspace needs to be changed.
316: *
317: * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
318: *
319: IF( ICOMPZ.EQ.0 ) THEN
320: CALL DSTERF( N, D, E, INFO )
321: GO TO 50
322: END IF
323: *
324: * If N is smaller than the minimum divide size (SMLSIZ+1), then
325: * solve the problem with another solver.
326: *
327: IF( N.LE.SMLSIZ ) THEN
328: *
329: CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
330: *
331: ELSE
332: *
333: * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
334: * use.
335: *
336: IF( ICOMPZ.EQ.1 ) THEN
337: STOREZ = 1 + N*N
338: ELSE
339: STOREZ = 1
340: END IF
341: *
342: IF( ICOMPZ.EQ.2 ) THEN
343: CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
344: END IF
345: *
346: * Scale.
347: *
348: ORGNRM = DLANST( 'M', N, D, E )
349: IF( ORGNRM.EQ.ZERO )
350: $ GO TO 50
351: *
352: EPS = DLAMCH( 'Epsilon' )
353: *
354: START = 1
355: *
356: * while ( START <= N )
357: *
358: 10 CONTINUE
359: IF( START.LE.N ) THEN
360: *
361: * Let FINISH be the position of the next subdiagonal entry
362: * such that E( FINISH ) <= TINY or FINISH = N if no such
363: * subdiagonal exists. The matrix identified by the elements
364: * between START and FINISH constitutes an independent
365: * sub-problem.
366: *
367: FINISH = START
368: 20 CONTINUE
369: IF( FINISH.LT.N ) THEN
370: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
371: $ SQRT( ABS( D( FINISH+1 ) ) )
372: IF( ABS( E( FINISH ) ).GT.TINY ) THEN
373: FINISH = FINISH + 1
374: GO TO 20
375: END IF
376: END IF
377: *
378: * (Sub) Problem determined. Compute its size and solve it.
379: *
380: M = FINISH - START + 1
381: IF( M.EQ.1 ) THEN
382: START = FINISH + 1
383: GO TO 10
384: END IF
385: IF( M.GT.SMLSIZ ) THEN
386: *
387: * Scale.
388: *
389: ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
390: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
391: $ INFO )
392: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
393: $ M-1, INFO )
394: *
395: IF( ICOMPZ.EQ.1 ) THEN
396: STRTRW = 1
397: ELSE
398: STRTRW = START
399: END IF
400: CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
401: $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
402: $ WORK( STOREZ ), IWORK, INFO )
403: IF( INFO.NE.0 ) THEN
404: INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
405: $ MOD( INFO, ( M+1 ) ) + START - 1
406: GO TO 50
407: END IF
408: *
409: * Scale back.
410: *
411: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
412: $ INFO )
413: *
414: ELSE
415: IF( ICOMPZ.EQ.1 ) THEN
416: *
417: * Since QR won't update a Z matrix which is larger than
418: * the length of D, we must solve the sub-problem in a
419: * workspace and then multiply back into Z.
420: *
421: CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
422: $ WORK( M*M+1 ), INFO )
423: CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
424: $ WORK( STOREZ ), N )
425: CALL DGEMM( 'N', 'N', N, M, M, ONE,
426: $ WORK( STOREZ ), N, WORK, M, ZERO,
427: $ Z( 1, START ), LDZ )
428: ELSE IF( ICOMPZ.EQ.2 ) THEN
429: CALL DSTEQR( 'I', M, D( START ), E( START ),
430: $ Z( START, START ), LDZ, WORK, INFO )
431: ELSE
432: CALL DSTERF( M, D( START ), E( START ), INFO )
433: END IF
434: IF( INFO.NE.0 ) THEN
435: INFO = START*( N+1 ) + FINISH
436: GO TO 50
437: END IF
438: END IF
439: *
440: START = FINISH + 1
441: GO TO 10
442: END IF
443: *
444: * endwhile
445: *
446: * If the problem split any number of times, then the eigenvalues
447: * will not be properly ordered. Here we permute the eigenvalues
448: * (and the associated eigenvectors) into ascending order.
449: *
450: IF( M.NE.N ) THEN
451: IF( ICOMPZ.EQ.0 ) THEN
452: *
453: * Use Quick Sort
454: *
455: CALL DLASRT( 'I', N, D, INFO )
456: *
457: ELSE
458: *
459: * Use Selection Sort to minimize swaps of eigenvectors
460: *
461: DO 40 II = 2, N
462: I = II - 1
463: K = I
464: P = D( I )
465: DO 30 J = II, N
466: IF( D( J ).LT.P ) THEN
467: K = J
468: P = D( J )
469: END IF
470: 30 CONTINUE
471: IF( K.NE.I ) THEN
472: D( K ) = D( I )
473: D( I ) = P
474: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
475: END IF
476: 40 CONTINUE
477: END IF
478: END IF
479: END IF
480: *
481: 50 CONTINUE
482: WORK( 1 ) = LWMIN
483: IWORK( 1 ) = LIWMIN
484: *
485: RETURN
486: *
487: * End of DSTEDC
488: *
489: 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