Annotation of rpl/lapack/lapack/zggesx.f, revision 1.8
1.8 ! bertrand 1: *> \brief <b> ZGGESX 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 ZGGESX + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggesx.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggesx.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggesx.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
! 22: * B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
! 23: * LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
! 24: * IWORK, LIWORK, BWORK, INFO )
! 25: *
! 26: * .. Scalar Arguments ..
! 27: * CHARACTER JOBVSL, JOBVSR, SENSE, SORT
! 28: * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
! 29: * $ SDIM
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * LOGICAL BWORK( * )
! 33: * INTEGER IWORK( * )
! 34: * DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
! 35: * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
! 36: * $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
! 37: * $ WORK( * )
! 38: * ..
! 39: * .. Function Arguments ..
! 40: * LOGICAL SELCTG
! 41: * EXTERNAL SELCTG
! 42: * ..
! 43: *
! 44: *
! 45: *> \par Purpose:
! 46: * =============
! 47: *>
! 48: *> \verbatim
! 49: *>
! 50: *> ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
! 51: *> (A,B), the generalized eigenvalues, the complex Schur form (S,T),
! 52: *> and, optionally, the left and/or right matrices of Schur vectors (VSL
! 53: *> and VSR). This gives the generalized Schur factorization
! 54: *>
! 55: *> (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
! 56: *>
! 57: *> where (VSR)**H is the conjugate-transpose of VSR.
! 58: *>
! 59: *> Optionally, it also orders the eigenvalues so that a selected cluster
! 60: *> of eigenvalues appears in the leading diagonal blocks of the upper
! 61: *> triangular matrix S and the upper triangular matrix T; computes
! 62: *> a reciprocal condition number for the average of the selected
! 63: *> eigenvalues (RCONDE); and computes a reciprocal condition number for
! 64: *> the right and left deflating subspaces corresponding to the selected
! 65: *> eigenvalues (RCONDV). The leading columns of VSL and VSR then form
! 66: *> an orthonormal basis for the corresponding left and right eigenspaces
! 67: *> (deflating subspaces).
! 68: *>
! 69: *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
! 70: *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
! 71: *> usually represented as the pair (alpha,beta), as there is a
! 72: *> reasonable interpretation for beta=0 or for both being zero.
! 73: *>
! 74: *> A pair of matrices (S,T) is in generalized complex Schur form if T is
! 75: *> upper triangular with non-negative diagonal and S is upper
! 76: *> triangular.
! 77: *> \endverbatim
! 78: *
! 79: * Arguments:
! 80: * ==========
! 81: *
! 82: *> \param[in] JOBVSL
! 83: *> \verbatim
! 84: *> JOBVSL is CHARACTER*1
! 85: *> = 'N': do not compute the left Schur vectors;
! 86: *> = 'V': compute the left Schur vectors.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] JOBVSR
! 90: *> \verbatim
! 91: *> JOBVSR is CHARACTER*1
! 92: *> = 'N': do not compute the right Schur vectors;
! 93: *> = 'V': compute the right Schur vectors.
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] SORT
! 97: *> \verbatim
! 98: *> SORT is CHARACTER*1
! 99: *> Specifies whether or not to order the eigenvalues on the
! 100: *> diagonal of the generalized Schur form.
! 101: *> = 'N': Eigenvalues are not ordered;
! 102: *> = 'S': Eigenvalues are ordered (see SELCTG).
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] SELCTG
! 106: *> \verbatim
! 107: *> SELCTG is procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
! 108: *> SELCTG must be declared EXTERNAL in the calling subroutine.
! 109: *> If SORT = 'N', SELCTG is not referenced.
! 110: *> If SORT = 'S', SELCTG is used to select eigenvalues to sort
! 111: *> to the top left of the Schur form.
! 112: *> Note that a selected complex eigenvalue may no longer satisfy
! 113: *> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
! 114: *> ordering may change the value of complex eigenvalues
! 115: *> (especially if the eigenvalue is ill-conditioned), in this
! 116: *> case INFO is set to N+3 see INFO below).
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] SENSE
! 120: *> \verbatim
! 121: *> SENSE is CHARACTER*1
! 122: *> Determines which reciprocal condition numbers are computed.
! 123: *> = 'N' : None are computed;
! 124: *> = 'E' : Computed for average of selected eigenvalues only;
! 125: *> = 'V' : Computed for selected deflating subspaces only;
! 126: *> = 'B' : Computed for both.
! 127: *> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[in] N
! 131: *> \verbatim
! 132: *> N is INTEGER
! 133: *> The order of the matrices A, B, VSL, and VSR. N >= 0.
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in,out] A
! 137: *> \verbatim
! 138: *> A is COMPLEX*16 array, dimension (LDA, N)
! 139: *> On entry, the first of the pair of matrices.
! 140: *> On exit, A has been overwritten by its generalized Schur
! 141: *> form S.
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[in] LDA
! 145: *> \verbatim
! 146: *> LDA is INTEGER
! 147: *> The leading dimension of A. LDA >= max(1,N).
! 148: *> \endverbatim
! 149: *>
! 150: *> \param[in,out] B
! 151: *> \verbatim
! 152: *> B is COMPLEX*16 array, dimension (LDB, N)
! 153: *> On entry, the second of the pair of matrices.
! 154: *> On exit, B has been overwritten by its generalized Schur
! 155: *> form T.
! 156: *> \endverbatim
! 157: *>
! 158: *> \param[in] LDB
! 159: *> \verbatim
! 160: *> LDB is INTEGER
! 161: *> The leading dimension of B. LDB >= max(1,N).
! 162: *> \endverbatim
! 163: *>
! 164: *> \param[out] SDIM
! 165: *> \verbatim
! 166: *> SDIM is INTEGER
! 167: *> If SORT = 'N', SDIM = 0.
! 168: *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
! 169: *> for which SELCTG is true.
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[out] ALPHA
! 173: *> \verbatim
! 174: *> ALPHA is COMPLEX*16 array, dimension (N)
! 175: *> \endverbatim
! 176: *>
! 177: *> \param[out] BETA
! 178: *> \verbatim
! 179: *> BETA is COMPLEX*16 array, dimension (N)
! 180: *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
! 181: *> generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
! 182: *> the diagonals of the complex Schur form (S,T). BETA(j) will
! 183: *> be non-negative real.
! 184: *>
! 185: *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
! 186: *> underflow, and BETA(j) may even be zero. Thus, the user
! 187: *> should avoid naively computing the ratio alpha/beta.
! 188: *> However, ALPHA will be always less than and usually
! 189: *> comparable with norm(A) in magnitude, and BETA always less
! 190: *> than and usually comparable with norm(B).
! 191: *> \endverbatim
! 192: *>
! 193: *> \param[out] VSL
! 194: *> \verbatim
! 195: *> VSL is COMPLEX*16 array, dimension (LDVSL,N)
! 196: *> If JOBVSL = 'V', VSL will contain the left Schur vectors.
! 197: *> Not referenced if JOBVSL = 'N'.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[in] LDVSL
! 201: *> \verbatim
! 202: *> LDVSL is INTEGER
! 203: *> The leading dimension of the matrix VSL. LDVSL >=1, and
! 204: *> if JOBVSL = 'V', LDVSL >= N.
! 205: *> \endverbatim
! 206: *>
! 207: *> \param[out] VSR
! 208: *> \verbatim
! 209: *> VSR is COMPLEX*16 array, dimension (LDVSR,N)
! 210: *> If JOBVSR = 'V', VSR will contain the right Schur vectors.
! 211: *> Not referenced if JOBVSR = 'N'.
! 212: *> \endverbatim
! 213: *>
! 214: *> \param[in] LDVSR
! 215: *> \verbatim
! 216: *> LDVSR is INTEGER
! 217: *> The leading dimension of the matrix VSR. LDVSR >= 1, and
! 218: *> if JOBVSR = 'V', LDVSR >= N.
! 219: *> \endverbatim
! 220: *>
! 221: *> \param[out] RCONDE
! 222: *> \verbatim
! 223: *> RCONDE is DOUBLE PRECISION array, dimension ( 2 )
! 224: *> If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
! 225: *> reciprocal condition numbers for the average of the selected
! 226: *> eigenvalues.
! 227: *> Not referenced if SENSE = 'N' or 'V'.
! 228: *> \endverbatim
! 229: *>
! 230: *> \param[out] RCONDV
! 231: *> \verbatim
! 232: *> RCONDV is DOUBLE PRECISION array, dimension ( 2 )
! 233: *> If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
! 234: *> reciprocal condition number for the selected deflating
! 235: *> subspaces.
! 236: *> Not referenced if SENSE = 'N' or 'E'.
! 237: *> \endverbatim
! 238: *>
! 239: *> \param[out] WORK
! 240: *> \verbatim
! 241: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 242: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 243: *> \endverbatim
! 244: *>
! 245: *> \param[in] LWORK
! 246: *> \verbatim
! 247: *> LWORK is INTEGER
! 248: *> The dimension of the array WORK.
! 249: *> If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
! 250: *> LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
! 251: *> LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2.
! 252: *> Note also that an error is only returned if
! 253: *> LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may
! 254: *> not be large enough.
! 255: *>
! 256: *> If LWORK = -1, then a workspace query is assumed; the routine
! 257: *> only calculates the bound on the optimal size of the WORK
! 258: *> array and the minimum size of the IWORK array, returns these
! 259: *> values as the first entries of the WORK and IWORK arrays, and
! 260: *> no error message related to LWORK or LIWORK is issued by
! 261: *> XERBLA.
! 262: *> \endverbatim
! 263: *>
! 264: *> \param[out] RWORK
! 265: *> \verbatim
! 266: *> RWORK is DOUBLE PRECISION array, dimension ( 8*N )
! 267: *> Real workspace.
! 268: *> \endverbatim
! 269: *>
! 270: *> \param[out] IWORK
! 271: *> \verbatim
! 272: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 273: *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
! 274: *> \endverbatim
! 275: *>
! 276: *> \param[in] LIWORK
! 277: *> \verbatim
! 278: *> LIWORK is INTEGER
! 279: *> The dimension of the array IWORK.
! 280: *> If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
! 281: *> LIWORK >= N+2.
! 282: *>
! 283: *> If LIWORK = -1, then a workspace query is assumed; the
! 284: *> routine only calculates the bound on the optimal size of the
! 285: *> WORK array and the minimum size of the IWORK array, returns
! 286: *> these values as the first entries of the WORK and IWORK
! 287: *> arrays, and no error message related to LWORK or LIWORK is
! 288: *> issued by XERBLA.
! 289: *> \endverbatim
! 290: *>
! 291: *> \param[out] BWORK
! 292: *> \verbatim
! 293: *> BWORK is LOGICAL array, dimension (N)
! 294: *> Not referenced if SORT = 'N'.
! 295: *> \endverbatim
! 296: *>
! 297: *> \param[out] INFO
! 298: *> \verbatim
! 299: *> INFO is INTEGER
! 300: *> = 0: successful exit
! 301: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 302: *> = 1,...,N:
! 303: *> The QZ iteration failed. (A,B) are not in Schur
! 304: *> form, but ALPHA(j) and BETA(j) should be correct for
! 305: *> j=INFO+1,...,N.
! 306: *> > N: =N+1: other than QZ iteration failed in ZHGEQZ
! 307: *> =N+2: after reordering, roundoff changed values of
! 308: *> some complex eigenvalues so that leading
! 309: *> eigenvalues in the Generalized Schur form no
! 310: *> longer satisfy SELCTG=.TRUE. This could also
! 311: *> be caused due to scaling.
! 312: *> =N+3: reordering failed in ZTGSEN.
! 313: *> \endverbatim
! 314: *
! 315: * Authors:
! 316: * ========
! 317: *
! 318: *> \author Univ. of Tennessee
! 319: *> \author Univ. of California Berkeley
! 320: *> \author Univ. of Colorado Denver
! 321: *> \author NAG Ltd.
! 322: *
! 323: *> \date November 2011
! 324: *
! 325: *> \ingroup complex16GEeigen
! 326: *
! 327: * =====================================================================
1.1 bertrand 328: SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
329: $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
330: $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
331: $ IWORK, LIWORK, BWORK, INFO )
332: *
1.8 ! bertrand 333: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 334: * -- LAPACK is a software package provided by Univ. of Tennessee, --
335: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 336: * November 2011
1.1 bertrand 337: *
338: * .. Scalar Arguments ..
339: CHARACTER JOBVSL, JOBVSR, SENSE, SORT
340: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
341: $ SDIM
342: * ..
343: * .. Array Arguments ..
344: LOGICAL BWORK( * )
345: INTEGER IWORK( * )
346: DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
347: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
348: $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
349: $ WORK( * )
350: * ..
351: * .. Function Arguments ..
352: LOGICAL SELCTG
353: EXTERNAL SELCTG
354: * ..
355: *
356: * =====================================================================
357: *
358: * .. Parameters ..
359: DOUBLE PRECISION ZERO, ONE
360: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
361: COMPLEX*16 CZERO, CONE
362: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
363: $ CONE = ( 1.0D+0, 0.0D+0 ) )
364: * ..
365: * .. Local Scalars ..
366: LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
367: $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV
368: INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
369: $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK,
370: $ LIWMIN, LWRK, MAXWRK, MINWRK
371: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
372: $ PR, SMLNUM
373: * ..
374: * .. Local Arrays ..
375: DOUBLE PRECISION DIF( 2 )
376: * ..
377: * .. External Subroutines ..
378: EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
379: $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
380: $ ZUNMQR
381: * ..
382: * .. External Functions ..
383: LOGICAL LSAME
384: INTEGER ILAENV
385: DOUBLE PRECISION DLAMCH, ZLANGE
386: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
387: * ..
388: * .. Intrinsic Functions ..
389: INTRINSIC MAX, SQRT
390: * ..
391: * .. Executable Statements ..
392: *
393: * Decode the input arguments
394: *
395: IF( LSAME( JOBVSL, 'N' ) ) THEN
396: IJOBVL = 1
397: ILVSL = .FALSE.
398: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
399: IJOBVL = 2
400: ILVSL = .TRUE.
401: ELSE
402: IJOBVL = -1
403: ILVSL = .FALSE.
404: END IF
405: *
406: IF( LSAME( JOBVSR, 'N' ) ) THEN
407: IJOBVR = 1
408: ILVSR = .FALSE.
409: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
410: IJOBVR = 2
411: ILVSR = .TRUE.
412: ELSE
413: IJOBVR = -1
414: ILVSR = .FALSE.
415: END IF
416: *
417: WANTST = LSAME( SORT, 'S' )
418: WANTSN = LSAME( SENSE, 'N' )
419: WANTSE = LSAME( SENSE, 'E' )
420: WANTSV = LSAME( SENSE, 'V' )
421: WANTSB = LSAME( SENSE, 'B' )
422: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
423: IF( WANTSN ) THEN
424: IJOB = 0
425: ELSE IF( WANTSE ) THEN
426: IJOB = 1
427: ELSE IF( WANTSV ) THEN
428: IJOB = 2
429: ELSE IF( WANTSB ) THEN
430: IJOB = 4
431: END IF
432: *
433: * Test the input arguments
434: *
435: INFO = 0
436: IF( IJOBVL.LE.0 ) THEN
437: INFO = -1
438: ELSE IF( IJOBVR.LE.0 ) THEN
439: INFO = -2
440: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
441: INFO = -3
442: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
443: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
444: INFO = -5
445: ELSE IF( N.LT.0 ) THEN
446: INFO = -6
447: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
448: INFO = -8
449: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
450: INFO = -10
451: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
452: INFO = -15
453: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
454: INFO = -17
455: END IF
456: *
457: * Compute workspace
458: * (Note: Comments in the code beginning "Workspace:" describe the
459: * minimal amount of workspace needed at that point in the code,
460: * as well as the preferred amount for good performance.
461: * NB refers to the optimal block size for the immediately
462: * following subroutine, as returned by ILAENV.)
463: *
464: IF( INFO.EQ.0 ) THEN
465: IF( N.GT.0) THEN
466: MINWRK = 2*N
467: MAXWRK = N*(1 + ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
468: MAXWRK = MAX( MAXWRK, N*( 1 +
469: $ ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) ) )
470: IF( ILVSL ) THEN
471: MAXWRK = MAX( MAXWRK, N*( 1 +
472: $ ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) ) )
473: END IF
474: LWRK = MAXWRK
475: IF( IJOB.GE.1 )
476: $ LWRK = MAX( LWRK, N*N/2 )
477: ELSE
478: MINWRK = 1
479: MAXWRK = 1
480: LWRK = 1
481: END IF
482: WORK( 1 ) = LWRK
483: IF( WANTSN .OR. N.EQ.0 ) THEN
484: LIWMIN = 1
485: ELSE
486: LIWMIN = N + 2
487: END IF
488: IWORK( 1 ) = LIWMIN
489: *
490: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
491: INFO = -21
492: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN
493: INFO = -24
494: END IF
495: END IF
496: *
497: IF( INFO.NE.0 ) THEN
498: CALL XERBLA( 'ZGGESX', -INFO )
499: RETURN
500: ELSE IF (LQUERY) THEN
501: RETURN
502: END IF
503: *
504: * Quick return if possible
505: *
506: IF( N.EQ.0 ) THEN
507: SDIM = 0
508: RETURN
509: END IF
510: *
511: * Get machine constants
512: *
513: EPS = DLAMCH( 'P' )
514: SMLNUM = DLAMCH( 'S' )
515: BIGNUM = ONE / SMLNUM
516: CALL DLABAD( SMLNUM, BIGNUM )
517: SMLNUM = SQRT( SMLNUM ) / EPS
518: BIGNUM = ONE / SMLNUM
519: *
520: * Scale A if max element outside range [SMLNUM,BIGNUM]
521: *
522: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
523: ILASCL = .FALSE.
524: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
525: ANRMTO = SMLNUM
526: ILASCL = .TRUE.
527: ELSE IF( ANRM.GT.BIGNUM ) THEN
528: ANRMTO = BIGNUM
529: ILASCL = .TRUE.
530: END IF
531: IF( ILASCL )
532: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
533: *
534: * Scale B if max element outside range [SMLNUM,BIGNUM]
535: *
536: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
537: ILBSCL = .FALSE.
538: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
539: BNRMTO = SMLNUM
540: ILBSCL = .TRUE.
541: ELSE IF( BNRM.GT.BIGNUM ) THEN
542: BNRMTO = BIGNUM
543: ILBSCL = .TRUE.
544: END IF
545: IF( ILBSCL )
546: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
547: *
548: * Permute the matrix to make it more nearly triangular
549: * (Real Workspace: need 6*N)
550: *
551: ILEFT = 1
552: IRIGHT = N + 1
553: IRWRK = IRIGHT + N
554: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
555: $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
556: *
557: * Reduce B to triangular form (QR decomposition of B)
558: * (Complex Workspace: need N, prefer N*NB)
559: *
560: IROWS = IHI + 1 - ILO
561: ICOLS = N + 1 - ILO
562: ITAU = 1
563: IWRK = ITAU + IROWS
564: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
565: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
566: *
567: * Apply the unitary transformation to matrix A
568: * (Complex Workspace: need N, prefer N*NB)
569: *
570: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
571: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
572: $ LWORK+1-IWRK, IERR )
573: *
574: * Initialize VSL
575: * (Complex Workspace: need N, prefer N*NB)
576: *
577: IF( ILVSL ) THEN
578: CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
579: IF( IROWS.GT.1 ) THEN
580: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
581: $ VSL( ILO+1, ILO ), LDVSL )
582: END IF
583: CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
584: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
585: END IF
586: *
587: * Initialize VSR
588: *
589: IF( ILVSR )
590: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
591: *
592: * Reduce to generalized Hessenberg form
593: * (Workspace: none needed)
594: *
595: CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
596: $ LDVSL, VSR, LDVSR, IERR )
597: *
598: SDIM = 0
599: *
600: * Perform QZ algorithm, computing Schur vectors if desired
601: * (Complex Workspace: need N)
602: * (Real Workspace: need N)
603: *
604: IWRK = ITAU
605: CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
606: $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
607: $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
608: IF( IERR.NE.0 ) THEN
609: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
610: INFO = IERR
611: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
612: INFO = IERR - N
613: ELSE
614: INFO = N + 1
615: END IF
616: GO TO 40
617: END IF
618: *
619: * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
620: * condition number(s)
621: *
622: IF( WANTST ) THEN
623: *
624: * Undo scaling on eigenvalues before SELCTGing
625: *
626: IF( ILASCL )
627: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
628: IF( ILBSCL )
629: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
630: *
631: * Select eigenvalues
632: *
633: DO 10 I = 1, N
634: BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
635: 10 CONTINUE
636: *
637: * Reorder eigenvalues, transform Generalized Schur vectors, and
638: * compute reciprocal condition numbers
639: * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
640: * otherwise, need 1 )
641: *
642: CALL ZTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
643: $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR,
644: $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK,
645: $ IERR )
646: *
647: IF( IJOB.GE.1 )
648: $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
649: IF( IERR.EQ.-21 ) THEN
650: *
651: * not enough complex workspace
652: *
653: INFO = -21
654: ELSE
655: IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
656: RCONDE( 1 ) = PL
657: RCONDE( 2 ) = PR
658: END IF
659: IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
660: RCONDV( 1 ) = DIF( 1 )
661: RCONDV( 2 ) = DIF( 2 )
662: END IF
663: IF( IERR.EQ.1 )
664: $ INFO = N + 3
665: END IF
666: *
667: END IF
668: *
669: * Apply permutation to VSL and VSR
670: * (Workspace: none needed)
671: *
672: IF( ILVSL )
673: $ CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
674: $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
675: *
676: IF( ILVSR )
677: $ CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
678: $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
679: *
680: * Undo scaling
681: *
682: IF( ILASCL ) THEN
683: CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
684: CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
685: END IF
686: *
687: IF( ILBSCL ) THEN
688: CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
689: CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
690: END IF
691: *
692: IF( WANTST ) THEN
693: *
694: * Check if reordering is correct
695: *
696: LASTSL = .TRUE.
697: SDIM = 0
698: DO 30 I = 1, N
699: CURSL = SELCTG( ALPHA( I ), BETA( I ) )
700: IF( CURSL )
701: $ SDIM = SDIM + 1
702: IF( CURSL .AND. .NOT.LASTSL )
703: $ INFO = N + 2
704: LASTSL = CURSL
705: 30 CONTINUE
706: *
707: END IF
708: *
709: 40 CONTINUE
710: *
711: WORK( 1 ) = MAXWRK
712: IWORK( 1 ) = LIWMIN
713: *
714: RETURN
715: *
716: * End of ZGGESX
717: *
718: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack