lbfgs.c

NIST Handwriting OCR Testbed

    r1 = -work[inmc];
    saxpy_(&n, &r1, &work[iycn + 1], &one, &work[1], &one);
  }

  for(i = 1; i <= n; i++)
    work[i] = diag[i] * work[i];

  for(i = 1; i <= bound; i++) {
    yr = sdot_(&n, &work[iypt + cp * n + 1], &one, &work[1], &one);
    beta = work[n + cp + 1] * yr;
    inmc = n + m + cp + 1;
    beta = work[inmc] - beta;
    iscn = ispt + cp * n;
    saxpy_(&n, &beta, &work[iscn + 1], &one, &work[1], &one);
    if(++cp == m)
      cp = 0;
  }

  /* -------------------------------
     Store the new search direction.
     ------------------------------- */
  for(i = 1; i <= n; i++)
    work[ispt + point * n + i] = work[i];

  /* -------------------------------------------------------------
     Obtain the one-dimensional minimizer of the function by using
     the line search routine mcsrch.
     ------------------------------------------------------------- */
 L165:
  nfev = 0;
  *stp = 1.;
  if(*iter == 1)
    *stp = stp1;
  for(i = 1; i <= n; i++)
    work[i] = g[i];

#ifdef SURVEY
  survey(n, &x[1], &work[ispt + point * n + 1], *stp);
  sprintf(str, " calling mcsrch with *stp %e\n\n", *stp);
  fsaso(str);
#endif

 L172:
  mcsrch(n, &x[1], f, &g[1], &work[ispt + point * n + 1], stp, ftol,
    gtol, xtol, stpmin, stpmax, maxfev, info, &nfev, &diag[1], fp_err);
  if(*info == -1) {
    *iflag = 1;
    return;
  }
  if(*info != 1)
    goto L190;
  nfun += nfev;

  /* -----------------------------------------
     Compute the new step and gradient change.
     ----------------------------------------- */
  npt = point * n;
  for(i = 1; i <= n; i++) {
    work[ispt + npt + i] = *stp * work[ispt + npt + i];
    work[iypt + npt + i] = g[i] - work[i];
  }
  if(++point == m)
    point = 0;

  /* -----------------
     Termination test.
     ---------------- */
  gnorm = snrm2_(&n, &g[1], &one);
  xnorm = snrm2_(&n, &x[1], &one);
  xnorm =  (int)max(1., xnorm);
  /* jlb change:  The original Fortran had
      IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
  here, but now the gnorm/xnorm termination is handled by
  optchk(). */

  if(iprint[1] >= 0)
    lb1(&iprint[1], *iter, nfun, gnorm, n, m, &x[1], f, &g[1],
      *stp, finish, fp_mon);
  optchk(FALSE, *itopt, &x[1], f, ierr, &pctmin_junk);
  if(*ierr != 0) {
    *iflag = 10;
    return;
  }
  if(finish) {
    *iflag = 0;
    return;
  }
  goto L80;

  /* -----------------------------------------
     End of main iteration loop.  Error exits.
     ----------------------------------------- */
 L190:
  *iflag = -1;
  if(fp_err != (FILE *)NULL)
    fprintf(fp_err, "\n *iflag == -1\n line search failed. see \
documentation of routine mcsrch\n error return of line search: \
*info = %d\n possible causes: function or gradient are incorrect or \
incorrect tolerances\n", *info);
  return;
 L195:
  *iflag = -2;
  if(fp_err != (FILE *)NULL)
    fprintf(fp_err, "\n *iflag == -2\n the %d-th diagonal element of \
the\n inverse Hessian approximation is not positive\n", i);
  return;
 L196:
  *iflag = -3;
  if(fp_err != (FILE *)NULL)
    fprintf(fp_err, "\n *iflag == -3\n improper input parameters \
(n or m is not positive)\n");
  return;
}

/* end of function lbfgs */

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

  /* ------------
     function lb1
     ------------

  Prints monitoring information.  The frequency and amount of output
  are controlled by iprint.

  Input args:
    iprint: Array of two ints, controlling printing.
    iter:
    nfun:
    gnorm:
    n:
    m:
    x:
    f:
    g:
    stp:
    finish:
    fp_mon;
*/

