lbfgs.c

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

    }
/*<   10  ITER= 0 >*/
L10:
    iter = 0;
/*<       IF(N.LE.0.OR.M.LE.0) GO TO 196 >*/
    if (*n <= 0 || *m <= 0) {
        goto L196;
    }
/*<       IF(GTOL.LE.1.D-04) THEN >*/
    if (lb3_1.gtol <= 1e-4) {
/*<         IF(LP.GT.0) WRITE(LP,245) >*/
/*
 245  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',
     .       / ' IT HAS BEEN RESET TO 9.D-01')
*/
        if (lb3_1.lp > 0) {
            printf("  GTOL IS LESS THAN OR EQUAL TO 1.D-04");
            printf("  IT HAS BEEN RESET TO 9.D-01");
        }
/*<         GTOL=9.D-01 >*/
        lb3_1.gtol = .9;
/*<       ENDIF >*/
    }
/*<       NFUN= 1 >*/
    nfun = 1;
/*<       POINT= 0 >*/
    point = 0;
/*<       FINISH= .FALSE. >*/
    finish = FALSE_;
/*<       IF(DIAGCO) THEN >*/
    if (*diagco) {
/*<          DO 30 I=1,N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<  30      IF (DIAG(I).LE.ZERO) GO TO 195 >*/
/* L30: */
            if (diag[i__] <= zero) {
                goto L195;
            }
        }
/*<       ELSE >*/
    } else {
/*<          DO 40 I=1,N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<  40      DIAG(I)= 1.0D0 >*/
/* L40: */
            diag[i__] = 1.;
        }
/*<       ENDIF >*/
    }

/*     THE WORK VECTOR W IS DIVIDED AS FOLLOWS: */
/*     --------------------------------------- */
/*     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND */
/*         OTHER TEMPORARY INFORMATION. */
/*     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO. */
/*     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED */
/*         IN THE FORMULA THAT COMPUTES H*G. */
/*     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH */
/*         STEPS. */
/*     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M */
/*         GRADIENT DIFFERENCES. */

/*     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A */
/*     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT. */

/*<       ISPT= N+2*M >*/
    ispt = *n + (*m << 1);
/*<       IYPT= ISPT+N*M      >*/
    iypt = ispt + *n * *m;
/*<       DO 50 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<  50   W(ISPT+I)= -G(I)*DIAG(I) >*/
/* L50: */
        w[ispt + i__] = -g[i__] * diag[i__];
    }
/*<       GNORM= DSQRT(DDOT(N,G,1,G,1)) >*/
    gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
/*<       STP1= ONE/GNORM >*/
    stp1 = one / gnorm;

/*     PARAMETERS FOR LINE SEARCH ROUTINE */

/*<       FTOL= 1.0D-4 >*/
    ftol = 1e-4;
/*<       MAXFEV= 20 >*/
    maxfev = 20;

/*<    >*/
    if (iprint[1] >= 0) {
        lb1_(&iprint[1], n, m, &x[1], f, &g[1], v3p_netlib_lbfgs_global_arg);
    }

/*    -------------------- */
/*     MAIN ITERATION LOOP */
/*    -------------------- */

/*<  80   ITER= ITER+1 >*/
L80:
    ++iter;
/*<       INFO=0 >*/
    info = 0;
/*<       BOUND=ITER-1 >*/
    bound = iter - 1;
/*<       IF(ITER.EQ.1) GO TO 165 >*/
    if (iter == 1) {
        goto L165;
    }
/*<       IF (ITER .GT. M)BOUND=M >*/
    if (iter > *m) {
        bound = *m;
    }

/*<          YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1) >*/
    ys = ddot_(n, &w[iypt + npt + 1], &c__1, &w[ispt + npt + 1], &c__1);
/*<       IF(.NOT.DIAGCO) THEN >*/
    if (! (*diagco)) {
/*<          YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1) >*/
        yy = ddot_(n, &w[iypt + npt + 1], &c__1, &w[iypt + npt + 1], &c__1);
/*<          DO 90 I=1,N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<    90    DIAG(I)= YS/YY >*/
/* L90: */
            diag[i__] = ys / yy;
        }
