dhgeqz.c

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

    z__ -= z_offset;
    --work;

    /* Function Body */
    if (lsame_(job, "E", (ftnlen)1, (ftnlen)1)) {
/*<          ILSCHR = .FALSE. >*/
        ilschr = FALSE_;
/*<          ISCHUR = 1 >*/
        ischur = 1;
/*<       ELSE IF( LSAME( JOB, 'S' ) ) THEN >*/
    } else if (lsame_(job, "S", (ftnlen)1, (ftnlen)1)) {
/*<          ILSCHR = .TRUE. >*/
        ilschr = TRUE_;
/*<          ISCHUR = 2 >*/
        ischur = 2;
/*<       ELSE >*/
    } else {
/*<          ISCHUR = 0 >*/
        ischur = 0;
/*<       END IF >*/
    }

/*<       IF( LSAME( COMPQ, 'N' ) ) THEN >*/
    if (lsame_(compq, "N", (ftnlen)1, (ftnlen)1)) {
/*<          ILQ = .FALSE. >*/
        ilq = FALSE_;
/*<          ICOMPQ = 1 >*/
        icompq = 1;
/*<       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN >*/
    } else if (lsame_(compq, "V", (ftnlen)1, (ftnlen)1)) {
/*<          ILQ = .TRUE. >*/
        ilq = TRUE_;
/*<          ICOMPQ = 2 >*/
        icompq = 2;
/*<       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN >*/
    } else if (lsame_(compq, "I", (ftnlen)1, (ftnlen)1)) {
/*<          ILQ = .TRUE. >*/
        ilq = TRUE_;
/*<          ICOMPQ = 3 >*/
        icompq = 3;
/*<       ELSE >*/
    } else {
/*<          ICOMPQ = 0 >*/
        icompq = 0;
/*<       END IF >*/
    }

/*<       IF( LSAME( COMPZ, 'N' ) ) THEN >*/
    if (lsame_(compz, "N", (ftnlen)1, (ftnlen)1)) {
/*<          ILZ = .FALSE. >*/
        ilz = FALSE_;
/*<          ICOMPZ = 1 >*/
        icompz = 1;
/*<       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN >*/
    } else if (lsame_(compz, "V", (ftnlen)1, (ftnlen)1)) {
/*<          ILZ = .TRUE. >*/
        ilz = TRUE_;
/*<          ICOMPZ = 2 >*/
        icompz = 2;
/*<       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN >*/
    } else if (lsame_(compz, "I", (ftnlen)1, (ftnlen)1)) {
/*<          ILZ = .TRUE. >*/
        ilz = TRUE_;
/*<          ICOMPZ = 3 >*/
        icompz = 3;
/*<       ELSE >*/
    } else {
/*<          ICOMPZ = 0 >*/
        icompz = 0;
/*<       END IF >*/
    }

/*     Check Argument Values */

/*<       INFO = 0 >*/
    *info = 0;
/*<       WORK( 1 ) = MAX( 1, N ) >*/
    work[1] = (doublereal) max(1,*n);
/*<       LQUERY = ( LWORK.EQ.-1 ) >*/
    lquery = *lwork == -1;
/*<       IF( ISCHUR.EQ.0 ) THEN >*/
    if (ischur == 0) {
/*<          INFO = -1 >*/
        *info = -1;
/*<       ELSE IF( ICOMPQ.EQ.0 ) THEN >*/
    } else if (icompq == 0) {
/*<          INFO = -2 >*/
        *info = -2;
/*<       ELSE IF( ICOMPZ.EQ.0 ) THEN >*/
    } else if (icompz == 0) {
/*<          INFO = -3 >*/
        *info = -3;
/*<       ELSE IF( N.LT.0 ) THEN >*/
    } else if (*n < 0) {
/*<          INFO = -4 >*/
        *info = -4;
/*<       ELSE IF( ILO.LT.1 ) THEN >*/
    } else if (*ilo < 1) {
/*<          INFO = -5 >*/
        *info = -5;
/*<       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN >*/
    } else if (*ihi > *n || *ihi < *ilo - 1) {
/*<          INFO = -6 >*/
        *info = -6;
/*<       ELSE IF( LDA.LT.N ) THEN >*/
    } else if (*lda < *n) {
/*<          INFO = -8 >*/
        *info = -8;
/*<       ELSE IF( LDB.LT.N ) THEN >*/
    } else if (*ldb < *n) {
/*<          INFO = -10 >*/
        *info = -10;
/*<       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN >*/
    } else if (*ldq < 1 || (ilq && *ldq < *n)) {
/*<          INFO = -15 >*/
        *info = -15;
/*<       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN >*/
    } else if (*ldz < 1 || (ilz && *ldz < *n)) {
/*<          INFO = -17 >*/
        *info = -17;
/*<       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN >*/
    } else if (*lwork < max(1,*n) && ! lquery) {
/*<          INFO = -19 >*/
        *info = -19;
/*<       END IF >*/
    }
/*<       IF( INFO.NE.0 ) THEN >*/
    if (*info != 0) {
/*<          CALL XERBLA( 'DHGEQZ', -INFO ) >*/
        i__1 = -(*info);
        xerbla_("DHGEQZ", &i__1, (ftnlen)6);
/*<          RETURN >*/
        return 0;
/*<       ELSE IF( LQUERY ) THEN >*/
    } else if (lquery) {
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Quick return if possible */

/*<       IF( N.LE.0 ) THEN >*/
    if (*n <= 0) {
/*<          WORK( 1 ) = DBLE( 1 ) >*/
        work[1] = 1.;
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Initialize Q and Z */

/*<    >*/
    if (icompq == 3) {
        dlaset_("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq, (ftnlen)4);
    }
/*<    >*/
    if (icompz == 3) {
        dlaset_("Full", n, n, &c_b12, &c_b13, &z__[z_offset], ldz, (ftnlen)4);
    }

/*     Machine Constants */

/*<       IN = IHI + 1 - ILO >*/
    in = *ihi + 1 - *ilo;
/*<       SAFMIN = DLAMCH( 'S' ) >*/
    safmin = dlamch_("S", (ftnlen)1);
/*<       SAFMAX = ONE / SAFMIN >*/
    safmax = 1. / safmin;
/*<       ULP = DLAMCH( 'E' )*DLAMCH( 'B' ) >*/
    ulp = dlamch_("E", (ftnlen)1) * dlamch_("B", (ftnlen)1);
/*<       ANORM = DLANHS( 'F', IN, A( ILO, ILO ), LDA, WORK ) >*/
    anorm = dlanhs_("F", &in, &a[*ilo + *ilo * a_dim1], lda, &work[1], (
            ftnlen)1);
/*<       BNORM = DLANHS( 'F', IN, B( ILO, ILO ), LDB, WORK ) >*/
    bnorm = dlanhs_("F", &in, &b[*ilo + *ilo * b_dim1], ldb, &work[1], (
            ftnlen)1);
/*<       ATOL = MAX( SAFMIN, ULP*ANORM ) >*/
/* Computing MAX */
    d__1 = safmin, d__2 = ulp * anorm;
    atol = max(d__1,d__2);
/*<       BTOL = MAX( SAFMIN, ULP*BNORM ) >*/
/* Computing MAX */
    d__1 = safmin, d__2 = ulp * bnorm;
    btol = max(d__1,d__2);
/*<       ASCALE = ONE / MAX( SAFMIN, ANORM ) >*/
    ascale = 1. / max(safmin,anorm);
/*<       BSCALE = ONE / MAX( SAFMIN, BNORM ) >*/
    bscale = 1. / max(safmin,bnorm);

/*     Set Eigenvalues IHI+1:N */

/*<       DO 30 J = IHI + 1, N >*/
    i__1 = *n;
    for (j = *ihi + 1; j <= i__1; ++j) {
/*<          IF( B( J, J ).LT.ZERO ) THEN >*/
        if (b[j + j * b_dim1] < 0.) {
/*<             IF( ILSCHR ) THEN >*/
            if (ilschr) {
/*<                DO 10 JR = 1, J >*/
                i__2 = j;
                for (jr = 1; jr <= i__2; ++jr) {
/*<                   A( JR, J ) = -A( JR, J ) >*/
                    a[jr + j * a_dim1] = -a[jr + j * a_dim1];
/*<                   B( JR, J ) = -B( JR, J ) >*/
                    b[jr + j * b_dim1] = -b[jr + j * b_dim1];
/*<    10          CONTINUE >*/
/* L10: */
                }
/*<             ELSE >*/
            } else {
/*<                A( J, J ) = -A( J, J ) >*/
                a[j + j * a_dim1] = -a[j + j * a_dim1];
/*<                B( J, J ) = -B( J, J ) >*/
                b[j + j * b_dim1] = -b[j + j * b_dim1];
/*<             END IF >*/
            }
/*<             IF( ILZ ) THEN >*/
            if (ilz) {
/*<                DO 20 JR = 1, N >*/
                i__2 = *n;
                for (jr = 1; jr <= i__2; ++jr) {
/*<                   Z( JR, J ) = -Z( JR, J ) >*/
                    z__[jr + j * z_dim1] = -z__[jr + j * z_dim1];
/*<    20          CONTINUE >*/
/* L20: */
                }
/*<             END IF >*/
            }
/*<          END IF >*/
        }
/*<          ALPHAR( J ) = A( J, J ) >*/
        alphar[j] = a[j + j * a_dim1];
/*<          ALPHAI( J ) = ZERO >*/
        alphai[j] = 0.;
/*<          BETA( J ) = B( J, J ) >*/
        beta[j] = b[j + j * b_dim1];
/*<    30 CONTINUE >*/
/* L30: */
    }

/*     If IHI < ILO, skip QZ steps */

/*<    >*/
    if (*ihi < *ilo) {
        goto L380;
    }

/*     MAIN QZ ITERATION LOOP */

/*     Initialize dynamic indices */

/*     Eigenvalues ILAST+1:N have been found. */
/*        Column operations modify rows IFRSTM:whatever. */
/*        Row operations modify columns whatever:ILASTM. */

/*     If only eigenvalues are being computed, then */
/*        IFRSTM is the row of the last splitting row above row ILAST; */
/*        this is always at least ILO. */
/*     IITER counts iterations since the last eigenvalue was found, */
/*        to tell when to use an extraordinary shift. */
/*     MAXIT is the maximum number of QZ sweeps allowed. */

/*<       ILAST = IHI >*/
    ilast = *ihi;
/*<       IF( ILSCHR ) THEN >*/
    if (ilschr) {
/*<          IFRSTM = 1 >*/
        ifrstm = 1;
/*<          ILASTM = N >*/
        ilastm = *n;
/*<       ELSE >*/
    } else {
/*<          IFRSTM = ILO >*/
        ifrstm = *ilo;
/*<          ILASTM = IHI >*/
        ilastm = *ihi;
/*<       END IF >*/
    }
/*<       IITER = 0 >*/
    iiter = 0;
/*<       ESHIFT = ZERO >*/
    eshift = 0.;
/*<       MAXIT = 30*( IHI-ILO+1 ) >*/
    maxit = (*ihi - *ilo + 1) * 30;

/*<       DO 360 JITER = 1, MAXIT >*/
    i__1 = maxit;
    for (jiter = 1; jiter <= i__1; ++jiter) {

/*        Split the matrix if possible. */

/*        Two tests: */
/*           1: A(j,j-1)=0  or  j=ILO */
/*           2: B(j,j)=0 */

/*<          IF( ILAST.EQ.ILO ) THEN >*/
        if (ilast == *ilo) {

/*           Special case: j=ILAST */

/*<             GO TO 80 >*/
            goto L80;
/*<          ELSE >*/
        } else {
/*<             IF( ABS( A( ILAST, ILAST-1 ) ).LE.ATOL ) THEN >*/
            if ((d__1 = a[ilast + (ilast - 1) * a_dim1], abs(d__1)) <= atol) {
/*<                A( ILAST, ILAST-1 ) = ZERO >*/
                a[ilast + (ilast - 1) * a_dim1] = 0.;
/*<                GO TO 80 >*/
                goto L80;
/*<             END IF >*/
            }
/*<          END IF >*/
        }

/*<          IF( ABS( B( ILAST, ILAST ) ).LE.BTOL ) THEN >*/
        if ((d__1 = b[ilast + ilast * b_dim1], abs(d__1)) <= btol) {
/*<             B( ILAST, ILAST ) = ZERO >*/
            b[ilast + ilast * b_dim1] = 0.;
/*<             GO TO 70 >*/
            goto L70;
/*<          END IF >*/
        }

/*        General case: j<ILAST */

/*<          DO 60 J = ILAST - 1, ILO, -1 >*/
        i__2 = *ilo;
        for (j = ilast - 1; j >= i__2; --j) {

/*           Test 1: for A(j,j-1)=0 or j=ILO */

/*<             IF( J.EQ.ILO ) THEN >*/
            if (j == *ilo) {
/*<                ILAZRO = .TRUE. >*/
                ilazro = TRUE_;
/*<             ELSE >*/
            } else {
/*<                IF( ABS( A( J, J-1 ) ).LE.ATOL ) THEN >*/
                if ((d__1 = a[j + (j - 1) * a_dim1], abs(d__1)) <= atol) {
/*<                   A( J, J-1 ) = ZERO >*/
                    a[j + (j - 1) * a_dim1] = 0.;
/*<                   ILAZRO = .TRUE. >*/
                    ilazro = TRUE_;