void
lb1(iprint, iter, nfun, gnorm, n, m, x, f, g, stp, finish, fp_mon)
int *iprint, iter, nfun, n, m, finish;
float gnorm, *x, f, *g, stp;
FILE *fp_mon;
{
  int i;
  static int one = 1;
  static float gbyx, xnorm;

  /* parameter adjustments */
  --g;
  --x;
  --iprint;

  /* jlb added next line (gtc changed use of sdot_ and sqrt, to a call
  of snrm2): */
  xnorm = snrm2_(&n, &x[1], &one);

  gbyx = 0.;
  if(xnorm > 0.)
    gbyx = gnorm / xnorm;
  if(iter == 0) {
    fprintf(fp_mon, "**********************************************\
***\n");
    fprintf(fp_mon, "  n = %d   number of corrections = %d\n\
       initial values:\n", n, m);
    fprintf(fp_mon, " f = %10.3e   gnorm = %10.3e\n", f, gnorm);
    if(iprint[2] >= 1) {
      fprintf(fp_mon, " vector x = ");
      for(i = 0; i < n; i++)
	fprintf(fp_mon, "  %10.3e", x[i]);
      fprintf(fp_mon, "\n");
      fprintf(fp_mon, " gradient vector g = ");
      for(i = 0; i < n; i++)
	fprintf(fp_mon, "  %10.3e", g[i]);
      fprintf(fp_mon, "\n");
    }
    fprintf(fp_mon, "**********************************************\
***\n");
    fprintf(fp_mon, "\n   i   nfn    func        gnorm       \
steplength\n");
  }
  else {
    if(iprint[1] == 0 && (iter != 1 && !finish))
      return;
    if(iprint[1] != 0) {
      if((iter - 1) % iprint[1] == 0 || finish) {
	if(iprint[2] > 1 && iter > 1)
	  fprintf(fp_mon, "\n   i   nfn    func        gnorm       \
steplength\n");
	fprintf(fp_mon, "%d %d    %10.3e  %10.3e  %10.3e  %10.3e\n",
          iter, nfun, f, gnorm, stp, gbyx);
      }
      else
	return;
    }
    else {
      if(iprint[2] > 1 && finish)
	fprintf(fp_mon, "\n   i   nfn    func        gnorm       \
steplength\n");
      fprintf(fp_mon, "%d %d    %10.3e  %10.3e  %10.3e  %10.3e\n",
          iter, nfun, f, gnorm, stp, gbyx);
    }      
    if(iprint[2] == 2 || iprint[2] == 3) {
      if(finish)
	fprintf(fp_mon, " final point x = ");
      else
	fprintf(fp_mon, " vector x = ");
      for(i = 0; i < n; i++)
	fprintf(fp_mon, "  %10.3e", x[i]);
      fprintf(fp_mon, "\n");
      if(iprint[2] == 3) {
	fprintf(fp_mon, " gradient vector g = ");
	for(i = 0; i < n; i++)
	  fprintf(fp_mon, "  %10.3e", g[i]);
	fprintf(fp_mon, "\n");
      }
    }
    if(finish)
      fprintf(fp_mon, "\n the minimization terminated without \
detecting errors.\n iflag = 0\n");
  }
  return;
}

/* end of function lb1 */

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

  /* ---------------
     function 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 mcsrch 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) <= f(x) + ftol*stp*(gradf(x)'s),

  and the curvature condition

        abs(gradf(x+stp*s)'s)) <= 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 calling statement is

     mcsrch(n, x, f, g, s, stp, ftol, gtol, xtol, stpmin, stpmax,
       maxfev, info, nfev, wa, fp_err);

  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.
      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.

    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.

      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.

    fp_err is a FILE pointer to which to write error messages.
    If it is (FILE *)NULL, then no error messages are produced.

  Subprograms called:

    mcstep (source in this file)

  Argonne National Laboratory.  MINPACK Project.  June 1983.
  Jorge J. More', David J. Thuente. */

void
mcsrch(n, x, f, g, s, stp, ftol, gtol, xtol, stpmin, stpmax, maxfev,
  info, nfev, wa, fp_err)
int n, maxfev, *info, *nfev;
float *x, f, *g, *s, *stp, ftol, gtol, xtol, stpmin, stpmax, *wa;
FILE *fp_err;
{
  static int j, infoc, stage1, brackt, one = 1;
  /* float r1; */
  static float dgxm, dgym, finit, width, stmin, stmax, width1, ftest1,
    dg, fm, gs, fx, fy, dginit, dgtest, dgm, dgx, dgy, fxm, fym, stx,
    sty;
  void mcstep();

#define XTRAPF 4.

  /* parameter adjustments */
  --wa;
  --s;
  --g;
  --x;

  if(*info == -1)
    goto L45;
  gs = sdot_(&n, &g[1], &one, &s[1], &one) / (snrm2_(&n, &g[1],
    &one) * snrm2_(&n, &s[1], &one));
  infoc = 1;

  /* Check the input parameters for errors. */
  if(n <= 0 || *stp <= 0. || ftol < 0. || gtol < 0. || xtol < 0. ||
    stpmin < 0. || stpmax < stpmin || maxfev <= 0)
    return;

  /* Compute the initial gradient in the search direction and check
  that s is a descent direction. */
  dginit = 0.;
  for(j = 1; j <= n; j++)
    dginit += g[j] * s[j];
  if(dginit >= 0.) {
    if(fp_err != (FILE *)NULL)
      fprintf(fp_err, "\n  the search direction is not a descent \
direction\n");
    return;
  }

  /* Initialize local variables. */
  brackt = FALSE;
  stage1 = TRUE;
  *nfev = 0;
  finit = f;
  dgtest = ftol * dginit;