/*<       ELSE >*/
    } else {
/*<          IFLAG=2 >*/
        *iflag = 2;
/*<          RETURN >*/
        return 0;
/*<       ENDIF >*/
    }
/*<  100  CONTINUE >*/
L100:
/*<       IF(DIAGCO) THEN >*/
    if (*diagco) {
/*<         DO 110 I=1,N >*/
        i__1 = *n;
        for (i__ = 1; i__ <= i__1; ++i__) {
/*<  110    IF (DIAG(I).LE.ZERO) GO TO 195 >*/
/* L110: */
            if (diag[i__] <= zero) {
                goto L195;
            }
        }
/*<       ENDIF >*/
    }

/*     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, */
/*     "Updating quasi-Newton matrices with limited storage", */
/*     Mathematics of Computation, Vol.24, No.151, pp. 773-782. */
/*     --------------------------------------------------------- */

/*<       CP= POINT >*/
    cp = point;
/*<       IF (POINT.EQ.0) CP=M >*/
    if (point == 0) {
        cp = *m;
    }
/*<       W(N+CP)= ONE/YS >*/
    w[*n + cp] = one / ys;
/*<       DO 112 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<  112  W(I)= -G(I) >*/
/* L112: */
        w[i__] = -g[i__];
    }
/*<       CP= POINT >*/
    cp = point;
/*<       DO 125 I= 1,BOUND >*/
    i__1 = bound;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          CP=CP-1 >*/
        --cp;
/*<          IF (CP.EQ. -1)CP=M-1 >*/
        if (cp == -1) {
            cp = *m - 1;
        }
/*<          SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1) >*/
        sq = ddot_(n, &w[ispt + cp * *n + 1], &c__1, &w[1], &c__1);
/*<          INMC=N+M+CP+1 >*/
        inmc = *n + *m + cp + 1;
/*<          IYCN=IYPT+CP*N >*/
        iycn = iypt + cp * *n;
/*<          W(INMC)= W(N+CP+1)*SQ >*/
        w[inmc] = w[*n + cp + 1] * sq;
/*<          CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1) >*/
        d__1 = -w[inmc];
        daxpy_(n, &d__1, &w[iycn + 1], &c__1, &w[1], &c__1);
/*<  125  CONTINUE >*/
/* L125: */
    }

/*<       DO 130 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<  130  W(I)=DIAG(I)*W(I) >*/
/* L130: */
        w[i__] = diag[i__] * w[i__];
    }

/*<       DO 145 I=1,BOUND >*/
    i__1 = bound;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          YR= DDOT(N,W(IYPT+CP*N+1),1,W,1) >*/
        yr = ddot_(n, &w[iypt + cp * *n + 1], &c__1, &w[1], &c__1);
/*<          BETA= W(N+CP+1)*YR >*/
        beta = w[*n + cp + 1] * yr;
/*<          INMC=N+M+CP+1 >*/
        inmc = *n + *m + cp + 1;
/*<          BETA= W(INMC)-BETA >*/
        beta = w[inmc] - beta;
/*<          ISCN=ISPT+CP*N >*/
        iscn = ispt + cp * *n;
/*<          CALL DAXPY(N,BETA,W(ISCN+1),1,W,1) >*/
        daxpy_(n, &beta, &w[iscn + 1], &c__1, &w[1], &c__1);
/*<          CP=CP+1 >*/
        ++cp;
/*<          IF (CP.EQ.M)CP=0 >*/
        if (cp == *m) {
            cp = 0;
        }
/*<  145  CONTINUE >*/
/* L145: */
    }

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

/*<        DO 160 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<  160   W(ISPT+POINT*N+I)= W(I) >*/
/* L160: */
        w[ispt + point * *n + i__] = w[i__];
    }

