lbfgs.c

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

/*                        STPMIN,STPMAX,INFOC) */

/*     WHERE */

/*       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP, */
/*         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED */
/*         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION */
/*         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE */
/*         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. */

/*       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP, */
/*         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF */
/*         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE */
/*         UPDATED APPROPRIATELY. */

/*       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP, */
/*         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP. */
/*         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE */
/*         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP. */

/*       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER */
/*         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED */
/*         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER */
/*         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE. */

/*       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER */
/*         AND UPPER BOUNDS FOR THE STEP. */

/*       INFOC IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: */
/*         IF INFOC = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED */
/*         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE */
/*         INFOC = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS. */

/*     SUBPROGRAMS CALLED */

/*       FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT */

/*     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 */
/*     JORGE J. MORE', DAVID J. THUENTE */

/*<       DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA >*/
/*<       INFOC = 0 >*/
    infoc = 0;

/*     CHECK THE INPUT PARAMETERS FOR ERRORS. */

/*<    >*/
    if ((brackt && (stp <= min(stx,sty) || stp >= max(stx,sty))) || *dx *
             (stp - stx) >= (float)0. || *stpmax < *stpmin) {
        return 0;
    }

/*     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. */

/*<       SGND = DP*(DX/ABS(DX)) >*/
    sgnd = *dp * (*dx / abs(*dx));

/*     FIRST CASE. A HIGHER FUNCTION VALUE. */
/*     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER */
/*     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, */
/*     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. */

/*<       IF (FP .GT. FX) THEN >*/
    if (*fp > fx) {
/*<          INFOC = 1 >*/
        infoc = 1;
/*<          MCBOUND = .TRUE. >*/
        mcbound = TRUE_;
/*<          THETA = 3*(FX - FP)/(STP - STX) + DX + DP >*/
        theta = (fx - *fp) * 3 / (stp - stx) + *dx + *dp;
/*<          S = MAX(ABS(THETA),ABS(DX),ABS(DP)) >*/
/* Computing MAX */
        d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                *dp);
        s = max(d__1,d__2);
/*<          GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) >*/
/* Computing 2nd power */
        d__1 = theta / s;
        gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
/*<          IF (STP .LT. STX) GAMMA = -GAMMA >*/
        if (stp < stx) {
            gamma = -gamma;
        }
/*<          P = (GAMMA - DX) + THETA >*/
        p = gamma - *dx + theta;
/*<          Q = ((GAMMA - DX) + GAMMA) + DP >*/
        q = gamma - *dx + gamma + *dp;
/*<          R = P/Q >*/
        r__ = p / q;
/*<          STPC = STX + R*(STP - STX) >*/
        stpc = stx + r__ * (stp - stx);
/*<          STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX) >*/
        stpq = stx + *dx / ((fx - *fp) / (stp - stx) + *dx) / 2 * (stp - 
                stx);
/*<          IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN >*/
        if ((d__1 = stpc - stx, abs(d__1)) < (d__2 = stpq - stx, abs(d__2)))
                 {
/*<             STPF = STPC >*/
            stpf = stpc;
/*<          ELSE >*/
        } else {
/*<            STPF = STPC + (STPQ - STPC)/2 >*/
            stpf = stpc + (stpq - stpc) / 2;
/*<            END IF >*/
        }
/*<          BRACKT = .TRUE. >*/
        brackt = TRUE_;

/*     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF */
/*     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC */
/*     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, */
/*     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. */

/*<       ELSE IF (SGND .LT. 0.0) THEN >*/
    } else if (sgnd < (float)0.) {
/*<          INFOC = 2 >*/
        infoc = 2;
/*<          MCBOUND = .FALSE. >*/
        mcbound = FALSE_;
/*<          THETA = 3*(FX - FP)/(STP - STX) + DX + DP >*/
        theta = (fx - *fp) * 3 / (stp - stx) + *dx + *dp;
/*<          S = MAX(ABS(THETA),ABS(DX),ABS(DP)) >*/
/* Computing MAX */
        d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                *dp);
        s = max(d__1,d__2);
