dgges.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/* lapack/double/dgges.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Table of constant values */

static integer c__1 = 1;
static integer c__0 = 0;
static integer c_n1 = -1;
static doublereal c_b33 = 0.;
static doublereal c_b34 = 1.;

/*<    >*/
/* Subroutine */ int dgges_(char *jobvsl, char *jobvsr, char *sort,
        logical (*delctg)(doublereal*,doublereal*,doublereal*),
        integer *n, doublereal *a, integer *lda, doublereal *b, 
        integer *ldb, integer *sdim, doublereal *alphar, doublereal *alphai, 
        doublereal *beta, doublereal *vsl, integer *ldvsl, doublereal *vsr, 
        integer *ldvsr, doublereal *work, integer *lwork, logical *bwork, 
        integer *info, ftnlen jobvsl_len, ftnlen jobvsr_len, ftnlen sort_len)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, vsl_dim1, vsl_offset, 
            vsr_dim1, vsr_offset, i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    integer i__, ip;
    doublereal dif[2];
    integer ihi, ilo;
    doublereal eps, anrm, bnrm;
    integer idum[1], ierr, itau, iwrk;
    doublereal pvsl, pvsr;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    integer ileft, icols;
    logical cursl, ilvsl, ilvsr;
    integer irows;
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *), dggbak_(
            char *, char *, integer *, integer *, integer *, doublereal *, 
            doublereal *, integer *, doublereal *, integer *, integer *, 
            ftnlen, ftnlen), dggbal_(char *, integer *, doublereal *, integer 
            *, doublereal *, integer *, integer *, integer *, doublereal *, 
            doublereal *, doublereal *, integer *, ftnlen);
    logical lst2sl;
    extern doublereal dlamch_(char *, ftnlen), dlange_(char *, integer *, 
            integer *, doublereal *, integer *, doublereal *, ftnlen);
    extern /* Subroutine */ int dgghrd_(char *, char *, integer *, integer *, 
            integer *, doublereal *, integer *, doublereal *, integer *, 
            doublereal *, integer *, doublereal *, integer *, integer *, 
            ftnlen, ftnlen), dlascl_(char *, integer *, integer *, doublereal 
            *, doublereal *, integer *, integer *, doublereal *, integer *, 
            integer *, ftnlen);
    logical ilascl, ilbscl;
    extern /* Subroutine */ int dgeqrf_(integer *, integer *, doublereal *, 
            integer *, doublereal *, doublereal *, integer *, integer *), 
            dlacpy_(char *, integer *, integer *, doublereal *, integer *, 
            doublereal *, integer *, ftnlen);
    doublereal safmin;
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
            doublereal *, doublereal *, doublereal *, integer *, ftnlen);
    doublereal safmax;
    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
    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 *, ftnlen, ftnlen, ftnlen);
    extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
            integer *, integer *, ftnlen, ftnlen);
    extern /* Subroutine */ int dtgsen_(integer *, logical *, logical *, 
            logical *, integer *, doublereal *, integer *, doublereal *, 
            integer *, doublereal *, doublereal *, doublereal *, doublereal *,
             integer *, doublereal *, integer *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, integer *, integer *, 
            integer *, integer *);
    integer ijobvl, iright, ijobvr;
    extern /* Subroutine */ int dorgqr_(integer *, integer *, integer *, 
            doublereal *, integer *, doublereal *, doublereal *, integer *, 
            integer *);
    doublereal anrmto, bnrmto;
    logical lastsl;
    extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, 
            integer *, doublereal *, integer *, doublereal *, doublereal *, 
            integer *, doublereal *, integer *, integer *, ftnlen, ftnlen);
    integer minwrk, maxwrk=0;
    doublereal smlnum;
    logical wantst, lquery;
    (void)jobvsl_len;
    (void)jobvsr_len;
    (void)sort_len;

/*  -- 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 */

/*     .. Scalar Arguments .. */
/*<       CHARACTER          JOBVSL, JOBVSR, SORT >*/
/*<       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM >*/
/*     .. */
/*     .. Array Arguments .. */
/*<       LOGICAL            BWORK( * ) >*/
/*<    >*/
/*     .. */
/*     .. Function Arguments .. */
/*<       LOGICAL            DELCTG >*/
/*<       EXTERNAL           DELCTG >*/
/*     .. */

/*  Purpose */
/*  ======= */

/*  DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B), */
/*  the generalized eigenvalues, the generalized real Schur form (S,T), */
/*  optionally, the left and/or right matrices of Schur vectors (VSL and */
/*  VSR). This gives the generalized Schur factorization */

/*           (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T ) */

/*  Optionally, it also orders the eigenvalues so that a selected cluster */
/*  of eigenvalues appears in the leading diagonal blocks of the upper */
/*  quasi-triangular matrix S and the upper triangular matrix T.The */
/*  leading columns of VSL and VSR then form an orthonormal basis for the */
/*  corresponding left and right eigenspaces (deflating subspaces). */

/*  (If only the generalized eigenvalues are needed, use the driver */
/*  DGGEV instead, which is faster.) */

/*  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w */
/*  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is */
/*  usually represented as the pair (alpha,beta), as there is a */
/*  reasonable interpretation for beta=0 or both being zero. */