/*     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION */
/*     BY USING THE LINE SEARCH ROUTINE MCSRCH */
/*     ---------------------------------------------------- */
/*<  165  NFEV=0 >*/
L165:
    nfev = 0;
/*<       STP=ONE >*/
/* awf changed initial step from ONE to be parametrized. */
    stp = lb3_1.stpinit;
/*<       IF (ITER.EQ.1) STP=STP1 >*/
    if (iter == 1) {
        stp = stp1;
    }
/*<       DO 170 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<  170  W(I)=G(I) >*/
/* L170: */
        w[i__] = g[i__];
    }
/*<  172  CONTINUE >*/
L172:
/*<    >*/
    mcsrch_(n, &x[1], f, &g[1], &w[ispt + point * *n + 1], xtol,
            &diag[1], v3p_netlib_lbfgs_global_arg);
/*<       IF (INFO .EQ. -1) THEN >*/
    if (info == -1) {
/*<         IFLAG=1 >*/
        *iflag = 1;
/*<         RETURN >*/
        return 0;
/*<       ENDIF >*/
    }
/*<       IF (INFO .NE. 1) GO TO 190 >*/
    if (info != 1) {
        goto L190;
    }
/*<       NFUN= NFUN + NFEV >*/
    nfun += nfev;

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

/*<       NPT=POINT*N >*/
    npt = point * *n;
/*<       DO 175 I=1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<       W(ISPT+NPT+I)= STP*W(ISPT+NPT+I) >*/
        w[ispt + npt + i__] = stp * w[ispt + npt + i__];
/*<  175  W(IYPT+NPT+I)= G(I)-W(I) >*/
/* L175: */
        w[iypt + npt + i__] = g[i__] - w[i__];
    }
/*<       POINT=POINT+1 >*/
    ++point;
/*<       IF (POINT.EQ.M)POINT=0 >*/
    if (point == *m) {
        point = 0;
    }

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

/*<       GNORM= DSQRT(DDOT(N,G,1,G,1)) >*/
    gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
/*<       XNORM= DSQRT(DDOT(N,X,1,X,1)) >*/
    xnorm = sqrt(ddot_(n, &x[1], &c__1, &x[1], &c__1));
/*<       XNORM= DMAX1(1.0D0,XNORM) >*/
    xnorm = max(1.,xnorm);
/*<       IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE. >*/
    if (gnorm / xnorm <= *eps) {
        finish = TRUE_;
    }

/*<    >*/
    if (iprint[1] >= 0) {
        lb1_(&iprint[1], n, m, &x[1], f, &g[1], v3p_netlib_lbfgs_global_arg);
    }
/*<       IF (FINISH) THEN >*/
    if (finish) {
/*<          IFLAG=0 >*/
        *iflag = 0;
/*<          RETURN >*/
        return 0;
/*<       ENDIF >*/
    }
/*<       GO TO 80 >*/
    goto L80;

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

/*<  190  IFLAG=-1 >*/
L190:
    *iflag = -1;
/*<       IF(LP.GT.0) WRITE(LP,200) INFO >*/
/*
 200  FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'
     .          ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'
     .          ' OF LINE SEARCH: INFO= ',I2,/
     .          ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
     .          ' OR INCORRECT TOLERANCES')
*/
    if (lb3_1.lp > 0) {
        printf(" IFLAG= -1  LINE SEARCH FAILED. SEE"
               " DOCUMENTATION OF ROUTINE MCSRCH ERROR RETURN"
               " OF LINE SEARCH: INFO= %ld"
               " POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT"
               " OR INCORRECT TOLERANCES", info);
    }
/*<       RETURN >*/
    return 0;
/*<  195  IFLAG=-2 >*/
L195:
    *iflag = -2;
/*<       IF(LP.GT.0) WRITE(LP,235) I >*/
/*
 235  FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
     .       ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
*/
    if (lb3_1.lp > 0) {
        printf("IFLAG= -2 THE %ld-TH DIAGONAL ELEMENT OF THE"
               " INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE", i__);
    }