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: *> \date June 2017
202: *
203: *> \ingroup complex16OTHERcomputational
204: *
205: *> \par Contributors:
206: * ==================
207: *>
208: *> Jeff Rutter, Computer Science Division, University of California
209: *> at Berkeley, USA
210: *
211: * =====================================================================
212: SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
213: $ LRWORK, IWORK, LIWORK, INFO )
214: *
215: * -- LAPACK computational routine (version 3.7.1) --
216: * -- LAPACK is a software package provided by Univ. of Tennessee, --
217: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
218: * June 2017
219: *
220: * .. Scalar Arguments ..
221: CHARACTER COMPZ
222: INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
223: * ..
224: * .. Array Arguments ..
225: INTEGER IWORK( * )
226: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
227: COMPLEX*16 WORK( * ), Z( LDZ, * )
228: * ..
229: *
230: * =====================================================================
231: *
232: * .. Parameters ..
233: DOUBLE PRECISION ZERO, ONE, TWO
234: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
235: * ..
236: * .. Local Scalars ..
237: LOGICAL LQUERY
238: INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
239: $ LRWMIN, LWMIN, M, SMLSIZ, START
240: DOUBLE PRECISION EPS, ORGNRM, P, TINY
241: * ..
242: * .. External Functions ..
243: LOGICAL LSAME
244: INTEGER ILAENV
245: DOUBLE PRECISION DLAMCH, DLANST
246: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
247: * ..
248: * .. External Subroutines ..
249: EXTERNAL DLASCL, DLASET, DSTEDC, DSTEQR, DSTERF, XERBLA,
250: $ ZLACPY, ZLACRM, ZLAED0, ZSTEQR, ZSWAP
251: * ..
252: * .. Intrinsic Functions ..
253: INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
254: * ..
255: * .. Executable Statements ..
256: *
257: * Test the input parameters.
258: *
259: INFO = 0
260: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
261: *
262: IF( LSAME( COMPZ, 'N' ) ) THEN
263: ICOMPZ = 0
264: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
265: ICOMPZ = 1
266: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
267: ICOMPZ = 2
268: ELSE
269: ICOMPZ = -1
270: END IF
271: IF( ICOMPZ.LT.0 ) THEN
272: INFO = -1
273: ELSE IF( N.LT.0 ) THEN
274: INFO = -2
275: ELSE IF( ( LDZ.LT.1 ) .OR.
276: $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
277: INFO = -6
278: END IF
279: *
280: IF( INFO.EQ.0 ) THEN
281: *
282: * Compute the workspace requirements
283: *
284: SMLSIZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 )
285: IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
286: LWMIN = 1
287: LIWMIN = 1
288: LRWMIN = 1
289: ELSE IF( N.LE.SMLSIZ ) THEN
290: LWMIN = 1
291: LIWMIN = 1
292: LRWMIN = 2*( N - 1 )
293: ELSE IF( ICOMPZ.EQ.1 ) THEN
294: LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) )
295: IF( 2**LGN.LT.N )
296: $ LGN = LGN + 1
297: IF( 2**LGN.LT.N )
298: $ LGN = LGN + 1
299: LWMIN = N*N
300: LRWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
301: LIWMIN = 6 + 6*N + 5*N*LGN
302: ELSE IF( ICOMPZ.EQ.2 ) THEN
303: LWMIN = 1
304: LRWMIN = 1 + 4*N + 2*N**2
305: LIWMIN = 3 + 5*N
306: END IF
307: WORK( 1 ) = LWMIN
308: RWORK( 1 ) = LRWMIN
309: IWORK( 1 ) = LIWMIN
310: *
311: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
312: INFO = -8
313: ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
314: INFO = -10
315: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
316: INFO = -12
317: END IF
318: END IF
319: *
320: IF( INFO.NE.0 ) THEN
321: CALL XERBLA( 'ZSTEDC', -INFO )
322: RETURN
323: ELSE IF( LQUERY ) THEN
324: RETURN
325: END IF
326: *
327: * Quick return if possible
328: *
329: IF( N.EQ.0 )
330: $ RETURN
331: IF( N.EQ.1 ) THEN
332: IF( ICOMPZ.NE.0 )
333: $ Z( 1, 1 ) = ONE
334: RETURN
335: END IF
336: *
337: * If the following conditional clause is removed, then the routine
338: * will use the Divide and Conquer routine to compute only the
339: * eigenvalues, which requires (3N + 3N**2) real workspace and
340: * (2 + 5N + 2N lg(N)) integer workspace.
341: * Since on many architectures DSTERF is much faster than any other
342: * algorithm for finding eigenvalues only, it is used here
343: * as the default. If the conditional clause is removed, then
344: * information on the size of workspace needs to be changed.
345: *
346: * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
347: *
348: IF( ICOMPZ.EQ.0 ) THEN
349: CALL DSTERF( N, D, E, INFO )
350: GO TO 70
351: END IF
352: *
353: * If N is smaller than the minimum divide size (SMLSIZ+1), then
354: * solve the problem with another solver.
355: *
356: IF( N.LE.SMLSIZ ) THEN
357: *
358: CALL ZSTEQR( COMPZ, N, D, E, Z, LDZ, RWORK, INFO )
359: *
360: ELSE
361: *
362: * If COMPZ = 'I', we simply call DSTEDC instead.
363: *
364: IF( ICOMPZ.EQ.2 ) THEN
365: CALL DLASET( 'Full', N, N, ZERO, ONE, RWORK, N )
366: LL = N*N + 1
367: CALL DSTEDC( 'I', N, D, E, RWORK, N,
368: $ RWORK( LL ), LRWORK-LL+1, IWORK, LIWORK, INFO )
369: DO 20 J = 1, N
370: DO 10 I = 1, N
371: Z( I, J ) = RWORK( ( J-1 )*N+I )
372: 10 CONTINUE
373: 20 CONTINUE
374: GO TO 70
375: END IF
376: *
377: * From now on, only option left to be handled is COMPZ = 'V',
378: * i.e. ICOMPZ = 1.
379: *
380: * Scale.
381: *
382: ORGNRM = DLANST( 'M', N, D, E )
383: IF( ORGNRM.EQ.ZERO )
384: $ GO TO 70
385: *
386: EPS = DLAMCH( 'Epsilon' )
387: *
388: START = 1
389: *
390: * while ( START <= N )
391: *
392: 30 CONTINUE
393: IF( START.LE.N ) THEN
394: *
395: * Let FINISH be the position of the next subdiagonal entry
396: * such that E( FINISH ) <= TINY or FINISH = N if no such
397: * subdiagonal exists. The matrix identified by the elements
398: * between START and FINISH constitutes an independent
399: * sub-problem.
400: *
401: FINISH = START
402: 40 CONTINUE
403: IF( FINISH.LT.N ) THEN
404: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
405: $ SQRT( ABS( D( FINISH+1 ) ) )
406: IF( ABS( E( FINISH ) ).GT.TINY ) THEN
407: FINISH = FINISH + 1
408: GO TO 40
409: END IF
410: END IF
411: *
412: * (Sub) Problem determined. Compute its size and solve it.
413: *
414: M = FINISH - START + 1
415: IF( M.GT.SMLSIZ ) THEN
416: *
417: * Scale.
418: *
419: ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
420: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
421: $ INFO )
422: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
423: $ M-1, INFO )
424: *
425: CALL ZLAED0( N, M, D( START ), E( START ), Z( 1, START ),
426: $ LDZ, WORK, N, RWORK, IWORK, INFO )
427: IF( INFO.GT.0 ) THEN
428: INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
429: $ MOD( INFO, ( M+1 ) ) + START - 1
430: GO TO 70
431: END IF
432: *
433: * Scale back.
434: *
435: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
436: $ INFO )
437: *
438: ELSE
439: CALL DSTEQR( 'I', M, D( START ), E( START ), RWORK, M,
440: $ RWORK( M*M+1 ), INFO )
441: CALL ZLACRM( N, M, Z( 1, START ), LDZ, RWORK, M, WORK, N,
442: $ RWORK( M*M+1 ) )
443: CALL ZLACPY( 'A', N, M, WORK, N, Z( 1, START ), LDZ )
444: IF( INFO.GT.0 ) THEN
445: INFO = START*( N+1 ) + FINISH
446: GO TO 70
447: END IF
448: END IF
449: *
450: START = FINISH + 1
451: GO TO 30
452: END IF
453: *
454: * endwhile
455: *
456: *
457: * Use Selection Sort to minimize swaps of eigenvectors
458: *
459: DO 60 II = 2, N
460: I = II - 1
461: K = I
462: P = D( I )
463: DO 50 J = II, N
464: IF( D( J ).LT.P ) THEN
465: K = J
466: P = D( J )
467: END IF
468: 50 CONTINUE
469: IF( K.NE.I ) THEN
470: D( K ) = D( I )
471: D( I ) = P
472: CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
473: END IF
474: 60 CONTINUE
475: END IF
476: *
477: 70 CONTINUE
478: WORK( 1 ) = LWMIN
479: RWORK( 1 ) = LRWMIN
480: IWORK( 1 ) = LIWMIN
481: *
482: RETURN
483: *
484: * End of ZSTEDC
485: *
486: 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