lbfgs.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,199 行 · 第 1/3 页

        if (cp == -1) {
            cp = *m - 1;
        }
        sq = ddot_(n, &w[ispt + cp * *n], &c__1, w, &c__1);
        inmc = *n + *m + cp;
        iycn = iypt + cp * *n;
        w[inmc] = w[*n + cp] * sq;
        d__1 = -w[inmc];
        daxpy_(n, &d__1, &w[iycn], &c__1, w, &c__1);
    }

    for (i = 0; i < *n; ++i) {
        w[i] *= diag[i];
    }

    for (i = 0; i < bound; ++i) {
        yr = ddot_(n, &w[iypt + cp * *n], &c__1, w, &c__1);
        beta = w[*n + cp] * yr;
        inmc = *n + *m + cp;
        beta = w[inmc] - beta;
        iscn = ispt + cp * *n;
        daxpy_(n, &beta, &w[iscn], &c__1, w, &c__1);
        ++cp;
        if (cp == *m) {
            cp = 0;
        }
    }

/*     STORE THE NEW SEARCH DIRECTION */
/*     ------------------------------ */

    for (i = 0; i < *n; ++i) {
        w[ispt + point * *n + i] = w[i];
    }

/*     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION */
/*     BY USING THE LINE SEARCH ROUTINE MCSRCH */
/*     ---------------------------------------------------- */
L165:
    nfev = 0;
/* awf changed initial step from ONE to be parametrized. */
    stp = lb3_1.stpawf;
    if (iter == 1) {
        stp = stp1;
    }
    for (i = 0; i < *n; ++i) {
        w[i] = g[i];
    }
L172:
    mcsrch_(n, x, f, g, &w[ispt + point * *n], &stp, &ftol, xtol, &maxfev, &info, &nfev, diag);
    if (info == -1) {
        *iflag = 1;
/*       Return, in order to get another sample of F and G. */
/*       Next call comes right back here. */
        return;
    }
    if (info != 1) {
        goto L190;
    }
    nfun += nfev;

/*     COMPUTE THE NEW STEP AND GRADIENT CHANGE */
/*     ----------------------------------------- */

    npt = point * *n;
    for (i = 0; i < *n; ++i) {
        w[ispt + npt + i] *= stp;
        w[iypt + npt + i] = g[i] - w[i];
    }
    ++point;
    if (point == *m) {
        point = 0;
    }

/*     TERMINATION TEST */
/*     ---------------- */

    gnorm = sqrt(ddot_(n, g, &c__1, g, &c__1));
    xnorm = sqrt(ddot_(n, x, &c__1, x, &c__1));
    xnorm = max(1.,xnorm);
    if (gnorm / xnorm <= *eps) {
        finish = TRUE_;
    }

    if (iprint[0] >= 0) {
        lb1_(iprint, &iter, &nfun, &gnorm, n, m, x, f, g, &stp, &finish);
    }
    if (finish) {
        *iflag = 0;
        return;
    }
    goto L80;

/*     ------------------------------------------------------------ */
/*     END OF MAIN ITERATION LOOP. ERROR EXITS. */
/*     ------------------------------------------------------------ */

L190:
    *iflag = -1;
    if (lb3_1.lp > 0) {
        lbptf_("IFLAG= -1. LINE SEARCH FAILED.");
        lbptf_(" SEE DOCUMENTATION OF ROUTINE MCSRCH");
        lbp1d_(" ERROR RETURN  OF LINE SEARCH: INFO=%d", &info);
        lbptf_(" POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE ");
        lbptf_(" INCORRECT OR INCORRECT TOLERANCES");
    }
    return;
L195:
    *iflag = -2;
    if (lb3_1.lp > 0) {
        lbp1d_("IFLAG=-2, THE %d-TH DIAGONAL ELEMENT OF THE", &i);
        lbptf_("INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE");
    }
    return;
L196:
    *iflag = -3;
    if (lb3_1.lp > 0) {
        lbptf_("IFLAG= -3, IMPROPER INPUT PARAMETERS.");
        lbptf_(" (N OR M ARE NOT POSITIVE)");
    }
    return;
} /* lbfgs_ */


