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