dhgeqz.c

来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 1,728 行 · 第 1/5 页

/* lapack/double/dhgeqz.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 doublereal c_b12 = 0.;
static doublereal c_b13 = 1.;
static integer c__1 = 1;
static integer c__3 = 3;

/*<    >*/
/* Subroutine */ int dhgeqz_(char *job, char *compq, char *compz, integer *n, 
        integer *ilo, integer *ihi, doublereal *a, integer *lda, doublereal *
        b, integer *ldb, doublereal *alphar, doublereal *alphai, doublereal *
        beta, doublereal *q, integer *ldq, doublereal *z__, integer *ldz, 
        doublereal *work, integer *lwork, integer *info, ftnlen job_len, 
        ftnlen compq_len, ftnlen compz_len)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 
            z_offset, i__1, i__2, i__3, i__4;
    doublereal d__1, d__2, d__3, d__4;

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

    /* Local variables */
    doublereal c__;
    integer j;
    doublereal s, t, v[3], s1, s2, u1, u2, a11, a12, a21, a22, b11, b22, c12, 
            c21;
    integer jc;
    doublereal an, bn, cl, cq, cr;
    integer in;
    doublereal u12, w11, w12, w21;
    integer jr;
    doublereal cz, w22, sl, wi, sr, vs, wr, b1a, b2a, a1i, a2i, b1i, b2i, a1r,
             a2r, b1r, b2r, wr2, ad11, ad12, ad21, ad22, c11i, c22i;
    integer jch;
    doublereal c11r, c22r, u12l;
    logical ilq=0;
    doublereal tau, sqi;
    logical ilz=0;
    doublereal ulp, sqr, szi, szr, ad11l, ad12l, ad21l, ad22l, ad32l, wabs, 
            atol, btol, temp;
    extern /* Subroutine */ int drot_(integer *, doublereal *, integer *, 
            doublereal *, integer *, doublereal *, doublereal *), dlag2_(
            doublereal *, integer *, doublereal *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *);
    doublereal temp2, s1inv, scale;
    extern logical lsame_(char *, char *, ftnlen, ftnlen);
    integer iiter, ilast, jiter;
    doublereal anorm, bnorm;
    integer maxit;
    doublereal tempi, tempr;
    extern doublereal dlapy2_(doublereal *, doublereal *), dlapy3_(doublereal 
            *, doublereal *, doublereal *);
    extern /* Subroutine */ int dlasv2_(doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *);
    logical ilazr2;
    doublereal ascale, bscale;
    extern doublereal dlamch_(char *, ftnlen);
    extern /* Subroutine */ int dlarfg_(integer *, doublereal *, doublereal *,
             integer *, doublereal *);
    extern doublereal dlanhs_(char *, integer *, doublereal *, integer *, 
            doublereal *, ftnlen);
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
            doublereal *, doublereal *, doublereal *, integer *, ftnlen);
    doublereal safmin;
    extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *);
    doublereal safmax;
    extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
    doublereal eshift;
    logical ilschr=0;
    integer icompq, ilastm, ischur;
    logical ilazro;
    integer icompz, ifirst, ifrstm, istart;
    logical ilpivt, lquery;
    (void)job_len;
    (void)compq_len;
    (void)compz_len;

/*  -- LAPACK 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          COMPQ, COMPZ, JOB >*/
/*<       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, LWORK, N >*/
/*     .. */
/*     .. Array Arguments .. */
/*<    >*/
/*     .. */

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

/*  DHGEQZ implements a single-/double-shift version of the QZ method for */
/*  finding the generalized eigenvalues */

/*  w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j)   of the equation */

/*       det( A - w(i) B ) = 0 */

/*  In addition, the pair A,B may be reduced to generalized Schur form: */
/*  B is upper triangular, and A is block upper triangular, where the */
/*  diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having */
/*  complex generalized eigenvalues (see the description of the argument */
/*  JOB.) */

/*  If JOB='S', then the pair (A,B) is simultaneously reduced to Schur */
/*  form by applying one orthogonal transformation (usually called Q) on */
/*  the left and another (usually called Z) on the right.  The 2-by-2 */
/*  upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks */
/*  of A will be reduced to positive diagonal matrices.  (I.e., */
/*  if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and */
/*  B(j+1,j+1) will be positive.) */

/*  If JOB='E', then at each iteration, the same transformations */
/*  are computed, but they are only applied to those parts of A and B */
/*  which are needed to compute ALPHAR, ALPHAI, and BETAR. */

/*  If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal */
/*  transformations used to reduce (A,B) are accumulated into the arrays */
/*  Q and Z s.t.: */

/*       Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)* */
/*       Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)* */

/*  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix */
/*       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), */
/*       pp. 241--256. */

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

/*  JOB     (input) CHARACTER*1 */
/*          = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will */
/*                 not necessarily be put into generalized Schur form. */
/*          = 'S': put A and B into generalized Schur form, as well */
/*                 as computing ALPHAR, ALPHAI, and BETA. */

/*  COMPQ   (input) CHARACTER*1 */
/*          = 'N': do not modify Q. */
/*          = 'V': multiply the array Q on the right by the transpose of */
/*                 the orthogonal transformation that is applied to the */
/*                 left side of A and B to reduce them to Schur form. */
/*          = 'I': like COMPQ='V', except that Q will be initialized to */
/*                 the identity first. */

/*  COMPZ   (input) CHARACTER*1 */
/*          = 'N': do not modify Z. */
/*          = 'V': multiply the array Z on the right by the orthogonal */
/*                 transformation that is applied to the right side of */
/*                 A and B to reduce them to Schur form. */
/*          = 'I': like COMPZ='V', except that Z will be initialized to */
/*                 the identity first. */

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

/*  ILO     (input) INTEGER */
/*  IHI     (input) INTEGER */
/*          It is assumed that A is already upper triangular in rows and */
/*          columns 1:ILO-1 and IHI+1:N. */
/*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */

/*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
/*          On entry, the N-by-N upper Hessenberg matrix A.  Elements */
/*          below the subdiagonal must be zero. */
/*          If JOB='S', then on exit A and B will have been */
/*             simultaneously reduced to generalized Schur form. */
/*          If JOB='E', then on exit A will have been destroyed. */
/*             The diagonal blocks will be correct, but the off-diagonal */
/*             portion will be meaningless. */

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

/*  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N) */
/*          On entry, the N-by-N upper triangular matrix B.  Elements */
/*          below the diagonal must be zero.  2-by-2 blocks in B */
/*          corresponding to 2-by-2 blocks in A will be reduced to */
/*          positive diagonal form.  (I.e., if A(j+1,j) is non-zero, */
/*          then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be */
/*          positive.) */
/*          If JOB='S', then on exit A and B will have been */
/*             simultaneously reduced to Schur form. */
/*          If JOB='E', then on exit B will have been destroyed. */
/*             Elements corresponding to diagonal blocks of A will be */
/*             correct, but the off-diagonal portion will be meaningless. */

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

/*  ALPHAR  (output) DOUBLE PRECISION array, dimension (N) */
/*          ALPHAR(1:N) will be set to real parts of the diagonal */
/*          elements of A that would result from reducing A and B to */
/*          Schur form and then further reducing them both to triangular */
/*          form using unitary transformations s.t. the diagonal of B */
/*          was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block */
/*          (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j). */
/*          Note that the (real or complex) values */
/*          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the */
/*          generalized eigenvalues of the matrix pencil A - wB. */

/*  ALPHAI  (output) DOUBLE PRECISION array, dimension (N) */
/*          ALPHAI(1:N) will be set to imaginary parts of the diagonal */
/*          elements of A that would result from reducing A and B to */
/*          Schur form and then further reducing them both to triangular */
/*          form using unitary transformations s.t. the diagonal of B */
/*          was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block */
/*          (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0. */
/*          Note that the (real or complex) values */
/*          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the */
/*          generalized eigenvalues of the matrix pencil A - wB. */