/*     LAST LINE OF SUBROUTINE LBFGS */


/*      SUBROUTINE LB1(IPRINT,ITER,NFUN,GNORM,N,M,X,F,G,STP,FINISH) */
/*      ** moved to c file */

/*   ---------------------------------------------------------- */

/*   These routines removed for insertion into TargetJr netlib */

/* awf       subroutine daxpy(n,da,dx,incx,dy,incy) */
/* awf c */
/* awf c     constant times a vector plus a vector. */
/* awf c     uses unrolled loops for increments equal to one. */
/* awf c     jack dongarra, linpack, 3/11/78. */
/* awf c */
/* awf       double precision dx(1),dy(1),da */
/* awf       integer i,incx,incy,ix,iy,m,mp1,n */
/* awf c */
/* awf       if(n.le.0)return */
/* awf       if (da .eq. 0.0d0) return */
/* awf       if(incx.eq.1.and.incy.eq.1)go to 20 */
/* awf c */
/* awf c        code for unequal increments or equal increments */
/* awf c          not equal to 1 */
/* awf c */
/* awf       ix = 1 */
/* awf       iy = 1 */
/* awf       if(incx.lt.0)ix = (-n+1)*incx + 1 */
/* awf       if(incy.lt.0)iy = (-n+1)*incy + 1 */
/* awf       do 10 i = 1,n */
/* awf         dy(iy) = dy(iy) + da*dx(ix) */
/* awf         ix = ix + incx */
/* awf         iy = iy + incy */
/* awf    10 continue */
/* awf       return */
/* awf c */
/* awf c        code for both increments equal to 1 */
/* awf c */
/* awf c */
/* awf c        clean-up loop */
/* awf c */
/* awf    20 m = mod(n,4) */
/* awf       if( m .eq. 0 ) go to 40 */
/* awf       do 30 i = 1,m */
/* awf         dy(i) = dy(i) + da*dx(i) */
/* awf    30 continue */
/* awf       if( n .lt. 4 ) return */
/* awf    40 mp1 = m + 1 */
/* awf       do 50 i = mp1,n,4 */
/* awf         dy(i) = dy(i) + da*dx(i) */
/* awf         dy(i + 1) = dy(i + 1) + da*dx(i + 1) */
/* awf         dy(i + 2) = dy(i + 2) + da*dx(i + 2) */
/* awf         dy(i + 3) = dy(i + 3) + da*dx(i + 3) */
/* awf    50 continue */
/* awf       return */
/* awf       end */
/* awf C */
/* awf C */
/* awf C   ---------------------------------------------------------- */
/* awf C */
/* awf       double precision function ddot(n,dx,incx,dy,incy) */
/* awf c */
/* awf c     forms the dot product of two vectors. */
/* awf c     uses unrolled loops for increments equal to one. */
/* awf c     jack dongarra, linpack, 3/11/78. */
/* awf c */
/* awf       double precision dx(1),dy(1),dtemp */
/* awf       integer i,incx,incy,ix,iy,m,mp1,n */
/* awf c */
/* awf       ddot = 0.0d0 */
/* awf       dtemp = 0.0d0 */
/* awf       if(n.le.0)return */
/* awf       if(incx.eq.1.and.incy.eq.1)go to 20 */
/* awf c */
/* awf c        code for unequal increments or equal increments */
/* awf c          not equal to 1 */
/* awf c */
/* awf       ix = 1 */
/* awf       iy = 1 */
/* awf       if(incx.lt.0)ix = (-n+1)*incx + 1 */
/* awf       if(incy.lt.0)iy = (-n+1)*incy + 1 */
/* awf       do 10 i = 1,n */
/* awf         dtemp = dtemp + dx(ix)*dy(iy) */
/* awf         ix = ix + incx */
/* awf         iy = iy + incy */
/* awf    10 continue */
/* awf       ddot = dtemp */
/* awf       return */
/* awf c */
/* awf c        code for both increments equal to 1 */
/* awf c */
/* awf c */
/* awf c        clean-up loop */
/* awf c */
/* awf    20 m = mod(n,5) */
/* awf       if( m .eq. 0 ) go to 40 */
/* awf       do 30 i = 1,m */
/* awf         dtemp = dtemp + dx(i)*dy(i) */
/* awf    30 continue */
/* awf       if( n .lt. 5 ) go to 60 */
/* awf    40 mp1 = m + 1 */
/* awf       do 50 i = mp1,n,5 */
/* awf         dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + */
/* awf      *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) */
/* awf    50 continue */
/* awf    60 ddot = dtemp */
/* awf       return */
/* awf       end */
/*    ------------------------------------------------------------------ */

