dggbal.c

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

    i__1 = *n - k + 1;
    dswap_(&i__1, &a[i__ + k * a_dim1], lda, &a[m + k * a_dim1], lda);
/*<       CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB ) >*/
    i__1 = *n - k + 1;
    dswap_(&i__1, &b[i__ + k * b_dim1], ldb, &b[m + k * b_dim1], ldb);

/*     Permute columns M and J */

/*<   170 CONTINUE >*/
L170:
/*<       RSCALE( M ) = J >*/
    rscale[m] = (doublereal) j;
/*<    >*/
    if (j == m) {
        goto L180;
    }
/*<       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 ) >*/
    dswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);
/*<       CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 ) >*/
    dswap_(&l, &b[j * b_dim1 + 1], &c__1, &b[m * b_dim1 + 1], &c__1);

/*<   180 CONTINUE >*/
L180:
/*<       GO TO ( 20, 90 )IFLOW >*/
    switch (iflow) {
        case 1:  goto L20;
        case 2:  goto L90;
    }

/*<   190 CONTINUE >*/
L190:
/*<       ILO = K >*/
    *ilo = k;
/*<       IHI = L >*/
    *ihi = l;

/*<    >*/
    if (*ilo == *ihi) {
        return 0;
    }

/*<    >*/
    if (lsame_(job, "P", (ftnlen)1, (ftnlen)1)) {
        return 0;
    }

/*     Balance the submatrix in rows ILO to IHI. */

/*<       NR = IHI - ILO + 1 >*/
    nr = *ihi - *ilo + 1;
/*<       DO 200 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          RSCALE( I ) = ZERO >*/
        rscale[i__] = 0.;
/*<          LSCALE( I ) = ZERO >*/
        lscale[i__] = 0.;

/*<          WORK( I ) = ZERO >*/
        work[i__] = 0.;
/*<          WORK( I+N ) = ZERO >*/
        work[i__ + *n] = 0.;
/*<          WORK( I+2*N ) = ZERO >*/
        work[i__ + (*n << 1)] = 0.;
/*<          WORK( I+3*N ) = ZERO >*/
        work[i__ + *n * 3] = 0.;
/*<          WORK( I+4*N ) = ZERO >*/
        work[i__ + (*n << 2)] = 0.;
/*<          WORK( I+5*N ) = ZERO >*/
        work[i__ + *n * 5] = 0.;
/*<   200 CONTINUE >*/
/* L200: */
    }

/*     Compute right side vector in resulting linear equations */

/*<       BASL = LOG10( SCLFAC ) >*/
    basl = d_lg10(&c_b34);
/*<       DO 240 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          DO 230 J = ILO, IHI >*/
        i__2 = *ihi;
        for (j = *ilo; j <= i__2; ++j) {
/*<             TB = B( I, J ) >*/
            tb = b[i__ + j * b_dim1];
/*<             TA = A( I, J ) >*/
            ta = a[i__ + j * a_dim1];
/*<    >*/
            if (ta == 0.) {
                goto L210;
            }
/*<             TA = LOG10( ABS( TA ) ) / BASL >*/
            d__1 = abs(ta);
            ta = d_lg10(&d__1) / basl;
/*<   210       CONTINUE >*/
L210:
/*<    >*/
            if (tb == 0.) {
                goto L220;
            }
/*<             TB = LOG10( ABS( TB ) ) / BASL >*/
            d__1 = abs(tb);
            tb = d_lg10(&d__1) / basl;
/*<   220       CONTINUE >*/
L220:
/*<             WORK( I+4*N ) = WORK( I+4*N ) - TA - TB >*/
            work[i__ + (*n << 2)] = work[i__ + (*n << 2)] - ta - tb;
/*<             WORK( J+5*N ) = WORK( J+5*N ) - TA - TB >*/
            work[j + *n * 5] = work[j + *n * 5] - ta - tb;
/*<   230    CONTINUE >*/
/* L230: */
        }
/*<   240 CONTINUE >*/
/* L240: */
    }