/*<          GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) >*/
/* Computing 2nd power */
        d__1 = theta / s;
        gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
/*<          IF (STP .GT. STX) GAMMA = -GAMMA >*/
        if (stp > stx) {
            gamma = -gamma;
        }
/*<          P = (GAMMA - DP) + THETA >*/
        p = gamma - *dp + theta;
/*<          Q = ((GAMMA - DP) + GAMMA) + DX >*/
        q = gamma - *dp + gamma + *dx;
/*<          R = P/Q >*/
        r__ = p / q;
/*<          STPC = STP + R*(STX - STP) >*/
        stpc = stp + r__ * (stx - stp);
/*<          STPQ = STP + (DP/(DP-DX))*(STX - STP) >*/
        stpq = stp + *dp / (*dp - *dx) * (stx - stp);
/*<          IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN >*/
        if ((d__1 = stpc - stp, abs(d__1)) > (d__2 = stpq - stp, abs(d__2)))
                 {
/*<             STPF = STPC >*/
            stpf = stpc;
/*<          ELSE >*/
        } else {
/*<             STPF = STPQ >*/
            stpf = stpq;
/*<             END IF >*/
        }
/*<          BRACKT = .TRUE. >*/
        brackt = TRUE_;

/*     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
/*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. */
/*     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY */
/*     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC */
/*     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE */
/*     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO */
/*     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP */
/*     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. */

/*<       ELSE IF (ABS(DP) .LT. ABS(DX)) THEN >*/
    } else if (abs(*dp) < abs(*dx)) {
/*<          INFOC = 3 >*/
        infoc = 3;
/*<          MCBOUND = .TRUE. >*/
        mcbound = TRUE_;
/*<          THETA = 3*(FX - FP)/(STP - STX) + DX + DP >*/
        theta = (fx - *fp) * 3 / (stp - stx) + *dx + *dp;
/*<          S = MAX(ABS(THETA),ABS(DX),ABS(DP)) >*/
/* Computing MAX */
        d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                *dp);
        s = max(d__1,d__2);

/*        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND */
/*        TO INFINITY IN THE DIRECTION OF THE STEP. */

/*<          GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S))) >*/
/* Computing MAX */
/* Computing 2nd power */
        d__3 = theta / s;
        d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s);
        gamma = s * sqrt((max(d__1,d__2)));
/*<          IF (STP .GT. STX) GAMMA = -GAMMA >*/
        if (stp > stx) {
            gamma = -gamma;
        }
/*<          P = (GAMMA - DP) + THETA >*/
        p = gamma - *dp + theta;
/*<          Q = (GAMMA + (DX - DP)) + GAMMA >*/
        q = gamma + (*dx - *dp) + gamma;
/*<          R = P/Q >*/
        r__ = p / q;
/*<          IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN >*/
        if (r__ < (float)0. && gamma != (float)0.) {
/*<             STPC = STP + R*(STX - STP) >*/
            stpc = stp + r__ * (stx - stp);
/*<          ELSE IF (STP .GT. STX) THEN >*/
        } else if (stp > stx) {
/*<             STPC = STPMAX >*/
            stpc = *stpmax;
/*<          ELSE >*/
        } else {
/*<             STPC = STPMIN >*/
            stpc = *stpmin;
/*<             END IF >*/
        }
/*<          STPQ = STP + (DP/(DP-DX))*(STX - STP) >*/
        stpq = stp + *dp / (*dp - *dx) * (stx - stp);
/*<          IF (BRACKT) THEN >*/
        if (brackt) {
/*<             IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN >*/
            if ((d__1 = stp - stpc, abs(d__1)) < (d__2 = stp - stpq, abs(
                    d__2))) {
/*<                STPF = STPC >*/
                stpf = stpc;
/*<             ELSE >*/
            } else {
/*<                STPF = STPQ >*/
                stpf = stpq;
/*<                END IF >*/
            }
/*<          ELSE >*/
        } else {
/*<             IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN >*/
            if ((d__1 = stp - stpc, abs(d__1)) > (d__2 = stp - stpq, abs(
                    d__2))) {
/*<                STPF = STPC >*/
                stpf = stpc;
/*<             ELSE >*/
            } else {
/*<                STPF = STPQ >*/
                stpf = stpq;
/*<                END IF >*/
            }
/*<             END IF >*/
        }

