lbfgs.c

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

/*         MAXFEV BY THE END OF AN ITERATION. */

/*       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: */

/*         INFO = 0  IMPROPER INPUT PARAMETERS. */

/*         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT. */

/*         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE */
/*                   DIRECTIONAL DERIVATIVE CONDITION HOLD. */

/*         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY */
/*                   IS AT MOST XTOL. */

/*         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV. */

/*         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN. */

/*         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX. */

/*         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS. */
/*                   THERE MAY NOT BE A STEP WHICH SATISFIES THE */
/*                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS. */
/*                   TOLERANCES MAY BE TOO SMALL. */

/*       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF */
/*         CALLS TO FCN. */

/*       WA IS A WORK ARRAY OF LENGTH N. */

/*     SUBPROGRAMS CALLED */

/*       MCSTEP */

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

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

/*     ********** */
/*<       INTEGER INFOC,J >*/
/*<       LOGICAL BRACKT,STAGE1 >*/
/*<    >*/
/*<       DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/ >*/
    /* Parameter adjustments */
    --wa;
    --s;
    --g;
    --x;

    /* Function Body */
/*<       IF(INFO.EQ.-1) GO TO 45 >*/
    if (info == -1) {
        goto L45;
    }
/*<       INFOC = 1 >*/
    infoc = 1;

/*     CHECK THE INPUT PARAMETERS FOR ERRORS. */

/*<    >*/ if (*n <= 0 || stp <= zero || ftol < zero || lb3_1.gtol < zero || *xtol 
            < zero || lb3_1.stpmin < zero || lb3_1.stpmax < lb3_1.stpmin ||
        maxfev <= 0) {
        return 0;
    }

/*     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION */
/*     AND CHECK THAT S IS A DESCENT DIRECTION. */

/*<       DGINIT = ZERO >*/
    dginit = zero;
/*<       DO 10 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<          DGINIT = DGINIT + G(J)*S(J) >*/
        dginit += g[j] * s[j];
/*<    10    CONTINUE >*/
/* L10: */
    }
/*<       IF (DGINIT .GE. ZERO) then >*/
    if (dginit >= zero) {
/*<          write(LP,15) >*/
        printf("  THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION\n");
/*<    15    FORMAT(/'  THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION') >*/
/*<          RETURN >*/
        return 0;
/*<          ENDIF >*/
    }

/*     INITIALIZE LOCAL VARIABLES. */

/*<       BRACKT = .FALSE. >*/
    brackt = FALSE_;
/*<       STAGE1 = .TRUE. >*/
    stage1 = TRUE_;
/*<       NFEV = 0 >*/
    nfev = 0;
/*<       FINIT = F >*/
    finit = *f;
/*<       DGTEST = FTOL*DGINIT >*/
    dgtest = ftol * dginit;
/*<       WIDTH = STPMAX - STPMIN >*/
    width = lb3_1.stpmax - lb3_1.stpmin;
/*<       WIDTH1 = WIDTH/P5 >*/
    width1 = width / p5;
/*<       DO 20 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<          WA(J) = X(J) >*/
        wa[j] = x[j];
/*<    20    CONTINUE >*/
/* L20: */
    }

/*     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, */
/*     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. */
/*     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, */
/*     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF */
/*     THE INTERVAL OF UNCERTAINTY. */
/*     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, */
/*     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. */

/*<       STX = ZERO >*/
    stx = zero;
/*<       FX = FINIT >*/
    fx = finit;
/*<       DGX = DGINIT >*/
    dgx = dginit;
/*<       STY = ZERO >*/
    sty = zero;
/*<       FY = FINIT >*/
    fy = finit;
/*<       DGY = DGINIT >*/
    dgy = dginit;

/*     START OF ITERATION. */

/*<    30 CONTINUE >*/
L30:

/*        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND */
/*        TO THE PRESENT INTERVAL OF UNCERTAINTY. */

/*<          IF (BRACKT) THEN >*/
    if (brackt) {
/*<             STMIN = MIN(STX,STY) >*/
        stmin = min(stx,sty);
/*<             STMAX = MAX(STX,STY) >*/
        stmax = max(stx,sty);
/*<          ELSE >*/
    } else {
/*<             STMIN = STX >*/
        stmin = stx;
/*<             STMAX = STP + XTRAPF*(STP - STX) >*/
        stmax = stp + xtrapf * (stp - stx);
/*<             END IF >*/
    }

/*        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. */

/*<          STP = MAX(STP,STPMIN) >*/
    stp = max(stp,lb3_1.stpmin);
/*<          STP = MIN(STP,STPMAX) >*/
    stp = min(stp,lb3_1.stpmax);

/*        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET */
/*        STP BE THE LOWEST POINT OBTAINED SO FAR. */

/*<    >*/
    if ((brackt && (stp <= stmin || stp >= stmax)) || nfev >= maxfev - 1 || 
            infoc == 0 || (brackt && stmax - stmin <= *xtol * stmax)) {
        stp = stx;
    }

/*        EVALUATE THE FUNCTION AND GRADIENT AT STP */
/*        AND COMPUTE THE DIRECTIONAL DERIVATIVE. */
/*        We return to main program to obtain F and G. */

/*<          DO 40 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<             X(J) = WA(J) + STP*S(J) >*/
        x[j] = wa[j] + stp * s[j];