/*<       COEF = ONE / DBLE( 2*NR ) >*/
    coef = 1. / (doublereal) (nr << 1);
/*<       COEF2 = COEF*COEF >*/
    coef2 = coef * coef;
/*<       COEF5 = HALF*COEF2 >*/
    coef5 = coef2 * .5;
/*<       NRP2 = NR + 2 >*/
    nrp2 = nr + 2;
/*<       BETA = ZERO >*/
    beta = 0.;
/*<       IT = 1 >*/
    it = 1;

/*     Start generalized conjugate gradient iteration */

/*<   250 CONTINUE >*/
L250:

/*<    >*/
    gamma = ddot_(&nr, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + (*n << 2)]
            , &c__1) + ddot_(&nr, &work[*ilo + *n * 5], &c__1, &work[*ilo + *
            n * 5], &c__1);

/*<       EW = ZERO >*/
    ew = 0.;
/*<       EWC = ZERO >*/
    ewc = 0.;
/*<       DO 260 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          EW = EW + WORK( I+4*N ) >*/
        ew += work[i__ + (*n << 2)];
/*<          EWC = EWC + WORK( I+5*N ) >*/
        ewc += work[i__ + *n * 5];
/*<   260 CONTINUE >*/
/* L260: */
    }

/*<       GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2 >*/
/* Computing 2nd power */
    d__1 = ew;
/* Computing 2nd power */
    d__2 = ewc;
/* Computing 2nd power */
    d__3 = ew - ewc;
    gamma = coef * gamma - coef2 * (d__1 * d__1 + d__2 * d__2) - coef5 * (
            d__3 * d__3);
/*<    >*/
    if (gamma == 0.) {
        goto L350;
    }
/*<    >*/
    if (it != 1) {
        beta = gamma / pgamma;
    }
/*<       T = COEF5*( EWC-THREE*EW ) >*/
    t = coef5 * (ewc - ew * 3.);
/*<       TC = COEF5*( EW-THREE*EWC ) >*/
    tc = coef5 * (ew - ewc * 3.);

/*<       CALL DSCAL( NR, BETA, WORK( ILO ), 1 ) >*/
    dscal_(&nr, &beta, &work[*ilo], &c__1);
/*<       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 ) >*/
    dscal_(&nr, &beta, &work[*ilo + *n], &c__1);

/*<       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 ) >*/
    daxpy_(&nr, &coef, &work[*ilo + (*n << 2)], &c__1, &work[*ilo + *n], &
            c__1);
/*<       CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 ) >*/
    daxpy_(&nr, &coef, &work[*ilo + *n * 5], &c__1, &work[*ilo], &c__1);

/*<       DO 270 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          WORK( I ) = WORK( I ) + TC >*/
        work[i__] += tc;
/*<          WORK( I+N ) = WORK( I+N ) + T >*/
        work[i__ + *n] += t;
/*<   270 CONTINUE >*/
/* L270: */
    }

/*     Apply matrix to vector */

/*<       DO 300 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          KOUNT = 0 >*/
        kount = 0;
/*<          SUM = ZERO >*/
        sum = 0.;
/*<          DO 290 J = ILO, IHI >*/
        i__2 = *ihi;
        for (j = *ilo; j <= i__2; ++j) {
/*<    >*/
            if (a[i__ + j * a_dim1] == 0.) {
                goto L280;
            }
/*<             KOUNT = KOUNT + 1 >*/
            ++kount;
/*<             SUM = SUM + WORK( J ) >*/
            sum += work[j];
/*<   280       CONTINUE >*/
L280:
/*<    >*/
            if (b[i__ + j * b_dim1] == 0.) {
                goto L290;
            }
/*<             KOUNT = KOUNT + 1 >*/
            ++kount;
/*<             SUM = SUM + WORK( J ) >*/
            sum += work[j];
/*<   290    CONTINUE >*/
L290:
            ;
        }
/*<          WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM >*/
        work[i__ + (*n << 1)] = (doublereal) kount * work[i__ + *n] + sum;
