lbfgs.c

Conditional Random Fields的训练识别工具

    dtemp += dx[ix] * dy[iy];
    ix += *incx;
    iy += *incy;
    /* L10: */
  }
  ret_val = dtemp;
  return ret_val;

  /*        code for both increments equal to 1 */


  /*        clean-up loop */

 L20:
  m = *n % 5;
  if (m == 0) {
    goto L40;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__) {
    dtemp += dx[i__] * dy[i__];
    /* L30: */
  }
  if (*n < 5) {
    goto L60;
  }
 L40:
  mp1 = m + 1;
  i__1 = *n;
  for (i__ = mp1; i__ <= i__1; i__ += 5) {
    dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
                                                                       i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
                                                                                                                               4] * dy[i__ + 4];
    /* L50: */
  }
 L60:
  ret_val = dtemp;
  return ret_val;
} /* ddot_ */

/* Subroutine */ static int 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 p5 = .5;
  static doublereal p66 = .66;
  static doublereal xtrapf = 4.;
  static doublereal zero = 0.;

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

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

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

  /* Function Body */
  if (*info == -1) {
    goto L45;
  }
  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;
  i__1 = *n;
  for (j = 1; j <= i__1; ++j) {
    dginit += g[j] * s[j];
    /* L10: */
  }
  if (dginit >= zero) {
    return 0;
  }

  /*     INITIALIZE LOCAL VARIABLES. */

  brackt = FALSE_;
  stage1 = TRUE_;
  *nfev = 0;
  finit = *f;
  dgtest = *ftol * dginit;
  width = lb3_1.stpmax - lb3_1.stpmin;
  width1 = width / p5;
  i__1 = *n;
  for (j = 1; j <= i__1; ++j) {
    wa[j] = x[j];
    /* 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;
  fx = finit;
  dgx = dginit;
  sty = zero;
  fy = finit;
  dgy = dginit;

  /*     START OF ITERATION. */

 L30:

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

  if (brackt) {
    stmin = min(stx,sty);
    stmax = max(stx,sty);
  } else {
    stmin = stx;
    stmax = *stp + xtrapf *(*stp - stx);
  }

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

  *stp = max(*stp,lb3_1.stpmin);
  *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. */

  i__1 = *n;
  for (j = 1; j <= i__1; ++j) {
    x[j] = wa[j] + *stp * s[j];
    /* L40: */
  }
  *info = -1;
  return 0;

 L45:
  *info = 0;
  ++(*nfev);
  dg = zero;
  i__1 = *n;
  for (j = 1; j <= i__1; ++j) {
    dg += g[j] * s[j];
    /* L50: */
  }
  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 >= *maxfev) {
    *info = 3;
  }
  if (brackt && stmax - stmin <= *xtol * stmax) {
    *info = 2;
  }
  if (*f <= ftest1 && abs(dg) <= lb3_1.gtol * (-dginit)) {
    *info = 1;
  }

  /*        CHECK FOR TERMINATION. */

  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 && *f <= fx && *f > ftest1) {

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

    fm = *f - *stp * dgtest;
    fxm = fx - stx * dgtest;
    fym = fy - sty * dgtest;
    dgm = dg - dgtest;
    dgxm = dgx - dgtest;
    dgym = dgy - dgtest;

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

    mcstep_(&stx, &fxm, &dgxm, &sty, &fym, &dgym, stp, &fm, &dgm, &brackt,
            &stmin, &stmax, &infoc);

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

    fx = fxm + stx * dgtest;
    fy = fym + sty * dgtest;
    dgx = dgxm + dgtest;
    dgy = dgym + dgtest;
  } else {

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

    mcstep_(&stx, &fx, &dgx, &sty, &fy, &dgy, stp, f, &dg, &brackt, &
            stmin, &stmax, &infoc);
  }

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

  if (brackt) {
    if ((d__1 = sty - stx, abs(d__1)) >= p66 * width1) {
      *stp = stx + p5 *(sty - stx);
    }
    width1 = width;
    width =(d__1 = sty - stx, abs(d__1));
  }

  /*        END OF ITERATION. */

  goto L30;

  /*     LAST LINE OF SUBROUTINE MCSRCH. */

} /* mcsrch_ */

