Annotation of rpl/lapack/lapack/zstedc.f, revision 1.8
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: *
! 202: *> \date November 2011
! 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.8 ! bertrand 216: * -- LAPACK computational routine (version 3.4.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.8 ! bertrand 219: * November 2011
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: *
457: * If the problem split any number of times, then the eigenvalues
458: * will not be properly ordered. Here we permute the eigenvalues
459: * (and the associated eigenvectors) into ascending order.
460: *
461: IF( M.NE.N ) THEN
462: *
463: * Use Selection Sort to minimize swaps of eigenvectors
464: *
465: DO 60 II = 2, N
466: I = II - 1
467: K = I
468: P = D( I )
469: DO 50 J = II, N
470: IF( D( J ).LT.P ) THEN
471: K = J
472: P = D( J )
473: END IF
474: 50 CONTINUE
475: IF( K.NE.I ) THEN
476: D( K ) = D( I )
477: D( I ) = P
478: CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
479: END IF
480: 60 CONTINUE
481: END IF
482: END IF
483: *
484: 70 CONTINUE
485: WORK( 1 ) = LWMIN
486: RWORK( 1 ) = LRWMIN
487: IWORK( 1 ) = LIWMIN
488: *
489: RETURN
490: *
491: * End of ZSTEDC
492: *
493: 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