/*<   300 CONTINUE >*/
/* L300: */
    }

/*<       DO 330 J = ILO, IHI >*/
    i__1 = *ihi;
    for (j = *ilo; j <= i__1; ++j) {
/*<          KOUNT = 0 >*/
        kount = 0;
/*<          SUM = ZERO >*/
        sum = 0.;
/*<          DO 320 I = ILO, IHI >*/
        i__2 = *ihi;
        for (i__ = *ilo; i__ <= i__2; ++i__) {
/*<    >*/
            if (a[i__ + j * a_dim1] == 0.) {
                goto L310;
            }
/*<             KOUNT = KOUNT + 1 >*/
            ++kount;
/*<             SUM = SUM + WORK( I+N ) >*/
            sum += work[i__ + *n];
/*<   310       CONTINUE >*/
L310:
/*<    >*/
            if (b[i__ + j * b_dim1] == 0.) {
                goto L320;
            }
/*<             KOUNT = KOUNT + 1 >*/
            ++kount;
/*<             SUM = SUM + WORK( I+N ) >*/
            sum += work[i__ + *n];
/*<   320    CONTINUE >*/
L320:
            ;
        }
/*<          WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM >*/
        work[j + *n * 3] = (doublereal) kount * work[j] + sum;
/*<   330 CONTINUE >*/
/* L330: */
    }

/*<    >*/
    sum = ddot_(&nr, &work[*ilo + *n], &c__1, &work[*ilo + (*n << 1)], &c__1) 
            + ddot_(&nr, &work[*ilo], &c__1, &work[*ilo + *n * 3], &c__1);
/*<       ALPHA = GAMMA / SUM >*/
    alpha = gamma / sum;

/*     Determine correction to current iteration */

/*<       CMAX = ZERO >*/
    cmax = 0.;
/*<       DO 340 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          COR = ALPHA*WORK( I+N ) >*/
        cor = alpha * work[i__ + *n];
/*<    >*/
        if (abs(cor) > cmax) {
            cmax = abs(cor);
        }
/*<          LSCALE( I ) = LSCALE( I ) + COR >*/
        lscale[i__] += cor;
/*<          COR = ALPHA*WORK( I ) >*/
        cor = alpha * work[i__];
/*<    >*/
        if (abs(cor) > cmax) {
            cmax = abs(cor);
        }
/*<          RSCALE( I ) = RSCALE( I ) + COR >*/
        rscale[i__] += cor;
/*<   340 CONTINUE >*/
/* L340: */
    }
/*<    >*/
    if (cmax < .5) {
        goto L350;
    }

/*<       CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 ) >*/
    d__1 = -alpha;
    daxpy_(&nr, &d__1, &work[*ilo + (*n << 1)], &c__1, &work[*ilo + (*n << 2)]
            , &c__1);
/*<       CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 ) >*/
    d__1 = -alpha;
    daxpy_(&nr, &d__1, &work[*ilo + *n * 3], &c__1, &work[*ilo + *n * 5], &
            c__1);

/*<       PGAMMA = GAMMA >*/
    pgamma = gamma;
/*<       IT = IT + 1 >*/
    ++it;
/*<    >*/
    if (it <= nrp2) {
        goto L250;
    }

/*     End generalized conjugate gradient iteration */

/*<   350 CONTINUE >*/
L350:
/*<       SFMIN = DLAMCH( 'S' ) >*/
    sfmin = dlamch_("S", (ftnlen)1);
/*<       SFMAX = ONE / SFMIN >*/
    sfmax = 1. / sfmin;
/*<       LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE ) >*/
    lsfmin = (integer) (d_lg10(&sfmin) / basl + 1.);
/*<       LSFMAX = INT( LOG10( SFMAX ) / BASL ) >*/
    lsfmax = (integer) (d_lg10(&sfmax) / basl);
/*<       DO 360 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA ) >*/
        i__2 = *n - *ilo + 1;
        irab = idamax_(&i__2, &a[i__ + *ilo * a_dim1], lda);
