lbfgs.c

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

/*<       RETURN >*/
    return 0;
/*<  196  IFLAG= -3 >*/
L196:
    *iflag = -3;
/*<       IF(LP.GT.0) WRITE(LP,240) >*/
/*
 240  FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
     .       ' ARE NOT POSITIVE)')
*/

    if (lb3_1.lp > 0) {
        printf(" IFLAG= -3 IMPROPER INPUT PARAMETERS (N OR M"
               " ARE NOT POSITIVE)");
    }

/*     FORMATS */
/*     ------- */

/*<  2 >*/
/*<  2 >*/
/*<  2 >*/
/*<  2 >*/
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* lbfgs_ */


/*     LAST LINE OF SUBROUTINE LBFGS */

static void write50(double* v, int n)
{
  int cols = 15;
  double vmax = 0;
  int i;
  double vmaxscale;
  for (i = 0; i < n; ++i)
    if (fabs(v[i]) > vmax)
      vmax = v[i];
  vmaxscale = log(fabs(vmax)) / log(10.0);
  vmaxscale = pow(10.0, ceil(vmaxscale) - 1);
  if (vmaxscale != 1.0)
    printf("  %e x\n", vmaxscale);

  for (i = 0; i < n; ++i) {
    if (i > 0 && i%cols == 0)
      printf("\n");
    printf(" %10.5f", v[i] / vmaxscale);
  }
  printf("\n");
}

/*<    >*/
/* Subroutine */ int lb1_(
  integer *iprint, integer *n, integer *m, doublereal *x, doublereal *f,
  doublereal *g, v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg
  )
{
/*     ------------------------------------------------------------- */
/*     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND */
/*     AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT. */
/*     ------------------------------------------------------------- */

/*<       INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M >*/
/*<       DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX >*/
/*<       LOGICAL FINISH >*/
/*<       COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX >*/

/*<       IF (ITER.EQ.0)THEN >*/
    /* Parameter adjustments */
    --iprint;
    --g;
    --x;

    /* Function Body */
    if (iter == 0) {
/*<            WRITE(MP,10) >*/
/*
  10   FORMAT('*************************************************')
*/
        printf("*************************************************\n");
/*<            WRITE(MP,20) N,M >*/
/*
 20   FORMAT('  N=',I5,'   NUMBER OF CORRECTIONS=',I2,
     .       /,  '       INITIAL VALUES')
*/
        printf("  N=%ld   NUMBER OF CORRECTIONS=%ld"
               "       INITIAL VALUES", *n, *m);
/*<            WRITE(MP,30)F,GNORM >*/
/*
 30   FORMAT(' F= ',1PD10.3,'   GNORM= ',1PD10.3)
*/
        printf(" F= %g   GNORM= %g\n", *f, gnorm);
/*<                  IF (IPRINT(2).GE.1)THEN >*/
        if (iprint[2] >= 1) {
/*<                      WRITE(MP,40) >*/
/*
 40   FORMAT(' VECTOR X= ')
*/
            printf(" VECTOR X= ");
/*<                      WRITE(MP,50) (X(I),I=1,N) >*/
/*
 50   FORMAT(6(2X,1PD10.3))
*/
            write50(x, *n);
/*<                      WRITE(MP,60) >*/
/*
 60   FORMAT(' GRADIENT VECTOR G= ')
*/
            printf(" GRADIENT VECTOR G= ");
/*<                      WRITE(MP,50) (G(I),I=1,N) >*/
/*
 50   FORMAT(6(2X,1PD10.3))
*/
            write50(g, *n);
/*<                   ENDIF >*/
        }
/*<            WRITE(MP,10) >*/
/*
  10   FORMAT('*************************************************')
*/
        printf("*************************************************\n");
/*<            WRITE(MP,70) >*/
/*
 70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
*/
        printf("   I   NFN    FUNC        GNORM       STEPLENGTH\n");
/*<       ELSE >*/
    } else {
/*<           IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN >*/
        if (iprint[1] == 0 && (iter != 1 && ! (finish))) {
            return 0;
        }
/*<               IF (IPRINT(1).NE.0)THEN >*/
        if (iprint[1] != 0) {
/*<                    IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN >*/
            if ((iter - 1) % iprint[1] == 0 || finish) {
/*<                          IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70) >*/
/*
 70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
*/
                if (iprint[2] > 1 && iter > 1) {
                    printf("   I   NFN    FUNC        GNORM       STEPLENGTH\n");
                }
/*<                          WRITE(MP,80)ITER,NFUN,F,GNORM,STP >*/
/*
 80   FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
*/
                printf("%4ld %4ld    %10.3f  %10.3f  %10.3f\n", iter, nfun, *f, gnorm, stp);
/*<                    ELSE >*/
            } else {
/*<                          RETURN >*/
                return 0;
/*<                    ENDIF >*/
            }
/*<               ELSE >*/
        } else {
/*<                    IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70) >*/
/*
 70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
*/
            if (iprint[2] > 1 && finish) {
                    printf("   I   NFN    FUNC        GNORM       STEPLENGTH\n");
            }
/*<                    WRITE(MP,80)ITER,NFUN,F,GNORM,STP >*/
/*
 80   FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
*/
            printf("%4ld %4ld    %10.3f  %10.3f  %10.3f\n", iter, nfun, *f, gnorm, stp);
/*<               ENDIF >*/
        }
/*<               IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN >*/
        if (iprint[2] == 2 || iprint[2] == 3) {
/*<                     IF (FINISH)THEN >*/
            if (finish) {
/*<                         WRITE(MP,90) >*/
/*
  90   FORMAT(' FINAL POINT X= ')
*/
                printf(" FINAL POINT X= ");
/*<                     ELSE >*/
            } else {
/*<                         WRITE(MP,40) >*/
/*
 40   FORMAT(' VECTOR X= ')
*/
                printf(" VECTOR X= ");
/*<                     ENDIF >*/
            }
/*<                       WRITE(MP,50)(X(I),I=1,N) >*/
/*
 50   FORMAT(6(2X,1PD10.3))
*/
            write50(x, *n);
/*<                   IF (IPRINT(2).EQ.3)THEN >*/
            if (iprint[2] == 3) {
/*<                       WRITE(MP,60) >*/
/*
 60   FORMAT(' GRADIENT VECTOR G= ')
*/
                printf(" GRADIENT VECTOR G= ");
/*<                       WRITE(MP,50)(G(I),I=1,N) >*/
/*
 50   FORMAT(6(2X,1PD10.3))
*/
                write50(g, *n);
/*<                   ENDIF >*/
            }
/*<               ENDIF >*/
        }
/*<             IF (FINISH) WRITE(MP,100) >*/
/*
 100  FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',
     .       /' IFLAG = 0')
*/
        if (finish) {
            printf(" THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.\n");
        }
/*<       ENDIF >*/
    }

/*<  10   FORMAT('*************************************************') >*/
/*<  2 >*/
/*<  30   FORMAT(' F= ',1PD10.3,'   GNORM= ',1PD10.3) >*/
/*<  40   FORMAT(' VECTOR X= ') >*/
/*<  50   FORMAT(6(2X,1PD10.3)) >*/
/*<  60   FORMAT(' GRADIENT VECTOR G= ') >*/
/*<  70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/) >*/
/*<  80   FORMAT(2(I4,1X),3X,3(1PD10.3,2X)) >*/
/*<  90   FORMAT(' FINAL POINT X= ') >*/
/*<  1 >*/

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

/*     ****** */


/*   ---------------------------------------------------------- */
/*     DATA */
/*   ---------------------------------------------------------- */

/*<       INTEGER LP,MP >*/
/*<       DOUBLE PRECISION GTOL,STPMIN,STPMAX >*/
/*<       COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX >*/
/*<       DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,9.0D-01,1.0D-20,1.0D+20/ >*/
/*<       END >*/

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

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

/*<       SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA) >*/
/* Subroutine */ int mcsrch_(
  integer *n, doublereal *x, doublereal *f,
  doublereal *g, doublereal *s,
  doublereal *xtol, doublereal *wa,
  v3p_netlib_lbfgs_global_t* v3p_netlib_lbfgs_global_arg
  )
{
    /* Initialized data */

    static doublereal p5 = .5;  /* constant */
    static doublereal p66 = .66;  /* constant */
    static doublereal xtrapf = 4.;  /* constant */
    static doublereal zero = 0.;  /* constant */

    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Local variables */
    integer j;
    extern /* Subroutine */ int mcstep_(
      doublereal *, doublereal *, doublereal *, doublereal *,
      doublereal *, doublereal *, v3p_netlib_lbfgs_global_t*);

/*<       INTEGER N,MAXFEV,INFO,NFEV >*/
/*<       DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX >*/
/*<       DOUBLE PRECISION X(N),G(N),S(N),WA(N) >*/
/*<       COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX >*/
/*<       SAVE >*/

/*                     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 */