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