/*<          RAB = ABS( A( I, IRAB+ILO-1 ) ) >*/
        rab = (d__1 = a[i__ + (irab + *ilo - 1) * a_dim1], abs(d__1));
/*<          IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDA ) >*/
        i__2 = *n - *ilo + 1;
        irab = idamax_(&i__2, &b[i__ + *ilo * b_dim1], lda);
/*<          RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) ) >*/
/* Computing MAX */
        d__2 = rab, d__3 = (d__1 = b[i__ + (irab + *ilo - 1) * b_dim1], abs(
                d__1));
        rab = max(d__2,d__3);
/*<          LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE ) >*/
        d__1 = rab + sfmin;
        lrab = (integer) (d_lg10(&d__1) / basl + 1.);
/*<          IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) ) >*/
        ir = (integer) (lscale[i__] + d_sign(&c_b70, &lscale[i__]));
/*<          IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB ) >*/
/* Computing MIN */
        i__2 = max(ir,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lrab;
        ir = min(i__2,i__3);
/*<          LSCALE( I ) = SCLFAC**IR >*/
        lscale[i__] = pow_di(&c_b34, &ir);
/*<          ICAB = IDAMAX( IHI, A( 1, I ), 1 ) >*/
        icab = idamax_(ihi, &a[i__ * a_dim1 + 1], &c__1);
/*<          CAB = ABS( A( ICAB, I ) ) >*/
        cab = (d__1 = a[icab + i__ * a_dim1], abs(d__1));
/*<          ICAB = IDAMAX( IHI, B( 1, I ), 1 ) >*/
        icab = idamax_(ihi, &b[i__ * b_dim1 + 1], &c__1);
/*<          CAB = MAX( CAB, ABS( B( ICAB, I ) ) ) >*/
/* Computing MAX */
        d__2 = cab, d__3 = (d__1 = b[icab + i__ * b_dim1], abs(d__1));
        cab = max(d__2,d__3);
/*<          LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE ) >*/
        d__1 = cab + sfmin;
        lcab = (integer) (d_lg10(&d__1) / basl + 1.);
/*<          JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) ) >*/
        jc = (integer) (rscale[i__] + d_sign(&c_b70, &rscale[i__]));
/*<          JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB ) >*/
/* Computing MIN */
        i__2 = max(jc,lsfmin), i__2 = min(i__2,lsfmax), i__3 = lsfmax - lcab;
        jc = min(i__2,i__3);
/*<          RSCALE( I ) = SCLFAC**JC >*/
        rscale[i__] = pow_di(&c_b34, &jc);
/*<   360 CONTINUE >*/
/* L360: */
    }

/*     Row scaling of matrices A and B */

/*<       DO 370 I = ILO, IHI >*/
    i__1 = *ihi;
    for (i__ = *ilo; i__ <= i__1; ++i__) {
/*<          CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA ) >*/
        i__2 = *n - *ilo + 1;
        dscal_(&i__2, &lscale[i__], &a[i__ + *ilo * a_dim1], lda);
/*<          CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB ) >*/
        i__2 = *n - *ilo + 1;
        dscal_(&i__2, &lscale[i__], &b[i__ + *ilo * b_dim1], ldb);
/*<   370 CONTINUE >*/
/* L370: */
    }

/*     Column scaling of matrices A and B */

/*<       DO 380 J = ILO, IHI >*/
    i__1 = *ihi;
    for (j = *ilo; j <= i__1; ++j) {
/*<          CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 ) >*/
        dscal_(ihi, &rscale[j], &a[j * a_dim1 + 1], &c__1);
/*<          CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 ) >*/
        dscal_(ihi, &rscale[j], &b[j * b_dim1 + 1], &c__1);
/*<   380 CONTINUE >*/
/* L380: */
    }

/*<       RETURN >*/
    return 0;

/*     End of DGGBAL */

/*<       END >*/
} /* dggbal_ */

#ifdef __cplusplus
        }
#endif