dtgsyl.c

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

    *info = 0;
/*<       NOTRAN = LSAME( TRANS, 'N' ) >*/
    notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
/*<       LQUERY = ( LWORK.EQ.-1 ) >*/
    lquery = *lwork == -1;

/*<       IF( ( IJOB.EQ.1 .OR. IJOB.EQ.2 ) .AND. NOTRAN ) THEN >*/
    if ((*ijob == 1 || *ijob == 2) && notran) {
/*<          LWMIN = MAX( 1, 2*M*N ) >*/
/* Computing MAX */
        i__1 = 1, i__2 = (*m << 1) * *n;
        lwmin = max(i__1,i__2);
/*<       ELSE >*/
    } else {
/*<          LWMIN = 1 >*/
        lwmin = 1;
/*<       END IF >*/
    }

/*<       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN >*/
    if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1)) {
/*<          INFO = -1 >*/
        *info = -1;
/*<       ELSE IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN >*/
    } else if (*ijob < 0 || *ijob > 4) {
/*<          INFO = -2 >*/
        *info = -2;
/*<       ELSE IF( M.LE.0 ) THEN >*/
    } else if (*m <= 0) {
/*<          INFO = -3 >*/
        *info = -3;
/*<       ELSE IF( N.LE.0 ) THEN >*/
    } else if (*n <= 0) {
/*<          INFO = -4 >*/
        *info = -4;
/*<       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN >*/
    } else if (*lda < max(1,*m)) {
/*<          INFO = -6 >*/
        *info = -6;
/*<       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN >*/
    } else if (*ldb < max(1,*n)) {
/*<          INFO = -8 >*/
        *info = -8;
/*<       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN >*/
    } else if (*ldc < max(1,*m)) {
/*<          INFO = -10 >*/
        *info = -10;
/*<       ELSE IF( LDD.LT.MAX( 1, M ) ) THEN >*/
    } else if (*ldd < max(1,*m)) {
/*<          INFO = -12 >*/
        *info = -12;
/*<       ELSE IF( LDE.LT.MAX( 1, N ) ) THEN >*/
    } else if (*lde < max(1,*n)) {
/*<          INFO = -14 >*/
        *info = -14;
/*<       ELSE IF( LDF.LT.MAX( 1, M ) ) THEN >*/
    } else if (*ldf < max(1,*m)) {
/*<          INFO = -16 >*/
        *info = -16;
/*<       ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN >*/
    } else if (*lwork < lwmin && ! lquery) {
/*<          INFO = -20 >*/
        *info = -20;
/*<       END IF >*/
    }

/*<       IF( INFO.EQ.0 ) THEN >*/
    if (*info == 0) {
/*<          WORK( 1 ) = LWMIN >*/
        work[1] = (doublereal) lwmin;
/*<       END IF >*/
    }