/* Subroutine */ static int mcstep_(stx, fx, dx, sty, fy, dy, stp, fp, dp, brackt,
                                    stpmin, stpmax, info)
     doublereal *stx, *fx, *dx, *sty, *fy, *dy, *stp, *fp, *dp;
     logical *brackt;
     doublereal *stpmin, *stpmax;
     integer *info;
{
  /* System generated locals */
  doublereal d__1, d__2, d__3;

  /* Builtin functions */
  double sqrt();

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

  *info = 0;

  /*     CHECK THE INPUT PARAMETERS FOR ERRORS. */

  if (*brackt && ((*stp <= min(*stx,*sty) || *stp >= max(*stx,*sty)) || *dx *
                  (*stp - *stx) >=(float)0. || *stpmax < *stpmin)) {
    return 0;
  }

  /*     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. */

  sgnd = *dp *(*dx / abs(*dx));

  /*     FIRST CASE. A HIGHER FUNCTION VALUE. */
  /*     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER */
  /*     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, */
  /*     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. */

  if (*fp > *fx) {
    *info = 1;
    bound = TRUE_;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                                                                          *dp);
    s = max(d__1,d__2);
    /* Computing 2nd power */
    d__1 = theta / s;
    gamma = s * sqrt(d__1 * d__1 - *dx / s *(*dp / s));
    if (*stp < *stx) {
      gamma = -gamma;
    }
    p = gamma - *dx + theta;
    q = gamma - *dx + gamma + *dp;
    r__ = p / q;
    stpc = *stx + r__ *(*stp - *stx);
    stpq = *stx + *dx /((*fx - *fp) / (*stp - *stx) + *dx) / 2 * (*stp -
                                                                  *stx);
    if ((d__1 = stpc - *stx, abs(d__1)) < (d__2 = stpq - *stx, abs(d__2)))
      {
        stpf = stpc;
      } else {
        stpf = stpc + (stpq - stpc) / 2;
      }
    *brackt = TRUE_;

    /*     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF */
    /*     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC */
    /*     STEP IS CLOSER TO STX THAN THE QUADRATIC(SECANT) STEP, */
    /*     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. */

  } else if (sgnd < (float)0.) {
    *info = 2;
    bound = FALSE_;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                                                                          *dp);
    s = max(d__1,d__2);
    /* Computing 2nd power */
    d__1 = theta / s;
    gamma = s * sqrt(d__1 * d__1 - *dx / s *(*dp / s));
    if (*stp > *stx) {
      gamma = -gamma;
    }
    p = gamma - *dp + theta;
    q = gamma - *dp + gamma + *dx;
    r__ = p / q;
    stpc = *stp + r__ *(*stx - *stp);
    stpq = *stp + *dp /(*dp - *dx) * (*stx - *stp);
    if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(d__2)))
      {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    *brackt = TRUE_;

    /*     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
    /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. */
    /*     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY */
    /*     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC */
    /*     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE */
    /*     EITHER STPMIN OR STPMAX. THE QUADRATIC(SECANT) STEP IS ALSO */
    /*     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP */
    /*     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. */

  } else if (abs(*dp) < abs(*dx)) {
    *info = 3;
    bound = TRUE_;
    theta =(*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
    /* Computing MAX */
    d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
                                                                          *dp);
    s = max(d__1,d__2);

    /*        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND */
    /*        TO INFINITY IN THE DIRECTION OF THE STEP. */

    /* Computing MAX */
    /* Computing 2nd power */
    d__3 = theta / s;
    d__1 = 0., d__2 = d__3 * d__3 - *dx / s *(*dp / s);
    gamma = s * sqrt((max(d__1,d__2)));
    if (*stp > *stx) {
      gamma = -gamma;
    }
    p = gamma - *dp + theta;
    q = gamma + (*dx - *dp) + gamma;
    r__ = p / q;
    if (r__ < (float)0. && gamma != (float)0.) {
      stpc = *stp + r__ *(*stx - *stp);
    } else if (*stp > *stx) {
      stpc = *stpmax;
    } else {
      stpc = *stpmin;
    }
    stpq = *stp + *dp /(*dp - *dx) * (*stx - *stp);
    if (*brackt) {
      if ((d__1 = *stp - stpc, abs(d__1)) < (d__2 = *stp - stpq, abs(
                                                                     d__2))) {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    } else {
      if ((d__1 = *stp - stpc, abs(d__1)) > (d__2 = *stp - stpq, abs(
                                                                     d__2))) {
        stpf = stpc;
      } else {
        stpf = stpq;
      }
    }

    /*     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
    /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES */
    /*     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP */
    /*     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. */

  } else {
    *info = 4;
    bound = FALSE_;
    if (*brackt) {
      theta =(*fp - *fy) * 3 / (*sty - *stp) + *dy + *dp;
      /* Computing MAX */
      d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 =
        abs(*dp);
      s = max(d__1,d__2);
      /* Computing 2nd power */
      d__1 = theta / s;
      gamma = s * sqrt(d__1 * d__1 - *dy / s *(*dp / s));
      if (*stp > *sty) {
        gamma = -gamma;
      }
      p = gamma - *dp + theta;
      q = gamma - *dp + gamma + *dy;
      r__ = p / q;
      stpc = *stp + r__ *(*sty - *stp);
      stpf = stpc;
    } else if (*stp > *stx) {
      stpf = *stpmax;
    } else {
      stpf = *stpmin;
    }
  }

  /*     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT */
  /*     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. */

  if (*fp > *fx) {
    *sty = *stp;
    *fy = *fp;
    *dy = *dp;
  } else {
    if (sgnd < (float)0.) {
      *sty = *stx;
      *fy = *fx;
      *dy = *dx;
    }
    *stx = *stp;
    *fx = *fp;
    *dx = *dp;
  }

  /*     COMPUTE THE NEW STEP AND SAFEGUARD IT. */

  stpf = min(*stpmax,stpf);
  stpf = max(*stpmin,stpf);
  *stp = stpf;
  if (*brackt && bound) {
    if (*sty > *stx) {
      /* Computing MIN */
      d__1 = *stx + (*sty - *stx) * (float).66;
      *stp = min(d__1,*stp);
    } else {
      /* Computing MAX */
      d__1 = *stx + (*sty - *stx) * (float).66;
      *stp = max(d__1,*stp);
    }
  }
  return 0;

  /*     LAST LINE OF SUBROUTINE MCSTEP. */

} /* mcstep_ */