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