/src/libsesctherm/CLAPACK/SRC/dggevx.c
C | 843 lines | 462 code | 67 blank | 314 comment | 139 complexity | 022a21feccb9e65cff46d1967accdcc1 MD5 | raw file
Possible License(s): GPL-2.0
- #include "blaswrap.h"
- #include "f2c.h"
- /* Subroutine */ int dggevx_(char *balanc, char *jobvl, char *jobvr, char *
- sense, integer *n, doublereal *a, integer *lda, doublereal *b,
- integer *ldb, doublereal *alphar, doublereal *alphai, doublereal *
- beta, doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,
- integer *ilo, integer *ihi, doublereal *lscale, doublereal *rscale,
- doublereal *abnrm, doublereal *bbnrm, doublereal *rconde, doublereal *
- rcondv, doublereal *work, integer *lwork, integer *iwork, logical *
- bwork, integer *info)
- {
- /* -- LAPACK driver routine (version 3.0) --
- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
- Courant Institute, Argonne National Lab, and Rice University
- June 30, 1999
- Purpose
- =======
- DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
- the generalized eigenvalues, and optionally, the left and/or right
- generalized eigenvectors.
- Optionally also, it computes a balancing transformation to improve
- the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
- LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
- the eigenvalues (RCONDE), and reciprocal condition numbers for the
- right eigenvectors (RCONDV).
- A generalized eigenvalue for a pair of matrices (A,B) is a scalar
- lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
- singular. It is usually represented as the pair (alpha,beta), as
- there is a reasonable interpretation for beta=0, and even for both
- being zero.
- The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
- of (A,B) satisfies
- A * v(j) = lambda(j) * B * v(j) .
- The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
- of (A,B) satisfies
- u(j)**H * A = lambda(j) * u(j)**H * B.
- where u(j)**H is the conjugate-transpose of u(j).
- Arguments
- =========
- BALANC (input) CHARACTER*1
- Specifies the balance option to be performed.
- = 'N': do not diagonally scale or permute;
- = 'P': permute only;
- = 'S': scale only;
- = 'B': both permute and scale.
- Computed reciprocal condition numbers will be for the
- matrices after permuting and/or balancing. Permuting does
- not change condition numbers (in exact arithmetic), but
- balancing does.
- JOBVL (input) CHARACTER*1
- = 'N': do not compute the left generalized eigenvectors;
- = 'V': compute the left generalized eigenvectors.
- JOBVR (input) CHARACTER*1
- = 'N': do not compute the right generalized eigenvectors;
- = 'V': compute the right generalized eigenvectors.
- SENSE (input) CHARACTER*1
- Determines which reciprocal condition numbers are computed.
- = 'N': none are computed;
- = 'E': computed for eigenvalues only;
- = 'V': computed for eigenvectors only;
- = 'B': computed for eigenvalues and eigenvectors.
- N (input) INTEGER
- The order of the matrices A, B, VL, and VR. N >= 0.
- A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
- On entry, the matrix A in the pair (A,B).
- On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
- or both, then A contains the first part of the real Schur
- form of the "balanced" versions of the input A and B.
- LDA (input) INTEGER
- The leading dimension of A. LDA >= max(1,N).
- B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
- On entry, the matrix B in the pair (A,B).
- On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
- or both, then B contains the second part of the real Schur
- form of the "balanced" versions of the input A and B.
- LDB (input) INTEGER
- The leading dimension of B. LDB >= max(1,N).
- ALPHAR (output) DOUBLE PRECISION array, dimension (N)
- ALPHAI (output) DOUBLE PRECISION array, dimension (N)
- BETA (output) DOUBLE PRECISION array, dimension (N)
- On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
- be the generalized eigenvalues. If ALPHAI(j) is zero, then
- the j-th eigenvalue is real; if positive, then the j-th and
- (j+1)-st eigenvalues are a complex conjugate pair, with
- ALPHAI(j+1) negative.
- Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
- may easily over- or underflow, and BETA(j) may even be zero.
- Thus, the user should avoid naively computing the ratio
- ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
- than and usually comparable with norm(A) in magnitude, and
- BETA always less than and usually comparable with norm(B).
- VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
- If JOBVL = 'V', the left eigenvectors u(j) are stored one
- after another in the columns of VL, in the same order as
- their eigenvalues. If the j-th eigenvalue is real, then
- u(j) = VL(:,j), the j-th column of VL. If the j-th and
- (j+1)-th eigenvalues form a complex conjugate pair, then
- u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
- Each eigenvector will be scaled so the largest component have
- abs(real part) + abs(imag. part) = 1.
- Not referenced if JOBVL = 'N'.
- LDVL (input) INTEGER
- The leading dimension of the matrix VL. LDVL >= 1, and
- if JOBVL = 'V', LDVL >= N.
- VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
- If JOBVR = 'V', the right eigenvectors v(j) are stored one
- after another in the columns of VR, in the same order as
- their eigenvalues. If the j-th eigenvalue is real, then
- v(j) = VR(:,j), the j-th column of VR. If the j-th and
- (j+1)-th eigenvalues form a complex conjugate pair, then
- v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
- Each eigenvector will be scaled so the largest component have
- abs(real part) + abs(imag. part) = 1.
- Not referenced if JOBVR = 'N'.
- LDVR (input) INTEGER
- The leading dimension of the matrix VR. LDVR >= 1, and
- if JOBVR = 'V', LDVR >= N.
- ILO,IHI (output) INTEGER
- ILO and IHI are integer values such that on exit
- A(i,j) = 0 and B(i,j) = 0 if i > j and
- j = 1,...,ILO-1 or i = IHI+1,...,N.
- If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
- LSCALE (output) DOUBLE PRECISION array, dimension (N)
- Details of the permutations and scaling factors applied
- to the left side of A and B. If PL(j) is the index of the
- row interchanged with row j, and DL(j) is the scaling
- factor applied to row j, then
- LSCALE(j) = PL(j) for j = 1,...,ILO-1
- = DL(j) for j = ILO,...,IHI
- = PL(j) for j = IHI+1,...,N.
- The order in which the interchanges are made is N to IHI+1,
- then 1 to ILO-1.
- RSCALE (output) DOUBLE PRECISION array, dimension (N)
- Details of the permutations and scaling factors applied
- to the right side of A and B. If PR(j) is the index of the
- column interchanged with column j, and DR(j) is the scaling
- factor applied to column j, then
- RSCALE(j) = PR(j) for j = 1,...,ILO-1
- = DR(j) for j = ILO,...,IHI
- = PR(j) for j = IHI+1,...,N
- The order in which the interchanges are made is N to IHI+1,
- then 1 to ILO-1.
- ABNRM (output) DOUBLE PRECISION
- The one-norm of the balanced matrix A.
- BBNRM (output) DOUBLE PRECISION
- The one-norm of the balanced matrix B.
- RCONDE (output) DOUBLE PRECISION array, dimension (N)
- If SENSE = 'E' or 'B', the reciprocal condition numbers of
- the selected eigenvalues, stored in consecutive elements of
- the array. For a complex conjugate pair of eigenvalues two
- consecutive elements of RCONDE are set to the same value.
- Thus RCONDE(j), RCONDV(j), and the j-th columns of VL and VR
- all correspond to the same eigenpair (but not in general the
- j-th eigenpair, unless all eigenpairs are selected).
- If SENSE = 'V', RCONDE is not referenced.
- RCONDV (output) DOUBLE PRECISION array, dimension (N)
- If SENSE = 'V' or 'B', the estimated reciprocal condition
- numbers of the selected eigenvectors, stored in consecutive
- elements of the array. For a complex eigenvector two
- consecutive elements of RCONDV are set to the same value. If
- the eigenvalues cannot be reordered to compute RCONDV(j),
- RCONDV(j) is set to 0; this can only occur when the true
- value would be very small anyway.
- If SENSE = 'E', RCONDV is not referenced.
- WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
- On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
- LWORK (input) INTEGER
- The dimension of the array WORK. LWORK >= max(1,6*N).
- If SENSE = 'E', LWORK >= 12*N.
- If SENSE = 'V' or 'B', LWORK >= 2*N*N+12*N+16.
- If LWORK = -1, then a workspace query is assumed; the routine
- only calculates the optimal size of the WORK array, returns
- this value as the first entry of the WORK array, and no error
- message related to LWORK is issued by XERBLA.
- IWORK (workspace) INTEGER array, dimension (N+6)
- If SENSE = 'E', IWORK is not referenced.
- BWORK (workspace) LOGICAL array, dimension (N)
- If SENSE = 'N', BWORK is not referenced.
- INFO (output) INTEGER
- = 0: successful exit
- < 0: if INFO = -i, the i-th argument had an illegal value.
- = 1,...,N:
- The QZ iteration failed. No eigenvectors have been
- calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
- should be correct for j=INFO+1,...,N.
- > N: =N+1: other than QZ iteration failed in DHGEQZ.
- =N+2: error return from DTGEVC.
- Further Details
- ===============
- Balancing a matrix pair (A,B) includes, first, permuting rows and
- columns to isolate eigenvalues, second, applying diagonal similarity
- transformation to the rows and columns to make the rows and columns
- as close in norm as possible. The computed reciprocal condition
- numbers correspond to the balanced matrix. Permuting rows and columns
- will not change the condition numbers (in exact arithmetic) but
- diagonal scaling will. For further explanation of balancing, see
- section 4.11.1.2 of LAPACK Users' Guide.
- An approximate error bound on the chordal distance between the i-th
- computed generalized eigenvalue w and the corresponding exact
- eigenvalue lambda is
- chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
- An approximate error bound for the angle between the i-th computed
- eigenvector VL(i) or VR(i) is given by
- EPS * norm(ABNRM, BBNRM) / DIF(i).
- For further explanation of the reciprocal condition numbers RCONDE
- and RCONDV, see section 4.11 of LAPACK User's Guide.
- =====================================================================
- Decode the input arguments
- Parameter adjustments */
- /* Table of constant values */
- static integer c__1 = 1;
- static integer c__0 = 0;
- static doublereal c_b47 = 0.;
- static doublereal c_b48 = 1.;
-
- /* System generated locals */
- integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
- vr_offset, i__1, i__2;
- doublereal d__1, d__2, d__3, d__4;
- /* Builtin functions */
- double sqrt(doublereal);
- /* Local variables */
- static logical pair;
- static doublereal anrm, bnrm;
- static integer ierr, itau;
- static doublereal temp;
- static logical ilvl, ilvr;
- static integer iwrk, iwrk1, i__, j, m;
- extern logical lsame_(char *, char *);
- static integer icols, irows;
- extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
- static integer jc;
- extern /* Subroutine */ int dggbak_(char *, char *, integer *, integer *,
- integer *, doublereal *, doublereal *, integer *, doublereal *,
- integer *, integer *), dggbal_(char *, integer *,
- doublereal *, integer *, doublereal *, integer *, integer *,
- integer *, doublereal *, doublereal *, doublereal *, integer *);
- static integer in;
- extern doublereal dlamch_(char *);
- static integer mm;
- extern doublereal dlange_(char *, integer *, integer *, doublereal *,
- integer *, doublereal *);
- static integer jr;
- extern /* Subroutine */ int dgghrd_(char *, char *, integer *, integer *,
- integer *, doublereal *, integer *, doublereal *, integer *,
- doublereal *, integer *, doublereal *, integer *, integer *), dlascl_(char *, integer *, integer *, doublereal
- *, doublereal *, integer *, integer *, doublereal *, integer *,
- integer *);
- static logical ilascl, ilbscl;
- extern /* Subroutine */ int dgeqrf_(integer *, integer *, doublereal *,
- integer *, doublereal *, doublereal *, integer *, integer *),
- dlacpy_(char *, integer *, integer *, doublereal *, integer *,
- doublereal *, integer *);
- static logical ldumma[1];
- static char chtemp[1];
- static doublereal bignum;
- extern /* Subroutine */ int dhgeqz_(char *, char *, char *, integer *,
- integer *, integer *, doublereal *, integer *, doublereal *,
- integer *, doublereal *, doublereal *, doublereal *, doublereal *,
- integer *, doublereal *, integer *, doublereal *, integer *,
- integer *), dlaset_(char *, integer *,
- integer *, doublereal *, doublereal *, doublereal *, integer *);
- static integer ijobvl;
- extern /* Subroutine */ int dtgevc_(char *, char *, logical *, integer *,
- doublereal *, integer *, doublereal *, integer *, doublereal *,
- integer *, doublereal *, integer *, integer *, integer *,
- doublereal *, integer *), dtgsna_(char *, char *,
- logical *, integer *, doublereal *, integer *, doublereal *,
- integer *, doublereal *, integer *, doublereal *, integer *,
- doublereal *, doublereal *, integer *, integer *, doublereal *,
- integer *, integer *, integer *), xerbla_(char *,
- integer *);
- extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
- integer *, integer *, ftnlen, ftnlen);
- static integer ijobvr;
- static logical wantsb;
- extern /* Subroutine */ int dorgqr_(integer *, integer *, integer *,
- doublereal *, integer *, doublereal *, doublereal *, integer *,
- integer *);
- static doublereal anrmto;
- static logical wantse;
- static doublereal bnrmto;
- extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *,
- integer *, doublereal *, integer *, doublereal *, doublereal *,
- integer *, doublereal *, integer *, integer *);
- static integer minwrk, maxwrk;
- static logical wantsn;
- static doublereal smlnum;
- static logical lquery, wantsv;
- static doublereal eps;
- static logical ilv;
- #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
- #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
- #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
- #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
- a_dim1 = *lda;
- a_offset = 1 + a_dim1 * 1;
- a -= a_offset;
- b_dim1 = *ldb;
- b_offset = 1 + b_dim1 * 1;
- b -= b_offset;
- --alphar;
- --alphai;
- --beta;
- vl_dim1 = *ldvl;
- vl_offset = 1 + vl_dim1 * 1;
- vl -= vl_offset;
- vr_dim1 = *ldvr;
- vr_offset = 1 + vr_dim1 * 1;
- vr -= vr_offset;
- --lscale;
- --rscale;
- --rconde;
- --rcondv;
- --work;
- --iwork;
- --bwork;
- /* Function Body */
- if (lsame_(jobvl, "N")) {
- ijobvl = 1;
- ilvl = FALSE_;
- } else if (lsame_(jobvl, "V")) {
- ijobvl = 2;
- ilvl = TRUE_;
- } else {
- ijobvl = -1;
- ilvl = FALSE_;
- }
- if (lsame_(jobvr, "N")) {
- ijobvr = 1;
- ilvr = FALSE_;
- } else if (lsame_(jobvr, "V")) {
- ijobvr = 2;
- ilvr = TRUE_;
- } else {
- ijobvr = -1;
- ilvr = FALSE_;
- }
- ilv = ilvl || ilvr;
- wantsn = lsame_(sense, "N");
- wantse = lsame_(sense, "E");
- wantsv = lsame_(sense, "V");
- wantsb = lsame_(sense, "B");
- /* Test the input arguments */
- *info = 0;
- lquery = *lwork == -1;
- if (! (lsame_(balanc, "N") || lsame_(balanc, "S") || lsame_(balanc, "P")
- || lsame_(balanc, "B"))) {
- *info = -1;
- } else if (ijobvl <= 0) {
- *info = -2;
- } else if (ijobvr <= 0) {
- *info = -3;
- } else if (! (wantsn || wantse || wantsb || wantsv)) {
- *info = -4;
- } else if (*n < 0) {
- *info = -5;
- } else if (*lda < max(1,*n)) {
- *info = -7;
- } else if (*ldb < max(1,*n)) {
- *info = -9;
- } else if (*ldvl < 1 || ilvl && *ldvl < *n) {
- *info = -14;
- } else if (*ldvr < 1 || ilvr && *ldvr < *n) {
- *info = -16;
- }
- /* Compute workspace
- (Note: Comments in the code beginning "Workspace:" describe the
- minimal amount of workspace needed at that point in the code,
- as well as the preferred amount for good performance.
- NB refers to the optimal block size for the immediately
- following subroutine, as returned by ILAENV. The workspace is
- computed assuming ILO = 1 and IHI = N, the worst case.) */
- minwrk = 1;
- if (*info == 0 && (*lwork >= 1 || lquery)) {
- maxwrk = *n * 5 + *n * ilaenv_(&c__1, "DGEQRF", " ", n, &c__1, n, &
- c__0, (ftnlen)6, (ftnlen)1);
- /* Computing MAX */
- i__1 = 1, i__2 = *n * 6;
- minwrk = max(i__1,i__2);
- if (wantse) {
- /* Computing MAX */
- i__1 = 1, i__2 = *n * 12;
- minwrk = max(i__1,i__2);
- } else if (wantsv || wantsb) {
- minwrk = (*n << 1) * *n + *n * 12 + 16;
- /* Computing MAX */
- i__1 = maxwrk, i__2 = (*n << 1) * *n + *n * 12 + 16;
- maxwrk = max(i__1,i__2);
- }
- work[1] = (doublereal) maxwrk;
- }
- if (*lwork < minwrk && ! lquery) {
- *info = -26;
- }
- if (*info != 0) {
- i__1 = -(*info);
- xerbla_("DGGEVX", &i__1);
- return 0;
- } else if (lquery) {
- return 0;
- }
- /* Quick return if possible */
- if (*n == 0) {
- return 0;
- }
- /* Get machine constants */
- eps = dlamch_("P");
- smlnum = dlamch_("S");
- bignum = 1. / smlnum;
- dlabad_(&smlnum, &bignum);
- smlnum = sqrt(smlnum) / eps;
- bignum = 1. / smlnum;
- /* Scale A if max element outside range [SMLNUM,BIGNUM] */
- anrm = dlange_("M", n, n, &a[a_offset], lda, &work[1]);
- ilascl = FALSE_;
- if (anrm > 0. && anrm < smlnum) {
- anrmto = smlnum;
- ilascl = TRUE_;
- } else if (anrm > bignum) {
- anrmto = bignum;
- ilascl = TRUE_;
- }
- if (ilascl) {
- dlascl_("G", &c__0, &c__0, &anrm, &anrmto, n, n, &a[a_offset], lda, &
- ierr);
- }
- /* Scale B if max element outside range [SMLNUM,BIGNUM] */
- bnrm = dlange_("M", n, n, &b[b_offset], ldb, &work[1]);
- ilbscl = FALSE_;
- if (bnrm > 0. && bnrm < smlnum) {
- bnrmto = smlnum;
- ilbscl = TRUE_;
- } else if (bnrm > bignum) {
- bnrmto = bignum;
- ilbscl = TRUE_;
- }
- if (ilbscl) {
- dlascl_("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, &b[b_offset], ldb, &
- ierr);
- }
- /* Permute and/or balance the matrix pair (A,B)
- (Workspace: need 6*N) */
- dggbal_(balanc, n, &a[a_offset], lda, &b[b_offset], ldb, ilo, ihi, &
- lscale[1], &rscale[1], &work[1], &ierr);
- /* Compute ABNRM and BBNRM */
- *abnrm = dlange_("1", n, n, &a[a_offset], lda, &work[1]);
- if (ilascl) {
- work[1] = *abnrm;
- dlascl_("G", &c__0, &c__0, &anrmto, &anrm, &c__1, &c__1, &work[1], &
- c__1, &ierr);
- *abnrm = work[1];
- }
- *bbnrm = dlange_("1", n, n, &b[b_offset], ldb, &work[1]);
- if (ilbscl) {
- work[1] = *bbnrm;
- dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, &c__1, &c__1, &work[1], &
- c__1, &ierr);
- *bbnrm = work[1];
- }
- /* Reduce B to triangular form (QR decomposition of B)
- (Workspace: need N, prefer N*NB ) */
- irows = *ihi + 1 - *ilo;
- if (ilv || ! wantsn) {
- icols = *n + 1 - *ilo;
- } else {
- icols = irows;
- }
- itau = 1;
- iwrk = itau + irows;
- i__1 = *lwork + 1 - iwrk;
- dgeqrf_(&irows, &icols, &b_ref(*ilo, *ilo), ldb, &work[itau], &work[iwrk],
- &i__1, &ierr);
- /* Apply the orthogonal transformation to A
- (Workspace: need N, prefer N*NB) */
- i__1 = *lwork + 1 - iwrk;
- dormqr_("L", "T", &irows, &icols, &irows, &b_ref(*ilo, *ilo), ldb, &work[
- itau], &a_ref(*ilo, *ilo), lda, &work[iwrk], &i__1, &ierr);
- /* Initialize VL and/or VR
- (Workspace: need N, prefer N*NB) */
- if (ilvl) {
- dlaset_("Full", n, n, &c_b47, &c_b48, &vl[vl_offset], ldvl)
- ;
- i__1 = irows - 1;
- i__2 = irows - 1;
- dlacpy_("L", &i__1, &i__2, &b_ref(*ilo + 1, *ilo), ldb, &vl_ref(*ilo
- + 1, *ilo), ldvl);
- i__1 = *lwork + 1 - iwrk;
- dorgqr_(&irows, &irows, &irows, &vl_ref(*ilo, *ilo), ldvl, &work[itau]
- , &work[iwrk], &i__1, &ierr);
- }
- if (ilvr) {
- dlaset_("Full", n, n, &c_b47, &c_b48, &vr[vr_offset], ldvr)
- ;
- }
- /* Reduce to generalized Hessenberg form
- (Workspace: none needed) */
- if (ilv || ! wantsn) {
- /* Eigenvectors requested -- work on whole matrix. */
- dgghrd_(jobvl, jobvr, n, ilo, ihi, &a[a_offset], lda, &b[b_offset],
- ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, &ierr);
- } else {
- dgghrd_("N", "N", &irows, &c__1, &irows, &a_ref(*ilo, *ilo), lda, &
- b_ref(*ilo, *ilo), ldb, &vl[vl_offset], ldvl, &vr[vr_offset],
- ldvr, &ierr);
- }
- /* Perform QZ algorithm (Compute eigenvalues, and optionally, the
- Schur forms and Schur vectors)
- (Workspace: need N) */
- if (ilv || ! wantsn) {
- *(unsigned char *)chtemp = 'S';
- } else {
- *(unsigned char *)chtemp = 'E';
- }
- dhgeqz_(chtemp, jobvl, jobvr, n, ilo, ihi, &a[a_offset], lda, &b[b_offset]
- , ldb, &alphar[1], &alphai[1], &beta[1], &vl[vl_offset], ldvl, &
- vr[vr_offset], ldvr, &work[1], lwork, &ierr);
- if (ierr != 0) {
- if (ierr > 0 && ierr <= *n) {
- *info = ierr;
- } else if (ierr > *n && ierr <= *n << 1) {
- *info = ierr - *n;
- } else {
- *info = *n + 1;
- }
- goto L130;
- }
- /* Compute Eigenvectors and estimate condition numbers if desired
- (Workspace: DTGEVC: need 6*N
- DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
- need N otherwise ) */
- if (ilv || ! wantsn) {
- if (ilv) {
- if (ilvl) {
- if (ilvr) {
- *(unsigned char *)chtemp = 'B';
- } else {
- *(unsigned char *)chtemp = 'L';
- }
- } else {
- *(unsigned char *)chtemp = 'R';
- }
- dtgevc_(chtemp, "B", ldumma, n, &a[a_offset], lda, &b[b_offset],
- ldb, &vl[vl_offset], ldvl, &vr[vr_offset], ldvr, n, &in, &
- work[1], &ierr);
- if (ierr != 0) {
- *info = *n + 2;
- goto L130;
- }
- }
- if (! wantsn) {
- /* compute eigenvectors (DTGEVC) and estimate condition
- numbers (DTGSNA). Note that the definition of the condition
- number is not invariant under transformation (u,v) to
- (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
- Schur form (S,T), Q and Z are orthogonal matrices. In order
- to avoid using extra 2*N*N workspace, we have to recalculate
- eigenvectors and estimate one condition numbers at a time. */
- pair = FALSE_;
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__) {
- if (pair) {
- pair = FALSE_;
- goto L20;
- }
- mm = 1;
- if (i__ < *n) {
- if (a_ref(i__ + 1, i__) != 0.) {
- pair = TRUE_;
- mm = 2;
- }
- }
- i__2 = *n;
- for (j = 1; j <= i__2; ++j) {
- bwork[j] = FALSE_;
- /* L10: */
- }
- if (mm == 1) {
- bwork[i__] = TRUE_;
- } else if (mm == 2) {
- bwork[i__] = TRUE_;
- bwork[i__ + 1] = TRUE_;
- }
- iwrk = mm * *n + 1;
- iwrk1 = iwrk + mm * *n;
- /* Compute a pair of left and right eigenvectors.
- (compute workspace: need up to 4*N + 6*N) */
- if (wantse || wantsb) {
- dtgevc_("B", "S", &bwork[1], n, &a[a_offset], lda, &b[
- b_offset], ldb, &work[1], n, &work[iwrk], n, &mm,
- &m, &work[iwrk1], &ierr);
- if (ierr != 0) {
- *info = *n + 2;
- goto L130;
- }
- }
- i__2 = *lwork - iwrk1 + 1;
- dtgsna_(sense, "S", &bwork[1], n, &a[a_offset], lda, &b[
- b_offset], ldb, &work[1], n, &work[iwrk], n, &rconde[
- i__], &rcondv[i__], &mm, &m, &work[iwrk1], &i__2, &
- iwork[1], &ierr);
- L20:
- ;
- }
- }
- }
- /* Undo balancing on VL and VR and normalization
- (Workspace: none needed) */
- if (ilvl) {
- dggbak_(balanc, "L", n, ilo, ihi, &lscale[1], &rscale[1], n, &vl[
- vl_offset], ldvl, &ierr);
- i__1 = *n;
- for (jc = 1; jc <= i__1; ++jc) {
- if (alphai[jc] < 0.) {
- goto L70;
- }
- temp = 0.;
- if (alphai[jc] == 0.) {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- /* Computing MAX */
- d__2 = temp, d__3 = (d__1 = vl_ref(jr, jc), abs(d__1));
- temp = max(d__2,d__3);
- /* L30: */
- }
- } else {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- /* Computing MAX */
- d__3 = temp, d__4 = (d__1 = vl_ref(jr, jc), abs(d__1)) + (
- d__2 = vl_ref(jr, jc + 1), abs(d__2));
- temp = max(d__3,d__4);
- /* L40: */
- }
- }
- if (temp < smlnum) {
- goto L70;
- }
- temp = 1. / temp;
- if (alphai[jc] == 0.) {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
- /* L50: */
- }
- } else {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- vl_ref(jr, jc) = vl_ref(jr, jc) * temp;
- vl_ref(jr, jc + 1) = vl_ref(jr, jc + 1) * temp;
- /* L60: */
- }
- }
- L70:
- ;
- }
- }
- if (ilvr) {
- dggbak_(balanc, "R", n, ilo, ihi, &lscale[1], &rscale[1], n, &vr[
- vr_offset], ldvr, &ierr);
- i__1 = *n;
- for (jc = 1; jc <= i__1; ++jc) {
- if (alphai[jc] < 0.) {
- goto L120;
- }
- temp = 0.;
- if (alphai[jc] == 0.) {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- /* Computing MAX */
- d__2 = temp, d__3 = (d__1 = vr_ref(jr, jc), abs(d__1));
- temp = max(d__2,d__3);
- /* L80: */
- }
- } else {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- /* Computing MAX */
- d__3 = temp, d__4 = (d__1 = vr_ref(jr, jc), abs(d__1)) + (
- d__2 = vr_ref(jr, jc + 1), abs(d__2));
- temp = max(d__3,d__4);
- /* L90: */
- }
- }
- if (temp < smlnum) {
- goto L120;
- }
- temp = 1. / temp;
- if (alphai[jc] == 0.) {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
- /* L100: */
- }
- } else {
- i__2 = *n;
- for (jr = 1; jr <= i__2; ++jr) {
- vr_ref(jr, jc) = vr_ref(jr, jc) * temp;
- vr_ref(jr, jc + 1) = vr_ref(jr, jc + 1) * temp;
- /* L110: */
- }
- }
- L120:
- ;
- }
- }
- /* Undo scaling if necessary */
- if (ilascl) {
- dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], n, &
- ierr);
- dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], n, &
- ierr);
- }
- if (ilbscl) {
- dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
- ierr);
- }
- L130:
- work[1] = (doublereal) maxwrk;
- return 0;
- /* End of DGGEVX */
- } /* dggevx_ */
- #undef vr_ref
- #undef vl_ref
- #undef b_ref
- #undef a_ref