/*     ************************** */
/*     LINE SEARCH ROUTINE MCSRCH */
/*     ************************** */

/* Subroutine */ void mcsrch_(n, x, f, g, s, stp, ftol, xtol, maxfev, info, nfev, wa)
integer *n;
doublereal *x, *f, *g, *s, *stp, *ftol, *xtol;
integer *maxfev, *info, *nfev;
doublereal *wa;
{
    /* Initialized data */

    static doublereal xtrapf = 4.;

    /* Local variables */
    static doublereal dgxm, dgym;
    static integer j, infoc;
    static doublereal finit, width, stmin, stmax;
    static logical stage1;
    static doublereal width1, ftest1, dg, fm, fx, fy;
    static logical brackt;
    static doublereal dginit, dgtest;
    static doublereal dgm, dgx, dgy, fxm, fym, stx, sty;


/*                     SUBROUTINE MCSRCH */

/*     A slight modification of the subroutine CSRCH of More' and Thuente. */
/*     The changes are to allow reverse communication, and do not affect */
/*     the performance of the routine. */

/*     THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES */
/*     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION. */

/*     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF */
/*     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF */
/*     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A */
/*     MINIMIZER OF THE MODIFIED FUNCTION */

/*          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S). */

/*     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION */
/*     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, */
/*     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT */
/*     CONTAINS A MINIMIZER OF F(X+STP*S). */

/*     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES */
/*     THE SUFFICIENT DECREASE CONDITION */

/*           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S), */

/*     AND THE CURVATURE CONDITION */

/*           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S). */

/*     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION */
/*     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES */
/*     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH */
/*     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING */
/*     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY */
/*     SATISFIES THE SUFFICIENT DECREASE CONDITION. */

/*     THE SUBROUTINE STATEMENT IS */

/*        SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA) */
/*     WHERE */

/*       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER */
/*         OF VARIABLES. */

/*       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE */
/*         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS */
/*         X + STP*S. */

/*       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F */
/*         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S. */

/*       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE */
/*         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT */
/*         OF F AT X + STP*S. */

/*       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE */
/*         SEARCH DIRECTION. */

/*       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN */
/*         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT */
/*         STP CONTAINS THE FINAL ESTIMATE. */

/*       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse */
/*         communication implementation GTOL is defined in a COMMON */
/*         statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE */
/*         CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE */
/*         SATISFIED. */

/*       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS */
/*         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY */
/*         IS AT MOST XTOL. */

/*       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH */
/*         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse */
/*         communication implementatin they are defined in a COMMON */
/*         statement). */

/*       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION */
/*         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST */
/*         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.  */
/*       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 */

/*     ********** */

    if (*info == -1) {
        goto L45;
    }
    infoc = 1;

/*     CHECK THE INPUT PARAMETERS FOR ERRORS. */

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

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

    dginit = 0.;
    for (j = 0; j < *n; ++j) {
        dginit += g[j] * s[j];
    }
    if (dginit >= 0.) {
        lbptf_("THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION");
        return;
    }

/*     INITIALIZE LOCAL VARIABLES. */

    brackt = FALSE_;
    stage1 = TRUE_;
    *nfev = 0;
    finit = *f;
    dgtest = *ftol * dginit;
    width = lb3_1.stpmax - lb3_1.stpmin;
    width1 = width / .5;
    for (j = 0; j < *n; ++j) {
        wa[j] = x[j];
    }