/*  A pair of matrices (S,T) is in generalized real Schur form if T is */
/*  upper triangular with non-negative diagonal and S is block upper */
/*  triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond */
/*  to real generalized eigenvalues, while 2-by-2 blocks of S will be */
/*  "standardized" by making the corresponding elements of T have the */
/*  form: */
/*          [  a  0  ] */
/*          [  0  b  ] */

/*  and the pair of corresponding 2-by-2 blocks in S and T will have a */
/*  complex conjugate pair of generalized eigenvalues. */


/*  Arguments */
/*  ========= */

/*  JOBVSL  (input) CHARACTER*1 */
/*          = 'N':  do not compute the left Schur vectors; */
/*          = 'V':  compute the left Schur vectors. */

/*  JOBVSR  (input) CHARACTER*1 */
/*          = 'N':  do not compute the right Schur vectors; */
/*          = 'V':  compute the right Schur vectors. */

/*  SORT    (input) CHARACTER*1 */
/*          Specifies whether or not to order the eigenvalues on the */
/*          diagonal of the generalized Schur form. */
/*          = 'N':  Eigenvalues are not ordered; */
/*          = 'S':  Eigenvalues are ordered (see DELZTG); */

/*  DELZTG  (input) LOGICAL FUNCTION of three DOUBLE PRECISION arguments */
/*          DELZTG must be declared EXTERNAL in the calling subroutine. */
/*          If SORT = 'N', DELZTG is not referenced. */
/*          If SORT = 'S', DELZTG is used to select eigenvalues to sort */
/*          to the top left of the Schur form. */
/*          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if */
/*          DELZTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either */
/*          one of a complex conjugate pair of eigenvalues is selected, */
/*          then both complex eigenvalues are selected. */

/*          Note that in the ill-conditioned case, a selected complex */
/*          eigenvalue may no longer satisfy DELZTG(ALPHAR(j),ALPHAI(j), */
/*          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2 */
/*          in this case. */

/*  N       (input) INTEGER */
/*          The order of the matrices A, B, VSL, and VSR.  N >= 0. */

/*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
/*          On entry, the first of the pair of matrices. */
/*          On exit, A has been overwritten by its generalized Schur */
/*          form S. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of A.  LDA >= max(1,N). */

/*  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N) */
/*          On entry, the second of the pair of matrices. */
/*          On exit, B has been overwritten by its generalized Schur */
/*          form T. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of B.  LDB >= max(1,N). */

/*  SDIM    (output) INTEGER */
/*          If SORT = 'N', SDIM = 0. */
/*          If SORT = 'S', SDIM = number of eigenvalues (after sorting) */
/*          for which DELZTG is true.  (Complex conjugate pairs for which */
/*          DELZTG is true for either eigenvalue count as 2.) */

/*  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.  ALPHAR(j) + ALPHAI(j)*i, */
/*          and  BETA(j),j=1,...,N are the diagonals of the complex Schur */
/*          form (S,T) that would result if the 2-by-2 diagonal blocks of */
/*          the real Schur form of (A,B) were further reduced to */
/*          triangular form using 2-by-2 complex unitary transformations. */
/*          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. */
/*          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). */

/*  VSL     (output) DOUBLE PRECISION array, dimension (LDVSL,N) */
/*          If JOBVSL = 'V', VSL will contain the left Schur vectors. */
/*          Not referenced if JOBVSL = 'N'. */

/*  LDVSL   (input) INTEGER */
/*          The leading dimension of the matrix VSL. LDVSL >=1, and */
/*          if JOBVSL = 'V', LDVSL >= N. */

/*  VSR     (output) DOUBLE PRECISION array, dimension (LDVSR,N) */
/*          If JOBVSR = 'V', VSR will contain the right Schur vectors. */
/*          Not referenced if JOBVSR = 'N'. */

/*  LDVSR   (input) INTEGER */
/*          The leading dimension of the matrix VSR. LDVSR >= 1, and */
/*          if JOBVSR = 'V', LDVSR >= N. */

/*  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 >= 8*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. */

/*  BWORK   (workspace) LOGICAL array, dimension (N) */
/*          Not referenced if SORT = 'N'. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
/*          = 1,...,N: */
/*                The QZ iteration failed.  (A,B) are not in Schur */
/*                form, 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: after reordering, roundoff changed values of */
/*                      some complex eigenvalues so that leading */
/*                      eigenvalues in the Generalized Schur form no */
/*                      longer satisfy DELZTG=.TRUE.  This could also */
/*                      be caused due to scaling. */
/*                =N+3: reordering failed in DTGSEN. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*<       DOUBLE PRECISION   ZERO, ONE >*/
/*<       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 ) >*/
/*     .. */
/*     .. Local Scalars .. */
/*<    >*/
/*<    >*/
/*<    >*/
/*     .. */
/*     .. Local Arrays .. */
/*<       INTEGER            IDUM( 1 ) >*/
/*<       DOUBLE PRECISION   DIF( 2 ) >*/
/*     .. */
/*     .. External Subroutines .. */
/*<    >*/
/*     .. */
/*     .. External Functions .. */