/*<    40       CONTINUE >*/
/* L40: */
    }
/*<          INFO=-1 >*/
    info = -1;
/*<          RETURN >*/
    return 0;

/*<    45    INFO=0 >*/
L45:
    info = 0;
/*<          NFEV = NFEV + 1 >*/
    ++(nfev);
/*<          DG = ZERO >*/
    dg = zero;
/*<          DO 50 J = 1, N >*/
    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
/*<             DG = DG + G(J)*S(J) >*/
        dg += g[j] * s[j];
/*<    50       CONTINUE >*/
/* L50: */
    }
/*<          FTEST1 = FINIT + STP*DGTEST >*/
    ftest1 = finit + stp * dgtest;

/*        TEST FOR CONVERGENCE. */

/*<    >*/
    if ((brackt && (stp <= stmin || stp >= stmax)) || infoc == 0) {
        info = 6;
    }
/*<    >*/
    if (stp == lb3_1.stpmax && *f <= ftest1 && dg <= dgtest) {
        info = 5;
    }
/*<    >*/
    if (stp == lb3_1.stpmin && (*f > ftest1 || dg >= dgtest)) {
        info = 4;
    }
/*<          IF (NFEV .GE. MAXFEV) INFO = 3 >*/
    if (nfev >= maxfev) {
        info = 3;
    }
/*<          IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2 >*/
    if (brackt && stmax - stmin <= *xtol * stmax) {
        info = 2;
    }
/*<          IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1 >*/
    if (*f <= ftest1 && abs(dg) <= lb3_1.gtol * (-dginit)) {
        info = 1;
    }

/*        CHECK FOR TERMINATION. */

/*<          IF (INFO .NE. 0) RETURN >*/
    if (info != 0) {
        return 0;
    }

/*        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED */
/*        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. */

/*<    >*/
    if (stage1 && *f <= ftest1 && dg >= min(ftol,lb3_1.gtol) * dginit) {
        stage1 = FALSE_;
    }

/*        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF */
/*        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED */
/*        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE */
/*        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN */
/*        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. */

/*<          IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN >*/
    if (stage1 && *f <= fx && *f > ftest1) {

/*           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. */

/*<             FM = F - STP*DGTEST >*/
        fm = *f - stp * dgtest;
/*<             FXM = FX - STX*DGTEST >*/
        fxm = fx - stx * dgtest;
/*<             FYM = FY - STY*DGTEST >*/
        fym = fy - sty * dgtest;
/*<             DGM = DG - DGTEST >*/
        dgm = dg - dgtest;
/*<             DGXM = DGX - DGTEST >*/
        dgxm = dgx - dgtest;
/*<             DGYM = DGY - DGTEST >*/
        dgym = dgy - dgtest;

/*           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY */
/*           AND TO COMPUTE THE NEW STEP. */

/*<    >*/
        mcstep_(&dgxm, &dgym, &fm, &dgm, &stmin, &stmax,
                v3p_netlib_lbfgs_global_arg);

/*           RESET THE FUNCTION AND GRADIENT VALUES FOR F. */

/*<             FX = FXM + STX*DGTEST >*/
        fx = fxm + stx * dgtest;
/*<             FY = FYM + STY*DGTEST >*/
        fy = fym + sty * dgtest;
/*<             DGX = DGXM + DGTEST >*/
        dgx = dgxm + dgtest;
/*<             DGY = DGYM + DGTEST >*/
        dgy = dgym + dgtest;
/*<          ELSE >*/
    } else {

/*           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY */
/*           AND TO COMPUTE THE NEW STEP. */

/*<    >*/
        mcstep_(&dgx, &dgy, f, &dg, &stmin, &stmax,
                v3p_netlib_lbfgs_global_arg);
/*<             END IF >*/
    }

/*        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE */
/*        INTERVAL OF UNCERTAINTY. */

/*<          IF (BRACKT) THEN >*/
    if (brackt) {
/*<    >*/
        if ((d__1 = sty - stx, abs(d__1)) >= p66 * width1) {
            stp = stx + p5 * (sty - stx);
        }
/*<             WIDTH1 = WIDTH >*/
        width1 = width;
/*<             WIDTH = ABS(STY-STX) >*/
        width = (d__1 = sty - stx, abs(d__1));
/*<             END IF >*/
    }

/*        END OF ITERATION. */

/*<          GO TO 30 >*/
    goto L30;

/*     LAST LINE OF SUBROUTINE MCSRCH. */

/*<       END >*/
} /* mcsrch_ */

/*<    >*/
/* Subroutine */ int mcstep_(
  doublereal *dx, doublereal *dy,
  doublereal *fp, doublereal *dp,
  doublereal *stpmin, doublereal *stpmax,
  v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg
  )
{
    /* System generated locals */
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    doublereal p, q, r__, s, sgnd, stpc, stpf, stpq, gamma, theta;
    logical mcbound;

/*<       INTEGER INFOC >*/
/*<       DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX >*/
/*<       LOGICAL BRACKT,MCBOUND >*/

/*     SUBROUTINE MCSTEP */

/*     THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR */
/*     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR */
/*     A MINIMIZER OF THE FUNCTION. */

/*     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION */
/*     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS */
/*     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE */
/*     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A */
/*     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY */
/*     WITH ENDPOINTS STX AND STY. */

/*     THE SUBROUTINE STATEMENT IS */

/*       SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, */