dgges.c

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

/*<    >*/
    if (ilvsr) {
        dlaset_("Full", n, n, &c_b33, &c_b34, &vsr[vsr_offset], ldvsr, (
                ftnlen)4);
    }

/*     Reduce to generalized Hessenberg form */
/*     (Workspace: none needed) */

/*<    >*/
    dgghrd_(jobvsl, jobvsr, n, &ilo, &ihi, &a[a_offset], lda, &b[b_offset], 
            ldb, &vsl[vsl_offset], ldvsl, &vsr[vsr_offset], ldvsr, &ierr, (
            ftnlen)1, (ftnlen)1);

/*     Perform QZ algorithm, computing Schur vectors if desired */
/*     (Workspace: need N) */

/*<       IWRK = ITAU >*/
    iwrk = itau;
/*<    >*/
    i__1 = *lwork + 1 - iwrk;
    dhgeqz_("S", jobvsl, jobvsr, n, &ilo, &ihi, &a[a_offset], lda, &b[
            b_offset], ldb, &alphar[1], &alphai[1], &beta[1], &vsl[vsl_offset]
            , ldvsl, &vsr[vsr_offset], ldvsr, &work[iwrk], &i__1, &ierr, (
            ftnlen)1, (ftnlen)1, (ftnlen)1);
/*<       IF( IERR.NE.0 ) THEN >*/
    if (ierr != 0) {
/*<          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN >*/
        if (ierr > 0 && ierr <= *n) {
/*<             INFO = IERR >*/
            *info = ierr;
/*<          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN >*/
        } else if (ierr > *n && ierr <= *n << 1) {
/*<             INFO = IERR - N >*/
            *info = ierr - *n;
/*<          ELSE >*/
        } else {
/*<             INFO = N + 1 >*/
            *info = *n + 1;
/*<          END IF >*/
        }
/*<          GO TO 50 >*/
        goto L50;
/*<       END IF >*/
    }

/*     Sort eigenvalues ALPHA/BETA if desired */
/*     (Workspace: need 4*N+16 ) */

/*<       SDIM = 0 >*/
    *sdim = 0;
/*<       IF( WANTST ) THEN >*/
    if (wantst) {

/*        Undo scaling on eigenvalues before DELZTGing */

/*<          IF( ILASCL ) THEN >*/
        if (ilascl) {
/*<    >*/
            dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], 
                    n, &ierr, (ftnlen)1);
/*<    >*/
            dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], 
                    n, &ierr, (ftnlen)1);
/*<          END IF >*/
        }
/*<    >*/
        if (ilbscl) {
            dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, 
                    &ierr, (ftnlen)1);
        }

/*        Select eigenvalues */

/*<          DO 10 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             BWORK( I ) = DELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) ) >*/
            bwork[i__] = (*delctg)(&alphar[i__], &alphai[i__], &beta[i__]);
/*<    10    CONTINUE >*/
/* L10: */
        }

/*<    >*/
        i__1 = *lwork - iwrk + 1;
        dtgsen_(&c__0, &ilvsl, &ilvsr, &bwork[1], n, &a[a_offset], lda, &b[
                b_offset], ldb, &alphar[1], &alphai[1], &beta[1], &vsl[
                vsl_offset], ldvsl, &vsr[vsr_offset], ldvsr, sdim, &pvsl, &
                pvsr, dif, &work[iwrk], &i__1, idum, &c__1, &ierr);
/*<    >*/
        if (ierr == 1) {
            *info = *n + 3;
        }

/*<       END IF >*/
    }

/*     Apply back-permutation to VSL and VSR */
/*     (Workspace: none needed) */

/*<    >*/
    if (ilvsl) {
        dggbak_("P", "L", n, &ilo, &ihi, &work[ileft], &work[iright], n, &vsl[
                vsl_offset], ldvsl, &ierr, (ftnlen)1, (ftnlen)1);
    }

/*<    >*/
    if (ilvsr) {
        dggbak_("P", "R", n, &ilo, &ihi, &work[ileft], &work[iright], n, &vsr[
                vsr_offset], ldvsr, &ierr, (ftnlen)1, (ftnlen)1);
    }

/*     Check if unscaling would cause over/underflow, if so, rescale */
/*     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of */
/*     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I) */

/*<       IF( ILASCL ) THEN >*/
    if (ilascl) {
/*<          DO 20 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             IF( ALPHAI( I ).NE.ZERO ) THEN >*/
            if (alphai[i__] != 0.) {
/*<    >*/
                if (alphar[i__] / safmax > anrmto / anrm || safmin / alphar[
                        i__] > anrm / anrmto) {
/*<                   WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) ) >*/
                    work[1] = (d__1 = a[i__ + i__ * a_dim1] / alphar[i__], 
                            abs(d__1));
/*<                   BETA( I ) = BETA( I )*WORK( 1 ) >*/
                    beta[i__] *= work[1];
/*<                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) >*/
                    alphar[i__] *= work[1];
/*<                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) >*/
                    alphai[i__] *= work[1];
/*<    >*/
                } else if (alphai[i__] / safmax > anrmto / anrm || safmin / 
                        alphai[i__] > anrm / anrmto) {
/*<                   WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) ) >*/
                    work[1] = (d__1 = a[i__ + (i__ + 1) * a_dim1] / alphai[
                            i__], abs(d__1));
/*<                   BETA( I ) = BETA( I )*WORK( 1 ) >*/
                    beta[i__] *= work[1];
/*<                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) >*/
                    alphar[i__] *= work[1];
/*<                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) >*/
                    alphai[i__] *= work[1];
/*<                END IF >*/
                }
/*<             END IF >*/
            }
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<       END IF >*/
    }

/*<       IF( ILBSCL ) THEN >*/
    if (ilbscl) {
/*<          DO 30 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             IF( ALPHAI( I ).NE.ZERO ) THEN >*/
            if (alphai[i__] != 0.) {
/*<    >*/
                if (beta[i__] / safmax > bnrmto / bnrm || safmin / beta[i__] 
                        > bnrm / bnrmto) {
/*<                   WORK( 1 ) = ABS( B( I, I ) / BETA( I ) ) >*/
                    work[1] = (d__1 = b[i__ + i__ * b_dim1] / beta[i__], abs(
                            d__1));
/*<                   BETA( I ) = BETA( I )*WORK( 1 ) >*/
                    beta[i__] *= work[1];
/*<                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 ) >*/
                    alphar[i__] *= work[1];
/*<                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 ) >*/
                    alphai[i__] *= work[1];
/*<                END IF >*/
                }
/*<             END IF >*/
            }
/*<    30    CONTINUE >*/
/* L30: */
        }
/*<       END IF >*/
    }

/*     Undo scaling */

/*<       IF( ILASCL ) THEN >*/
    if (ilascl) {
/*<          CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) >*/
        dlascl_("H", &c__0, &c__0, &anrmto, &anrm, n, n, &a[a_offset], lda, &
                ierr, (ftnlen)1);
/*<          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR ) >*/
        dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphar[1], n, &
                ierr, (ftnlen)1);
/*<          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR ) >*/
        dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, &alphai[1], n, &
                ierr, (ftnlen)1);
/*<       END IF >*/
    }

/*<       IF( ILBSCL ) THEN >*/
    if (ilbscl) {
/*<          CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) >*/
        dlascl_("U", &c__0, &c__0, &bnrmto, &bnrm, n, n, &b[b_offset], ldb, &
                ierr, (ftnlen)1);
/*<          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) >*/
        dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, &beta[1], n, &
                ierr, (ftnlen)1);
/*<       END IF >*/
    }

/*<       IF( WANTST ) THEN >*/
    if (wantst) {

/*        Check if reordering is correct */

/*<          LASTSL = .TRUE. >*/
        lastsl = TRUE_;
/*<          LST2SL = .TRUE. >*/
        lst2sl = TRUE_;
/*<          SDIM = 0 >*/
        *sdim = 0;
/*<          IP = 0 >*/
        ip = 0;
/*<          DO 40 I = 1, N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             CURSL = DELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) ) >*/
            cursl = (*delctg)(&alphar[i__], &alphai[i__], &beta[i__]);
/*<             IF( ALPHAI( I ).EQ.ZERO ) THEN >*/
            if (alphai[i__] == 0.) {
/*<    >*/
                if (cursl) {
                    ++(*sdim);
                }
/*<                IP = 0 >*/
                ip = 0;
/*<    >*/
                if (cursl && ! lastsl) {
                    *info = *n + 2;
                }
/*<             ELSE >*/
            } else {
/*<                IF( IP.EQ.1 ) THEN >*/
                if (ip == 1) {

/*                 Last eigenvalue of conjugate pair */

/*<                   CURSL = CURSL .OR. LASTSL >*/
                    cursl = cursl || lastsl;
/*<                   LASTSL = CURSL >*/
                    lastsl = cursl;
/*<    >*/
                    if (cursl) {
                        *sdim += 2;
                    }
/*<                   IP = -1 >*/
                    ip = -1;
/*<    >*/
                    if (cursl && ! lst2sl) {
                        *info = *n + 2;
                    }
/*<                ELSE >*/
                } else {

/*                 First eigenvalue of conjugate pair */

/*<                   IP = 1 >*/
                    ip = 1;
/*<                END IF >*/
                }
/*<             END IF >*/
            }
/*<             LST2SL = LASTSL >*/
            lst2sl = lastsl;
/*<             LASTSL = CURSL >*/
            lastsl = cursl;
/*<    40    CONTINUE >*/
/* L40: */
        }

/*<       END IF >*/
    }

/*<    50 CONTINUE >*/
L50:

/*<       WORK( 1 ) = MAXWRK >*/
    work[1] = (doublereal) maxwrk;

/*<       RETURN >*/
    return 0;

/*     End of DGGES */

/*<       END >*/
} /* dgges_ */

#ifdef __cplusplus
        }
#endif