/*<       IF( INFO.NE.0 ) THEN >*/
    if (*info != 0) {
/*<          CALL XERBLA( 'DTGSYL', -INFO ) >*/
        i__1 = -(*info);
        xerbla_("DTGSYL", &i__1, (ftnlen)6);
/*<          RETURN >*/
        return 0;
/*<       ELSE IF( LQUERY ) THEN >*/
    } else if (lquery) {
/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Determine optimal block sizes MB and NB */

/*<       MB = ILAENV( 2, 'DTGSYL', TRANS, M, N, -1, -1 ) >*/
    mb = ilaenv_(&c__2, "DTGSYL", trans, m, n, &c_n1, &c_n1, (ftnlen)6, (
            ftnlen)1);
/*<       NB = ILAENV( 5, 'DTGSYL', TRANS, M, N, -1, -1 ) >*/
    nb = ilaenv_(&c__5, "DTGSYL", trans, m, n, &c_n1, &c_n1, (ftnlen)6, (
            ftnlen)1);

/*<       ISOLVE = 1 >*/
    isolve = 1;
/*<       IFUNC = 0 >*/
    ifunc = 0;
/*<       IF( IJOB.GE.3 .AND. NOTRAN ) THEN >*/
    if (*ijob >= 3 && notran) {
/*<          IFUNC = IJOB - 2 >*/
        ifunc = *ijob - 2;
/*<          DO 10 J = 1, N >*/
        i__1 = *n;
        for (j = 1; j <= i__1; ++j) {
/*<             CALL DCOPY( M, ZERO, 0, C( 1, J ), 1 ) >*/
            dcopy_(m, &c_b14, &c__0, &c__[j * c_dim1 + 1], &c__1);
/*<             CALL DCOPY( M, ZERO, 0, F( 1, J ), 1 ) >*/
            dcopy_(m, &c_b14, &c__0, &f[j * f_dim1 + 1], &c__1);
/*<    10    CONTINUE >*/
/* L10: */
        }
/*<       ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN >*/
    } else if (*ijob >= 1 && notran) {
/*<          ISOLVE = 2 >*/
        isolve = 2;
/*<       END IF >*/
    }

/*<    >*/
    if ((mb <= 1 && nb <= 1) || (mb >= *m && nb >= *n)) {

/*<          DO 30 IROUND = 1, ISOLVE >*/
        i__1 = isolve;
        for (iround = 1; iround <= i__1; ++iround) {

/*           Use unblocked Level 2 solver */

/*<             DSCALE = ZERO >*/
            dscale = 0.;
/*<             DSUM = ONE >*/
            dsum = 1.;
/*<             PQ = 0 >*/
            pq = 0;
/*<    >*/
            dtgsy2_(trans, &ifunc, m, n, &a[a_offset], lda, &b[b_offset], ldb,
                     &c__[c_offset], ldc, &d__[d_offset], ldd, &e[e_offset], 
                    lde, &f[f_offset], ldf, scale, &dsum, &dscale, &iwork[1], 
                    &pq, info, (ftnlen)1);
/*<             IF( DSCALE.NE.ZERO ) THEN >*/
            if (dscale != 0.) {
/*<                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN >*/
                if (*ijob == 1 || *ijob == 3) {
/*<                   DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) >*/
                    *dif = sqrt((doublereal) ((*m << 1) * *n)) / (dscale * 
                            sqrt(dsum));
/*<                ELSE >*/
                } else {
/*<                   DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) >*/
                    *dif = sqrt((doublereal) pq) / (dscale * sqrt(dsum));
/*<                END IF >*/
                }
/*<             END IF >*/
            }

/*<             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN >*/
            if (isolve == 2 && iround == 1) {
/*<                IFUNC = IJOB >*/
                ifunc = *ijob;
/*<                SCALE2 = SCALE >*/
                scale2 = *scale;
/*<                CALL DLACPY( 'F', M, N, C, LDC, WORK, M ) >*/
                dlacpy_("F", m, n, &c__[c_offset], ldc, &work[1], m, (ftnlen)
                        1);
/*<                CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) >*/
                dlacpy_("F", m, n, &f[f_offset], ldf, &work[*m * *n + 1], m, (
                        ftnlen)1);
/*<                DO 20 J = 1, N >*/
                i__2 = *n;
                for (j = 1; j <= i__2; ++j) {
/*<                   CALL DCOPY( M, ZERO, 0, C( 1, J ), 1 ) >*/
                    dcopy_(m, &c_b14, &c__0, &c__[j * c_dim1 + 1], &c__1);
/*<                   CALL DCOPY( M, ZERO, 0, F( 1, J ), 1 ) >*/
                    dcopy_(m, &c_b14, &c__0, &f[j * f_dim1 + 1], &c__1);
/*<    20          CONTINUE >*/
/* L20: */
                }
/*<             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN >*/
            } else if (isolve == 2 && iround == 2) {
/*<                CALL DLACPY( 'F', M, N, WORK, M, C, LDC ) >*/
                dlacpy_("F", m, n, &work[1], m, &c__[c_offset], ldc, (ftnlen)
                        1);
/*<                CALL DLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) >*/
                dlacpy_("F", m, n, &work[*m * *n + 1], m, &f[f_offset], ldf, (
                        ftnlen)1);
/*<                SCALE = SCALE2 >*/
                *scale = scale2;
/*<             END IF >*/
            }
/*<    30    CONTINUE >*/
/* L30: */
        }

/*<          RETURN >*/
        return 0;
/*<       END IF >*/
    }

/*     Determine block structure of A */

/*<       P = 0 >*/
    p = 0;
/*<       I = 1 >*/
    i__ = 1;
/*<    40 CONTINUE >*/
L40:
/*<    >*/
    if (i__ > *m) {
        goto L50;
    }
/*<       P = P + 1 >*/
    ++p;
/*<       IWORK( P ) = I >*/
    iwork[p] = i__;
/*<       I = I + MB >*/
    i__ += mb;
/*<    >*/
    if (i__ >= *m) {
        goto L50;
    }
/*<    >*/
    if (a[i__ + (i__ - 1) * a_dim1] != 0.) {
        ++i__;
    }
/*<       GO TO 40 >*/
    goto L40;
/*<    50 CONTINUE >*/
L50:

/*<       IWORK( P+1 ) = M + 1 >*/
    iwork[p + 1] = *m + 1;
/*<    >*/
    if (iwork[p] == iwork[p + 1]) {
        --p;
    }

/*     Determine block structure of B */

/*<       Q = P + 1 >*/
    q = p + 1;
/*<       J = 1 >*/
    j = 1;
/*<    60 CONTINUE >*/
L60:
/*<    >*/
    if (j > *n) {
        goto L70;
    }
/*<       Q = Q + 1 >*/
    ++q;
/*<       IWORK( Q ) = J >*/
    iwork[q] = j;
/*<       J = J + NB >*/
    j += nb;
/*<    >*/
    if (j >= *n) {
        goto L70;
    }
/*<    >*/
    if (b[j + (j - 1) * b_dim1] != 0.) {
        ++j;
    }
/*<       GO TO 60 >*/
    goto L60;
/*<    70 CONTINUE >*/
L70:

/*<       IWORK( Q+1 ) = N + 1 >*/
    iwork[q + 1] = *n + 1;
/*<    >*/
    if (iwork[q] == iwork[q + 1]) {
        --q;
    }

/*<       IF( NOTRAN ) THEN >*/
    if (notran) {

/*<          DO 150 IROUND = 1, ISOLVE >*/
        i__1 = isolve;
        for (iround = 1; iround <= i__1; ++iround) {

/*           Solve (I, J)-subsystem */
/*               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) */
/*               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) */
/*           for I = P, P - 1,..., 1; J = 1, 2,..., Q */

/*<             DSCALE = ZERO >*/
            dscale = 0.;
/*<             DSUM = ONE >*/
            dsum = 1.;
/*<             PQ = 0 >*/
            pq = 0;
/*<             SCALE = ONE >*/
            *scale = 1.;
/*<             DO 130 J = P + 2, Q >*/
            i__2 = q;
            for (j = p + 2; j <= i__2; ++j) {
/*<                JS = IWORK( J ) >*/
                js = iwork[j];
/*<                JE = IWORK( J+1 ) - 1 >*/
                je = iwork[j + 1] - 1;
/*<                NB = JE - JS + 1 >*/
                nb = je - js + 1;
/*<                DO 120 I = P, 1, -1 >*/
                for (i__ = p; i__ >= 1; --i__) {
/*<                   IS = IWORK( I ) >*/
                    is = iwork[i__];
/*<                   IE = IWORK( I+1 ) - 1 >*/
                    ie = iwork[i__ + 1] - 1;
/*<                   MB = IE - IS + 1 >*/
                    mb = ie - is + 1;