dlag2.c

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

/*<          AS11 = A11 - S2*B11 >*/
        as11 = a11 - s2 * b11;
/*<          SS = A21*( BINV11*BINV22 ) >*/
        ss = a21 * (binv11 * binv22);
/*<          ABI22 = -SS*B12 >*/
        abi22 = -ss * b12;
/*<          PP = HALF*( AS11*BINV11+ABI22 ) >*/
        pp = (as11 * binv11 + abi22) * .5;
/*<          SHIFT = S2 >*/
        shift = s2;
/*<       END IF >*/
    }
/*<       QQ = SS*AS12 >*/
    qq = ss * as12;
/*<       IF( ABS( PP*RTMIN ).GE.ONE ) THEN >*/
    if ((d__1 = pp * rtmin, abs(d__1)) >= 1.) {
/*<          DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN >*/
/* Computing 2nd power */
        d__1 = rtmin * pp;
        discr = d__1 * d__1 + qq * *safmin;
/*<          R = SQRT( ABS( DISCR ) )*RTMAX >*/
        r__ = sqrt((abs(discr))) * rtmax;
/*<       ELSE >*/
    } else {
/*<          IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN >*/
/* Computing 2nd power */
        d__1 = pp;
        if (d__1 * d__1 + abs(qq) <= *safmin) {
/*<             DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX >*/
/* Computing 2nd power */
            d__1 = rtmax * pp;
            discr = d__1 * d__1 + qq * safmax;
/*<             R = SQRT( ABS( DISCR ) )*RTMIN >*/
            r__ = sqrt((abs(discr))) * rtmin;
/*<          ELSE >*/
        } else {
/*<             DISCR = PP**2 + QQ >*/
/* Computing 2nd power */
            d__1 = pp;
            discr = d__1 * d__1 + qq;
/*<             R = SQRT( ABS( DISCR ) ) >*/
            r__ = sqrt((abs(discr)));
/*<          END IF >*/
        }
/*<       END IF >*/
    }

/*     Note: the test of R in the following IF is to cover the case when */
/*           DISCR is small and negative and is flushed to zero during */
/*           the calculation of R.  On machines which have a consistent */
/*           flush-to-zero threshhold and handle numbers above that */
/*           threshhold correctly, it would not be necessary. */

/*<       IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN >*/
    if (discr >= 0. || r__ == 0.) {
/*<          SUM = PP + SIGN( R, PP ) >*/
        sum = pp + d_sign(&r__, &pp);
/*<          DIFF = PP - SIGN( R, PP ) >*/
        diff = pp - d_sign(&r__, &pp);
/*<          WBIG = SHIFT + SUM >*/
        wbig = shift + sum;

/*        Compute smaller eigenvalue */

/*<          WSMALL = SHIFT + DIFF >*/
        wsmall = shift + diff;
/*<          IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN >*/
/* Computing MAX */
        d__1 = abs(wsmall);
        if (abs(wbig) * .5 > max(d__1,*safmin)) {
/*<             WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 ) >*/
            wdet = (a11 * a22 - a12 * a21) * (binv11 * binv22);
/*<             WSMALL = WDET / WBIG >*/
            wsmall = wdet / wbig;
/*<          END IF >*/
        }

/*        Choose (real) eigenvalue closest to 2,2 element of A*B**(-1) */
/*        for WR1. */

/*<          IF( PP.GT.ABI22 ) THEN >*/
        if (pp > abi22) {
/*<             WR1 = MIN( WBIG, WSMALL ) >*/
            *wr1 = min(wbig,wsmall);
/*<             WR2 = MAX( WBIG, WSMALL ) >*/
            *wr2 = max(wbig,wsmall);
/*<          ELSE >*/
        } else {
/*<             WR1 = MAX( WBIG, WSMALL ) >*/
            *wr1 = max(wbig,wsmall);
/*<             WR2 = MIN( WBIG, WSMALL ) >*/
            *wr2 = min(wbig,wsmall);
/*<          END IF >*/
        }
/*<          WI = ZERO >*/
        *wi = 0.;
/*<       ELSE >*/
    } else {

/*        Complex eigenvalues */

/*<          WR1 = SHIFT + PP >*/
        *wr1 = shift + pp;
/*<          WR2 = WR1 >*/
        *wr2 = *wr1;
/*<          WI = R >*/
        *wi = r__;
/*<       END IF >*/
    }

/*     Further scaling to avoid underflow and overflow in computing */
/*     SCALE1 and overflow in computing w*B. */

/*     This scale factor (WSCALE) is bounded from above using C1 and C2, */
/*     and from below using C3 and C4. */
/*        C1 implements the condition  s A  must never overflow. */
/*        C2 implements the condition  w B  must never overflow. */
/*        C3, with C2, */
/*           implement the condition that s A - w B must never overflow. */
/*        C4 implements the condition  s    should not underflow. */
/*        C5 implements the condition  max(s,|w|) should be at least 2. */