/*     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
/*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES */
/*     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP */
/*     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. */

/*<       ELSE >*/
    } else {
/*<          INFOC = 4 >*/
        infoc = 4;
/*<          MCBOUND = .FALSE. >*/
        mcbound = FALSE_;
/*<          IF (BRACKT) THEN >*/
        if (brackt) {
/*<             THETA = 3*(FP - FY)/(STY - STP) + DY + DP >*/
            theta = (*fp - fy) * 3 / (sty - stp) + *dy + *dp;
/*<             S = MAX(ABS(THETA),ABS(DY),ABS(DP)) >*/
/* Computing MAX */
            d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 = 
                    abs(*dp);
            s = max(d__1,d__2);
/*<             GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S)) >*/
/* Computing 2nd power */
            d__1 = theta / s;
            gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s));
/*<             IF (STP .GT. STY) GAMMA = -GAMMA >*/
            if (stp > sty) {
                gamma = -gamma;
            }
/*<             P = (GAMMA - DP) + THETA >*/
            p = gamma - *dp + theta;
/*<             Q = ((GAMMA - DP) + GAMMA) + DY >*/
            q = gamma - *dp + gamma + *dy;
/*<             R = P/Q >*/
            r__ = p / q;
/*<             STPC = STP + R*(STY - STP) >*/
            stpc = stp + r__ * (sty - stp);
/*<             STPF = STPC >*/
            stpf = stpc;
/*<          ELSE IF (STP .GT. STX) THEN >*/
        } else if (stp > stx) {
/*<             STPF = STPMAX >*/
            stpf = *stpmax;
/*<          ELSE >*/
        } else {
/*<             STPF = STPMIN >*/
            stpf = *stpmin;
/*<             END IF >*/
        }
/*<          END IF >*/
    }

/*     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT */
/*     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. */

/*<       IF (FP .GT. FX) THEN >*/
    if (*fp > fx) {
/*<          STY = STP >*/
        sty = stp;
/*<          FY = FP >*/
        fy = *fp;
/*<          DY = DP >*/
        *dy = *dp;
/*<       ELSE >*/
    } else {
/*<          IF (SGND .LT. 0.0) THEN >*/
        if (sgnd < (float)0.) {
/*<             STY = STX >*/
            sty = stx;
/*<             FY = FX >*/
            fy = fx;
/*<             DY = DX >*/
            *dy = *dx;
/*<             END IF >*/
        }
/*<          STX = STP >*/
        stx = stp;
/*<          FX = FP >*/
        fx = *fp;
/*<          DX = DP >*/
        *dx = *dp;
/*<          END IF >*/
    }

/*     COMPUTE THE NEW STEP AND SAFEGUARD IT. */

/*<       STPF = MIN(STPMAX,STPF) >*/
    stpf = min(*stpmax,stpf);
/*<       STPF = MAX(STPMIN,STPF) >*/
    stpf = max(*stpmin,stpf);
/*<       STP = STPF >*/
    stp = stpf;
/*<       IF (BRACKT .AND. MCBOUND) THEN >*/
    if (brackt && mcbound) {
/*<          IF (STY .GT. STX) THEN >*/
        if (sty > stx) {
/*<             STP = MIN(STX+0.66*(STY-STX),STP) >*/
/* Computing MIN */
            d__1 = stx + (sty - stx) * (float).66;
            stp = min(d__1,stp);
/*<          ELSE >*/
        } else {
/*<             STP = MAX(STX+0.66*(STY-STX),STP) >*/
/* Computing MAX */
            d__1 = stx + (sty - stx) * (float).66;
            stp = max(d__1,stp);
/*<             END IF >*/
        }
/*<          END IF >*/
    }
/*<       RETURN >*/
    return 0;

/*     LAST LINE OF SUBROUTINE MCSTEP. */

/*<       END >*/
} /* mcstep_ */

#ifdef __cplusplus
        }
#endif