Annotation of rpl/lapack/lapack/zgeesx.f, revision 1.9
1.9 ! bertrand 1: *> \brief <b> ZGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGEESX + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeesx.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeesx.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeesx.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
! 22: * VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
! 23: * BWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER JOBVS, SENSE, SORT
! 27: * INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
! 28: * DOUBLE PRECISION RCONDE, RCONDV
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * LOGICAL BWORK( * )
! 32: * DOUBLE PRECISION RWORK( * )
! 33: * COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
! 34: * ..
! 35: * .. Function Arguments ..
! 36: * LOGICAL SELECT
! 37: * EXTERNAL SELECT
! 38: * ..
! 39: *
! 40: *
! 41: *> \par Purpose:
! 42: * =============
! 43: *>
! 44: *> \verbatim
! 45: *>
! 46: *> ZGEESX computes for an N-by-N complex nonsymmetric matrix A, the
! 47: *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
! 48: *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
! 49: *>
! 50: *> Optionally, it also orders the eigenvalues on the diagonal of the
! 51: *> Schur form so that selected eigenvalues are at the top left;
! 52: *> computes a reciprocal condition number for the average of the
! 53: *> selected eigenvalues (RCONDE); and computes a reciprocal condition
! 54: *> number for the right invariant subspace corresponding to the
! 55: *> selected eigenvalues (RCONDV). The leading columns of Z form an
! 56: *> orthonormal basis for this invariant subspace.
! 57: *>
! 58: *> For further explanation of the reciprocal condition numbers RCONDE
! 59: *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
! 60: *> these quantities are called s and sep respectively).
! 61: *>
! 62: *> A complex matrix is in Schur form if it is upper triangular.
! 63: *> \endverbatim
! 64: *
! 65: * Arguments:
! 66: * ==========
! 67: *
! 68: *> \param[in] JOBVS
! 69: *> \verbatim
! 70: *> JOBVS is CHARACTER*1
! 71: *> = 'N': Schur vectors are not computed;
! 72: *> = 'V': Schur vectors are computed.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] SORT
! 76: *> \verbatim
! 77: *> SORT is CHARACTER*1
! 78: *> Specifies whether or not to order the eigenvalues on the
! 79: *> diagonal of the Schur form.
! 80: *> = 'N': Eigenvalues are not ordered;
! 81: *> = 'S': Eigenvalues are ordered (see SELECT).
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] SELECT
! 85: *> \verbatim
! 86: *> SELECT is procedure) LOGICAL FUNCTION of one COMPLEX*16 argument
! 87: *> SELECT must be declared EXTERNAL in the calling subroutine.
! 88: *> If SORT = 'S', SELECT is used to select eigenvalues to order
! 89: *> to the top left of the Schur form.
! 90: *> If SORT = 'N', SELECT is not referenced.
! 91: *> An eigenvalue W(j) is selected if SELECT(W(j)) is true.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] SENSE
! 95: *> \verbatim
! 96: *> SENSE is CHARACTER*1
! 97: *> Determines which reciprocal condition numbers are computed.
! 98: *> = 'N': None are computed;
! 99: *> = 'E': Computed for average of selected eigenvalues only;
! 100: *> = 'V': Computed for selected right invariant subspace only;
! 101: *> = 'B': Computed for both.
! 102: *> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] N
! 106: *> \verbatim
! 107: *> N is INTEGER
! 108: *> The order of the matrix A. N >= 0.
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[in,out] A
! 112: *> \verbatim
! 113: *> A is COMPLEX*16 array, dimension (LDA, N)
! 114: *> On entry, the N-by-N matrix A.
! 115: *> On exit, A is overwritten by its Schur form T.
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[in] LDA
! 119: *> \verbatim
! 120: *> LDA is INTEGER
! 121: *> The leading dimension of the array A. LDA >= max(1,N).
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[out] SDIM
! 125: *> \verbatim
! 126: *> SDIM is INTEGER
! 127: *> If SORT = 'N', SDIM = 0.
! 128: *> If SORT = 'S', SDIM = number of eigenvalues for which
! 129: *> SELECT is true.
! 130: *> \endverbatim
! 131: *>
! 132: *> \param[out] W
! 133: *> \verbatim
! 134: *> W is COMPLEX*16 array, dimension (N)
! 135: *> W contains the computed eigenvalues, in the same order
! 136: *> that they appear on the diagonal of the output Schur form T.
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[out] VS
! 140: *> \verbatim
! 141: *> VS is COMPLEX*16 array, dimension (LDVS,N)
! 142: *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
! 143: *> vectors.
! 144: *> If JOBVS = 'N', VS is not referenced.
! 145: *> \endverbatim
! 146: *>
! 147: *> \param[in] LDVS
! 148: *> \verbatim
! 149: *> LDVS is INTEGER
! 150: *> The leading dimension of the array VS. LDVS >= 1, and if
! 151: *> JOBVS = 'V', LDVS >= N.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[out] RCONDE
! 155: *> \verbatim
! 156: *> RCONDE is DOUBLE PRECISION
! 157: *> If SENSE = 'E' or 'B', RCONDE contains the reciprocal
! 158: *> condition number for the average of the selected eigenvalues.
! 159: *> Not referenced if SENSE = 'N' or 'V'.
! 160: *> \endverbatim
! 161: *>
! 162: *> \param[out] RCONDV
! 163: *> \verbatim
! 164: *> RCONDV is DOUBLE PRECISION
! 165: *> If SENSE = 'V' or 'B', RCONDV contains the reciprocal
! 166: *> condition number for the selected right invariant subspace.
! 167: *> Not referenced if SENSE = 'N' or 'E'.
! 168: *> \endverbatim
! 169: *>
! 170: *> \param[out] WORK
! 171: *> \verbatim
! 172: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 173: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 174: *> \endverbatim
! 175: *>
! 176: *> \param[in] LWORK
! 177: *> \verbatim
! 178: *> LWORK is INTEGER
! 179: *> The dimension of the array WORK. LWORK >= max(1,2*N).
! 180: *> Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
! 181: *> where SDIM is the number of selected eigenvalues computed by
! 182: *> this routine. Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
! 183: *> that an error is only returned if LWORK < max(1,2*N), but if
! 184: *> SENSE = 'E' or 'V' or 'B' this may not be large enough.
! 185: *> For good performance, LWORK must generally be larger.
! 186: *>
! 187: *> If LWORK = -1, then a workspace query is assumed; the routine
! 188: *> only calculates upper bound on the optimal size of the
! 189: *> array WORK, returns this value as the first entry of the WORK
! 190: *> array, and no error message related to LWORK is issued by
! 191: *> XERBLA.
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[out] RWORK
! 195: *> \verbatim
! 196: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 197: *> \endverbatim
! 198: *>
! 199: *> \param[out] BWORK
! 200: *> \verbatim
! 201: *> BWORK is LOGICAL array, dimension (N)
! 202: *> Not referenced if SORT = 'N'.
! 203: *> \endverbatim
! 204: *>
! 205: *> \param[out] INFO
! 206: *> \verbatim
! 207: *> INFO is INTEGER
! 208: *> = 0: successful exit
! 209: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 210: *> > 0: if INFO = i, and i is
! 211: *> <= N: the QR algorithm failed to compute all the
! 212: *> eigenvalues; elements 1:ILO-1 and i+1:N of W
! 213: *> contain those eigenvalues which have converged; if
! 214: *> JOBVS = 'V', VS contains the transformation which
! 215: *> reduces A to its partially converged Schur form.
! 216: *> = N+1: the eigenvalues could not be reordered because some
! 217: *> eigenvalues were too close to separate (the problem
! 218: *> is very ill-conditioned);
! 219: *> = N+2: after reordering, roundoff changed values of some
! 220: *> complex eigenvalues so that leading eigenvalues in
! 221: *> the Schur form no longer satisfy SELECT=.TRUE. This
! 222: *> could also be caused by underflow due to scaling.
! 223: *> \endverbatim
! 224: *
! 225: * Authors:
! 226: * ========
! 227: *
! 228: *> \author Univ. of Tennessee
! 229: *> \author Univ. of California Berkeley
! 230: *> \author Univ. of Colorado Denver
! 231: *> \author NAG Ltd.
! 232: *
! 233: *> \date November 2011
! 234: *
! 235: *> \ingroup complex16GEeigen
! 236: *
! 237: * =====================================================================
1.1 bertrand 238: SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
239: $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
240: $ BWORK, INFO )
241: *
1.9 ! bertrand 242: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 243: * -- LAPACK is a software package provided by Univ. of Tennessee, --
244: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 245: * November 2011
1.1 bertrand 246: *
247: * .. Scalar Arguments ..
248: CHARACTER JOBVS, SENSE, SORT
249: INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
250: DOUBLE PRECISION RCONDE, RCONDV
251: * ..
252: * .. Array Arguments ..
253: LOGICAL BWORK( * )
254: DOUBLE PRECISION RWORK( * )
255: COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
256: * ..
257: * .. Function Arguments ..
258: LOGICAL SELECT
259: EXTERNAL SELECT
260: * ..
261: *
262: * =====================================================================
263: *
264: * .. Parameters ..
265: DOUBLE PRECISION ZERO, ONE
266: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
267: * ..
268: * .. Local Scalars ..
1.5 bertrand 269: LOGICAL LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
270: $ WANTSV, WANTVS
1.1 bertrand 271: INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
272: $ ITAU, IWRK, LWRK, MAXWRK, MINWRK
273: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
274: * ..
275: * .. Local Arrays ..
276: DOUBLE PRECISION DUM( 1 )
277: * ..
278: * .. External Subroutines ..
279: EXTERNAL DLABAD, DLASCL, XERBLA, ZCOPY, ZGEBAK, ZGEBAL,
280: $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
281: * ..
282: * .. External Functions ..
283: LOGICAL LSAME
284: INTEGER ILAENV
285: DOUBLE PRECISION DLAMCH, ZLANGE
286: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
287: * ..
288: * .. Intrinsic Functions ..
289: INTRINSIC MAX, SQRT
290: * ..
291: * .. Executable Statements ..
292: *
293: * Test the input arguments
294: *
295: INFO = 0
296: WANTVS = LSAME( JOBVS, 'V' )
297: WANTST = LSAME( SORT, 'S' )
298: WANTSN = LSAME( SENSE, 'N' )
299: WANTSE = LSAME( SENSE, 'E' )
300: WANTSV = LSAME( SENSE, 'V' )
301: WANTSB = LSAME( SENSE, 'B' )
1.5 bertrand 302: LQUERY = ( LWORK.EQ.-1 )
303: *
1.1 bertrand 304: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
305: INFO = -1
306: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
307: INFO = -2
308: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
309: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
310: INFO = -4
311: ELSE IF( N.LT.0 ) THEN
312: INFO = -5
313: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
314: INFO = -7
315: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
316: INFO = -11
317: END IF
318: *
319: * Compute workspace
320: * (Note: Comments in the code beginning "Workspace:" describe the
321: * minimal amount of real workspace needed at that point in the
322: * code, as well as the preferred amount for good performance.
323: * CWorkspace refers to complex workspace, and RWorkspace to real
324: * workspace. NB refers to the optimal block size for the
325: * immediately following subroutine, as returned by ILAENV.
326: * HSWORK refers to the workspace preferred by ZHSEQR, as
327: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
328: * the worst case.
329: * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
330: * depends on SDIM, which is computed by the routine ZTRSEN later
331: * in the code.)
332: *
333: IF( INFO.EQ.0 ) THEN
334: IF( N.EQ.0 ) THEN
335: MINWRK = 1
336: LWRK = 1
337: ELSE
338: MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
339: MINWRK = 2*N
340: *
341: CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
342: $ WORK, -1, IEVAL )
343: HSWORK = WORK( 1 )
344: *
345: IF( .NOT.WANTVS ) THEN
346: MAXWRK = MAX( MAXWRK, HSWORK )
347: ELSE
348: MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
349: $ ' ', N, 1, N, -1 ) )
350: MAXWRK = MAX( MAXWRK, HSWORK )
351: END IF
352: LWRK = MAXWRK
353: IF( .NOT.WANTSN )
354: $ LWRK = MAX( LWRK, ( N*N )/2 )
355: END IF
356: WORK( 1 ) = LWRK
357: *
1.5 bertrand 358: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
1.1 bertrand 359: INFO = -15
360: END IF
361: END IF
362: *
363: IF( INFO.NE.0 ) THEN
364: CALL XERBLA( 'ZGEESX', -INFO )
365: RETURN
1.5 bertrand 366: ELSE IF( LQUERY ) THEN
367: RETURN
1.1 bertrand 368: END IF
369: *
370: * Quick return if possible
371: *
372: IF( N.EQ.0 ) THEN
373: SDIM = 0
374: RETURN
375: END IF
376: *
377: * Get machine constants
378: *
379: EPS = DLAMCH( 'P' )
380: SMLNUM = DLAMCH( 'S' )
381: BIGNUM = ONE / SMLNUM
382: CALL DLABAD( SMLNUM, BIGNUM )
383: SMLNUM = SQRT( SMLNUM ) / EPS
384: BIGNUM = ONE / SMLNUM
385: *
386: * Scale A if max element outside range [SMLNUM,BIGNUM]
387: *
388: ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
389: SCALEA = .FALSE.
390: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
391: SCALEA = .TRUE.
392: CSCALE = SMLNUM
393: ELSE IF( ANRM.GT.BIGNUM ) THEN
394: SCALEA = .TRUE.
395: CSCALE = BIGNUM
396: END IF
397: IF( SCALEA )
398: $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
399: *
400: *
401: * Permute the matrix to make it more nearly triangular
402: * (CWorkspace: none)
403: * (RWorkspace: need N)
404: *
405: IBAL = 1
406: CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
407: *
408: * Reduce to upper Hessenberg form
409: * (CWorkspace: need 2*N, prefer N+N*NB)
410: * (RWorkspace: none)
411: *
412: ITAU = 1
413: IWRK = N + ITAU
414: CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
415: $ LWORK-IWRK+1, IERR )
416: *
417: IF( WANTVS ) THEN
418: *
419: * Copy Householder vectors to VS
420: *
421: CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS )
422: *
423: * Generate unitary matrix in VS
424: * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
425: * (RWorkspace: none)
426: *
427: CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
428: $ LWORK-IWRK+1, IERR )
429: END IF
430: *
431: SDIM = 0
432: *
433: * Perform QR iteration, accumulating Schur vectors in VS if desired
434: * (CWorkspace: need 1, prefer HSWORK (see comments) )
435: * (RWorkspace: none)
436: *
437: IWRK = ITAU
438: CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
439: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
440: IF( IEVAL.GT.0 )
441: $ INFO = IEVAL
442: *
443: * Sort eigenvalues if desired
444: *
445: IF( WANTST .AND. INFO.EQ.0 ) THEN
446: IF( SCALEA )
447: $ CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
448: DO 10 I = 1, N
449: BWORK( I ) = SELECT( W( I ) )
450: 10 CONTINUE
451: *
452: * Reorder eigenvalues, transform Schur vectors, and compute
453: * reciprocal condition numbers
454: * (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
455: * otherwise, need none )
456: * (RWorkspace: none)
457: *
458: CALL ZTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
459: $ RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
460: $ ICOND )
461: IF( .NOT.WANTSN )
462: $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
463: IF( ICOND.EQ.-14 ) THEN
464: *
465: * Not enough complex workspace
466: *
467: INFO = -15
468: END IF
469: END IF
470: *
471: IF( WANTVS ) THEN
472: *
473: * Undo balancing
474: * (CWorkspace: none)
475: * (RWorkspace: need N)
476: *
477: CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
478: $ IERR )
479: END IF
480: *
481: IF( SCALEA ) THEN
482: *
483: * Undo scaling for the Schur form of A
484: *
485: CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
486: CALL ZCOPY( N, A, LDA+1, W, 1 )
487: IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
488: DUM( 1 ) = RCONDV
489: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
490: RCONDV = DUM( 1 )
491: END IF
492: END IF
493: *
494: WORK( 1 ) = MAXWRK
495: RETURN
496: *
497: * End of ZGEESX
498: *
499: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgeesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack