dtgex2.c

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

    eps = dlamch_("P", (ftnlen)1);
/*<       SMLNUM = DLAMCH( 'S' ) / EPS >*/
    smlnum = dlamch_("S", (ftnlen)1) / eps;
/*<       DSCALE = ZERO >*/
    dscale = 0.;
/*<       DSUM = ONE >*/
    dsum = 1.;
/*<       CALL DLACPY( 'Full', M, M, S, LDST, WORK, M ) >*/
    dlacpy_("Full", &m, &m, s, &c__4, &work[1], &m, (ftnlen)4);
/*<       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) >*/
    i__1 = m * m;
    dlassq_(&i__1, &work[1], &c__1, &dscale, &dsum);
/*<       CALL DLACPY( 'Full', M, M, T, LDST, WORK, M ) >*/
    dlacpy_("Full", &m, &m, t, &c__4, &work[1], &m, (ftnlen)4);
/*<       CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) >*/
    i__1 = m * m;
    dlassq_(&i__1, &work[1], &c__1, &dscale, &dsum);
/*<       DNORM = DSCALE*SQRT( DSUM ) >*/
    dnorm = dscale * sqrt(dsum);
/*<       THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) >*/
/* Computing MAX */
    d__1 = eps * 10. * dnorm;
    thresh = max(d__1,smlnum);

/*<       IF( M.EQ.2 ) THEN >*/
    if (m == 2) {

/*        CASE 1: Swap 1-by-1 and 1-by-1 blocks. */

/*        Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks */
/*        using Givens rotations and perform the swap tentatively. */

/*<          F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 ) >*/
        f = s[5] * t[0] - t[5] * s[0];
/*<          G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 ) >*/
        g = s[5] * t[4] - t[5] * s[4];
/*<          SB = ABS( T( 2, 2 ) ) >*/
        sb = abs(t[5]);
/*<          SA = ABS( S( 2, 2 ) ) >*/
        sa = abs(s[5]);
/*<          CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM ) >*/
        dlartg_(&f, &g, &ir[4], ir, &ddum);
/*<          IR( 2, 1 ) = -IR( 1, 2 ) >*/
        ir[1] = -ir[4];
/*<          IR( 2, 2 ) = IR( 1, 1 ) >*/
        ir[5] = ir[0];
/*<    >*/
        drot_(&c__2, s, &c__1, &s[4], &c__1, ir, &ir[1]);
/*<    >*/
        drot_(&c__2, t, &c__1, &t[4], &c__1, ir, &ir[1]);
/*<          IF( SA.GE.SB ) THEN >*/
        if (sa >= sb) {
/*<    >*/
            dlartg_(s, &s[1], li, &li[1], &ddum);
/*<          ELSE >*/
        } else {
/*<    >*/
            dlartg_(t, &t[1], li, &li[1], &ddum);
/*<          END IF >*/
        }
/*<    >*/
        drot_(&c__2, s, &c__4, &s[1], &c__4, li, &li[1]);
/*<    >*/
        drot_(&c__2, t, &c__4, &t[1], &c__4, li, &li[1]);
/*<          LI( 2, 2 ) = LI( 1, 1 ) >*/
        li[5] = li[0];
/*<          LI( 1, 2 ) = -LI( 2, 1 ) >*/
        li[4] = -li[1];

/*        Weak stability test: */
/*           |S21| + |T21| <= O(EPS * F-norm((S, T))) */

/*<          WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) ) >*/
        ws = abs(s[1]) + abs(t[1]);
/*<          WEAK = WS.LE.THRESH >*/
        weak = ws <= thresh;
/*<    >*/
        if (! weak) {
            goto L70;
        }

/*<          IF( WANDS ) THEN >*/
        if (TRUE_) {

/*           Strong stability test: */
/*             F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A,B))) */

/*<    >*/
            dlacpy_("Full", &m, &m, &a[*j1 + *j1 * a_dim1], lda, &work[m * m 
                    + 1], &m, (ftnlen)4);
/*<    >*/
            dgemm_("N", "N", &m, &m, &m, &c_b38, li, &c__4, s, &c__4, &c_b3, &
                    work[1], &m, (ftnlen)1, (ftnlen)1);
/*<    >*/
            dgemm_("N", "T", &m, &m, &m, &c_b44, &work[1], &m, ir, &c__4, &
                    c_b38, &work[m * m + 1], &m, (ftnlen)1, (ftnlen)1);
/*<             DSCALE = ZERO >*/
            dscale = 0.;
/*<             DSUM = ONE >*/
            dsum = 1.;
/*<             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) >*/
            i__1 = m * m;
            dlassq_(&i__1, &work[m * m + 1], &c__1, &dscale, &dsum);

/*<    >*/
            dlacpy_("Full", &m, &m, &b[*j1 + *j1 * b_dim1], ldb, &work[m * m 
                    + 1], &m, (ftnlen)4);
/*<    >*/
            dgemm_("N", "N", &m, &m, &m, &c_b38, li, &c__4, t, &c__4, &c_b3, &
                    work[1], &m, (ftnlen)1, (ftnlen)1);
/*<    >*/
            dgemm_("N", "T", &m, &m, &m, &c_b44, &work[1], &m, ir, &c__4, &
                    c_b38, &work[m * m + 1], &m, (ftnlen)1, (ftnlen)1);
/*<             CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) >*/
            i__1 = m * m;
            dlassq_(&i__1, &work[m * m + 1], &c__1, &dscale, &dsum);
/*<             SS = DSCALE*SQRT( DSUM ) >*/
            ss = dscale * sqrt(dsum);
/*<             DTRONG = SS.LE.THRESH >*/
            dtrong = ss <= thresh;
/*<    >*/
            if (! dtrong) {
                goto L70;
            }
/*<          END IF >*/
        }

/*        Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and */
/*               (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)). */

/*<    >*/
        i__1 = *j1 + 1;
        drot_(&i__1, &a[*j1 * a_dim1 + 1], &c__1, &a[(*j1 + 1) * a_dim1 + 1], 
                &c__1, ir, &ir[1]);
/*<    >*/
        i__1 = *j1 + 1;
        drot_(&i__1, &b[*j1 * b_dim1 + 1], &c__1, &b[(*j1 + 1) * b_dim1 + 1], 
                &c__1, ir, &ir[1]);
/*<    >*/
        i__1 = *n - *j1 + 1;
        drot_(&i__1, &a[*j1 + *j1 * a_dim1], lda, &a[*j1 + 1 + *j1 * a_dim1], 
                lda, li, &li[1]);
/*<    >*/
        i__1 = *n - *j1 + 1;
        drot_(&i__1, &b[*j1 + *j1 * b_dim1], ldb, &b[*j1 + 1 + *j1 * b_dim1], 
                ldb, li, &li[1]);

/*        Set  N1-by-N2 (2,1) - blocks to ZERO. */

/*<          A( J1+1, J1 ) = ZERO >*/
        a[*j1 + 1 + *j1 * a_dim1] = 0.;
/*<          B( J1+1, J1 ) = ZERO >*/
        b[*j1 + 1 + *j1 * b_dim1] = 0.;

/*        Accumulate transformations into Q and Z if requested. */

/*<    >*/
        if (*wantz) {
            drot_(n, &z__[*j1 * z_dim1 + 1], &c__1, &z__[(*j1 + 1) * z_dim1 + 
                    1], &c__1, ir, &ir[1]);
        }
/*<    >*/
        if (*wantq) {
            drot_(n, &q[*j1 * q_dim1 + 1], &c__1, &q[(*j1 + 1) * q_dim1 + 1], 
                    &c__1, li, &li[1]);
        }

/*        Exit with INFO = 0 if swap was successfully performed. */

/*<          RETURN >*/
        return 0;

/*<       ELSE >*/
    } else {

/*        CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2 */
/*                and 2-by-2 blocks. */

/*        Solve the generalized Sylvester equation */
/*                 S11 * R - L * S22 = SCALE * S12 */
/*                 T11 * R - L * T22 = SCALE * T12 */
/*        for R and L. Solutions in LI and IR. */

/*<          CALL DLACPY( 'Full', N1, N2, T( 1, N1+1 ), LDST, LI, LDST ) >*/
        dlacpy_("Full", n1, n2, &t[((*n1 + 1) << 2) - 4], &c__4, li, &c__4, (
                ftnlen)4);
/*<    >*/
        dlacpy_("Full", n1, n2, &s[((*n1 + 1) << 2) - 4], &c__4, &ir[*n2 + 1 + (
                (*n1 + 1) << 2) - 5], &c__4, (ftnlen)4);
/*<    >*/
        dtgsy2_("N", &c__0, n1, n2, s, &c__4, &s[(*n1 + 1) + ((*n1 + 1) << 2) - 5]
                , &c__4, &ir[*n2 + 1 + ((*n1 + 1) << 2) - 5], &c__4, t, &c__4, &
                t[(*n1 + 1) + ((*n1 + 1) << 2) - 5], &c__4, li, &c__4, &scale, &
                dsum, &dscale, iwork, &idum, &linfo, (ftnlen)1);

/*        Compute orthogonal matrix QL: */

/*                    QL' * LI = [ TL ] */
/*                               [ 0  ] */
/*        where */
/*                    LI =  [      -L              ] */
/*                          [ SCALE * identity(N2) ] */

/*<          DO 10 I = 1, N2 >*/
        i__1 = *n2;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             CALL DSCAL( N1, -ONE, LI( 1, I ), 1 ) >*/
            dscal_(n1, &c_b44, &li[(i__ << 2) - 4], &c__1);
/*<             LI( N1+I, I ) = SCALE >*/
            li[*n1 + i__ + (i__ << 2) - 5] = scale;
/*<    10    CONTINUE >*/
/* L10: */
        }
/*<          CALL DGEQR2( M, N2, LI, LDST, TAUL, WORK, LINFO ) >*/
        dgeqr2_(&m, n2, li, &c__4, taul, &work[1], &linfo);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }
/*<          CALL DORG2R( M, M, N2, LI, LDST, TAUL, WORK, LINFO ) >*/
        dorg2r_(&m, &m, n2, li, &c__4, taul, &work[1], &linfo);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }

/*        Compute orthogonal matrix RQ: */

/*                    IR * RQ' =   [ 0  TR], */

/*         where IR = [ SCALE * identity(N1), R ] */

/*<          DO 20 I = 1, N1 >*/
        i__1 = *n1;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             IR( N2+I, I ) = SCALE >*/
            ir[*n2 + i__ + (i__ << 2) - 5] = scale;
/*<    20    CONTINUE >*/
/* L20: */
        }
/*<          CALL DGERQ2( N1, M, IR( N2+1, 1 ), LDST, TAUR, WORK, LINFO ) >*/
        dgerq2_(n1, &m, &ir[*n2], &c__4, taur, &work[1], &linfo);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }
/*<          CALL DORGR2( M, M, N1, IR, LDST, TAUR, WORK, LINFO ) >*/
        dorgr2_(&m, &m, n1, ir, &c__4, taur, &work[1], &linfo);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }

/*        Perform the swapping tentatively: */

/*<    >*/
        dgemm_("T", "N", &m, &m, &m, &c_b38, li, &c__4, s, &c__4, &c_b3, &
                work[1], &m, (ftnlen)1, (ftnlen)1);
/*<    >*/
        dgemm_("N", "T", &m, &m, &m, &c_b38, &work[1], &m, ir, &c__4, &c_b3, 
                s, &c__4, (ftnlen)1, (ftnlen)1);
/*<    >*/
        dgemm_("T", "N", &m, &m, &m, &c_b38, li, &c__4, t, &c__4, &c_b3, &
                work[1], &m, (ftnlen)1, (ftnlen)1);
/*<    >*/
        dgemm_("N", "T", &m, &m, &m, &c_b38, &work[1], &m, ir, &c__4, &c_b3, 
                t, &c__4, (ftnlen)1, (ftnlen)1);
/*<          CALL DLACPY( 'F', M, M, S, LDST, SCPY, LDST ) >*/
        dlacpy_("F", &m, &m, s, &c__4, scpy, &c__4, (ftnlen)1);
/*<          CALL DLACPY( 'F', M, M, T, LDST, TCPY, LDST ) >*/
        dlacpy_("F", &m, &m, t, &c__4, tcpy, &c__4, (ftnlen)1);
/*<          CALL DLACPY( 'F', M, M, IR, LDST, IRCOP, LDST ) >*/
        dlacpy_("F", &m, &m, ir, &c__4, ircop, &c__4, (ftnlen)1);
/*<          CALL DLACPY( 'F', M, M, LI, LDST, LICOP, LDST ) >*/
        dlacpy_("F", &m, &m, li, &c__4, licop, &c__4, (ftnlen)1);

/*        Triangularize the B-part by an RQ factorization. */
/*        Apply transformation (from left) to A-part, giving S. */

/*<          CALL DGERQ2( M, M, T, LDST, TAUR, WORK, LINFO ) >*/
        dgerq2_(&m, &m, t, &c__4, taur, &work[1], &linfo);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }
/*<    >*/
        dormr2_("R", "T", &m, &m, &m, t, &c__4, taur, s, &c__4, &work[1], &
                linfo, (ftnlen)1, (ftnlen)1);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }
/*<    >*/
        dormr2_("L", "N", &m, &m, &m, t, &c__4, taur, ir, &c__4, &work[1], &
                linfo, (ftnlen)1, (ftnlen)1);
/*<    >*/
        if (linfo != 0) {
            goto L70;
        }

/*        Compute F-norm(S21) in BRQA21. (T21 is 0.) */

/*<          DSCALE = ZERO >*/
        dscale = 0.;
/*<          DSUM = ONE >*/
        dsum = 1.;
/*<          DO 30 I = 1, N2 >*/
        i__1 = *n2;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<             CALL DLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM ) >*/
            dlassq_(n1, &s[*n2 + 1 + (i__ << 2) - 5], &c__1, &dscale, &dsum);
/*<    30    CONTINUE >*/
/* L30: */
        }
/*<          BRQA21 = DSCALE*SQRT( DSUM ) >*/
        brqa21 = dscale * sqrt(dsum);

/*        Triangularize the B-part by a QR factorization. */
/*        Apply transformation (from right) to A-part, giving S. */

/*<          CALL DGEQR2( M, M, TCPY, LDST, TAUL, WORK, LINFO ) >*/
        dgeqr2_(&m, &m, tcpy, &c__4, taul, &work[1], &linfo);