/*  BETA    (output) DOUBLE PRECISION array, dimension (N) */
/*          BETA(1:N) will be set to the (real) diagonal elements of B */
/*          that would result from reducing A and B to Schur form and */
/*          then further reducing them both to triangular form using */
/*          unitary transformations s.t. the diagonal of B was */
/*          non-negative real.  Thus, if A(j,j) is in a 1-by-1 block */
/*          (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j). */
/*          Note that the (real or complex) values */
/*          (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the */
/*          generalized eigenvalues of the matrix pencil A - wB. */
/*          (Note that BETA(1:N) will always be non-negative, and no */
/*          BETAI is necessary.) */

/*  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N) */
/*          If COMPQ='N', then Q will not be referenced. */
/*          If COMPQ='V' or 'I', then the transpose of the orthogonal */
/*             transformations which are applied to A and B on the left */
/*             will be applied to the array Q on the right. */

/*  LDQ     (input) INTEGER */
/*          The leading dimension of the array Q.  LDQ >= 1. */
/*          If COMPQ='V' or 'I', then LDQ >= N. */

/*  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N) */
/*          If COMPZ='N', then Z will not be referenced. */
/*          If COMPZ='V' or 'I', then the orthogonal transformations */
/*             which are applied to A and B on the right will be applied */
/*             to the array Z on the right. */

/*  LDZ     (input) INTEGER */
/*          The leading dimension of the array Z.  LDZ >= 1. */
/*          If COMPZ='V' or 'I', then LDZ >= 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 >= max(1,N). */

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

/*  INFO    (output) INTEGER */
/*          = 0: successful exit */
/*          < 0: if INFO = -i, the i-th argument had an illegal value */
/*          = 1,...,N: the QZ iteration did not converge.  (A,B) is not */
/*                     in Schur form, but ALPHAR(i), ALPHAI(i), and */
/*                     BETA(i), i=INFO+1,...,N should be correct. */
/*          = N+1,...,2*N: the shift calculation failed.  (A,B) is not */
/*                     in Schur form, but ALPHAR(i), ALPHAI(i), and */
/*                     BETA(i), i=INFO-N+1,...,N should be correct. */
/*          > 2*N:     various "impossible" errors. */

/*  Further Details */
/*  =============== */

/*  Iteration counters: */

/*  JITER  -- counts iterations. */
/*  IITER  -- counts iterations run since ILAST was last */
/*            changed.  This is therefore reset only when a 1-by-1 or */
/*            2-by-2 block deflates off the bottom. */

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

/*     .. Parameters .. */
/*    $                     SAFETY = 1.0E+0 ) */
/*<       DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY >*/
/*<    >*/
/*     .. */
/*     .. Local Scalars .. */
/*<    >*/
/*<    >*/
/*<    >*/
/*     .. */
/*     .. Local Arrays .. */
/*<       DOUBLE PRECISION   V( 3 ) >*/
/*     .. */
/*     .. External Functions .. */
/*<       LOGICAL            LSAME >*/
/*<       DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2, DLAPY3 >*/
/*<       EXTERNAL           LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3 >*/
/*     .. */
/*     .. External Subroutines .. */
/*<    >*/
/*     .. */
/*     .. Intrinsic Functions .. */
/*<       INTRINSIC          ABS, DBLE, MAX, MIN, SQRT >*/
/*     .. */
/*     .. Executable Statements .. */

/*     Decode JOB, COMPQ, COMPZ */

/*<       IF( LSAME( JOB, 'E' ) ) THEN >*/
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;
    --alphar;
    --alphai;
    --beta;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1;
    q -= q_offset;
    z_dim1 = *ldz;
    z_offset = 1 + z_dim1;