/*<       C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) ) >*/
    c1 = bsize * (*safmin * max(1.,ascale));
/*<       C2 = SAFMIN*MAX( ONE, BNORM ) >*/
    c2 = *safmin * max(1.,bnorm);
/*<       C3 = BSIZE*SAFMIN >*/
    c3 = bsize * *safmin;
/*<       IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN >*/
    if (ascale <= 1. && bsize <= 1.) {
/*<          C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE ) >*/
/* Computing MIN */
        d__1 = 1., d__2 = ascale / *safmin * bsize;
        c4 = min(d__1,d__2);
/*<       ELSE >*/
    } else {
/*<          C4 = ONE >*/
        c4 = 1.;
/*<       END IF >*/
    }
/*<       IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN >*/
    if (ascale <= 1. || bsize <= 1.) {
/*<          C5 = MIN( ONE, ASCALE*BSIZE ) >*/
/* Computing MIN */
        d__1 = 1., d__2 = ascale * bsize;
        c5 = min(d__1,d__2);
/*<       ELSE >*/
    } else {
/*<          C5 = ONE >*/
        c5 = 1.;
/*<       END IF >*/
    }

/*     Scale first eigenvalue */

/*<       WABS = ABS( WR1 ) + ABS( WI ) >*/
    wabs = abs(*wr1) + abs(*wi);
/*<    >*/
/* Computing MAX */
/* Computing MIN */
    d__3 = c4, d__4 = max(wabs,c5) * .5;
    d__1 = max(*safmin,c1), d__2 = (wabs * c2 + c3) * 1.0000100000000001, 
            d__1 = max(d__1,d__2), d__2 = min(d__3,d__4);
    wsize = max(d__1,d__2);
/*<       IF( WSIZE.NE.ONE ) THEN >*/
    if (wsize != 1.) {
/*<          WSCALE = ONE / WSIZE >*/
        wscale = 1. / wsize;
/*<          IF( WSIZE.GT.ONE ) THEN >*/
        if (wsize > 1.) {
/*<    >*/
            *scale1 = max(ascale,bsize) * wscale * min(ascale,bsize);
/*<          ELSE >*/
        } else {
/*<    >*/
            *scale1 = min(ascale,bsize) * wscale * max(ascale,bsize);
/*<          END IF >*/
        }
/*<          WR1 = WR1*WSCALE >*/
        *wr1 *= wscale;
/*<          IF( WI.NE.ZERO ) THEN >*/
        if (*wi != 0.) {
/*<             WI = WI*WSCALE >*/
            *wi *= wscale;
/*<             WR2 = WR1 >*/
            *wr2 = *wr1;
/*<             SCALE2 = SCALE1 >*/
            *scale2 = *scale1;
/*<          END IF >*/
        }
/*<       ELSE >*/
    } else {
/*<          SCALE1 = ASCALE*BSIZE >*/
        *scale1 = ascale * bsize;
/*<          SCALE2 = SCALE1 >*/
        *scale2 = *scale1;
/*<       END IF >*/
    }

/*     Scale second eigenvalue (if real) */

/*<       IF( WI.EQ.ZERO ) THEN >*/
    if (*wi == 0.) {
/*<    >*/
/* Computing MAX */
/* Computing MIN */
/* Computing MAX */
        d__5 = abs(*wr2);
        d__3 = c4, d__4 = max(d__5,c5) * .5;
        d__1 = max(*safmin,c1), d__2 = (abs(*wr2) * c2 + c3) * 
                1.0000100000000001, d__1 = max(d__1,d__2), d__2 = min(d__3,
                d__4);
        wsize = max(d__1,d__2);
/*<          IF( WSIZE.NE.ONE ) THEN >*/
        if (wsize != 1.) {
/*<             WSCALE = ONE / WSIZE >*/
            wscale = 1. / wsize;
/*<             IF( WSIZE.GT.ONE ) THEN >*/
            if (wsize > 1.) {
/*<    >*/
                *scale2 = max(ascale,bsize) * wscale * min(ascale,bsize);
/*<             ELSE >*/
            } else {
/*<    >*/
                *scale2 = min(ascale,bsize) * wscale * max(ascale,bsize);
/*<             END IF >*/
            }
/*<             WR2 = WR2*WSCALE >*/
            *wr2 *= wscale;
/*<          ELSE >*/
        } else {
/*<             SCALE2 = ASCALE*BSIZE >*/
            *scale2 = ascale * bsize;
/*<          END IF >*/
        }
/*<       END IF >*/
    }

/*     End of DLAG2 */

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

#ifdef __cplusplus
        }
#endif