/*  MeCab -- Yet Another Part-of-Speech and Morphological Analyzer  $Id: lbfgs.c 1528 2006-08-07 02:39:50Z taku $;  lbfgs.c was ported from the FORTRAN code of lbfgs.m to C  using f2c converter  http://www.ece.northwestern.edu/~nocedal/lbfgs.html  Software for Large-scale Unconstrained Optimization  L-BFGS is a limited-memory quasi-Newton code for unconstrained optimization.  The code has been developed at the Optimization Technology Center,  a joint venture of Argonne National Laboratory and Northwestern University.  Authors  Jorge Nocedal  References  - J. Nocedal. Updating Quasi-Newton Matrices with Limited Storage(1980),  Mathematics of Computation 35, pp. 773-782.  - D.C. Liu and J. Nocedal. On the Limited Memory Method for  Large Scale Optimization(1989),  Mathematical Programming B, 45, 3, pp. 503-528.*/#include <math.h>typedef int integer;typedef unsigned int uinteger;typedef char *address;typedef short int shortint;typedef float real;typedef double doublereal;typedef int logical;typedef short int shortlogical;typedef char logical1;typedef char integer1;#define TRUE_ (1)#define FALSE_ (0)#define abs(x) ((x) >= 0 ? (x) : -(x))#define min(a,b) ((a) <= (b) ? (a) : (b))#define max(a,b) ((a) >= (b) ? (a) : (b))/* Common Block Declarations */struct lb3_1_ {  integer mp, lp;  doublereal gtol, stpmin, stpmax;};/* Table of constant values */static struct lb3_1_ lb3_1;static integer c__1 = 1;static doublereal ddot_();static int daxpy_();static int mcsrch_();static int mcstep_();/*     ---------------------------------------------------------------------- *//*     This file contains the LBFGS algorithm and supporting routines *//*     **************** *//*     LBFGS SUBROUTINE *//*     **************** *//* Subroutine */ int lbfgs(n, m, x, f, g, diagco, diag, iprint, eps, xtol, w, iflag)     integer *n, *m;     doublereal *x, *f, *g;     integer *diagco;     doublereal *diag;     integer *iprint;     doublereal *eps, *xtol, *w;     integer *iflag;{  /* Initialized data */  static doublereal one = 1.0;  static doublereal zero = 0.0;  /* System generated locals */  integer i__1;  doublereal d__1;  /* Builtin functions */  double sqrt();  /* Local variables */  static doublereal beta;  static integer inmc;  static integer info, iscn, nfev, iycn, iter;  static doublereal ftol;  static integer nfun, ispt, iypt, i__, bound;  static doublereal gnorm;  static integer point;  static doublereal xnorm;  static integer cp;  static doublereal sq, yr, ys;  static logical finish;  static doublereal yy;  static integer maxfev;  static integer npt;  static doublereal stp, stp1;  lb3_1.mp = 6;  lb3_1.lp = 6;  lb3_1.gtol = .9;  lb3_1.stpmin = 1e-20;  lb3_1.stpmax = 1e20;  /* Parameter adjustments */  --diag;  --g;  --x;  --w;  --iprint;  /* Function Body */  /*     INITIALIZE */  /*     ---------- */  if (*iflag == 0) {    goto L10;  }  switch((int)*iflag) {  case 1:  goto L172;  case 2:  goto L100;  } L10:  iter = 0;  if (*n <= 0 || *m <= 0) {    goto L196;  }  if (lb3_1.gtol <= 1e-4) {    if (lb3_1.lp > 0) {}    lb3_1.gtol = .9;  }  nfun = 1;  point = 0;  finish = FALSE_;  if (*diagco != 0) {    i__1 = *n;    for (i__ = 1; i__ <= i__1; ++i__) {      /* L30: */      if (diag[i__] <= zero) {        goto L195;      }    }  } else {    i__1 = *n;    for (i__ = 1; i__ <= i__1; ++i__) {      /* L40: */      diag[i__] = 1.;    }  }  ispt = *n + (*m << 1);  iypt = ispt + *n * *m;  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    /* L50: */    w[ispt + i__] = -g[i__] * diag[i__];  }  gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));  stp1 = one / gnorm;  /*     PARAMETERS FOR LINE SEARCH ROUTINE */  ftol = 1e-4;  maxfev = 20;  /*    if (iprint[1] >= 0) {        lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &        finish);        } */  /*    -------------------- */  /*     MAIN ITERATION LOOP */  /*    -------------------- */ L80:  ++iter;  info = 0;  bound = iter - 1;  if (iter == 1) {    goto L165;  }  if (iter > *m) {    bound = *m;  }  ys = ddot_(n, &w[iypt + npt + 1], &c__1, &w[ispt + npt + 1], &c__1);  if (*diagco == 0) {    yy = ddot_(n, &w[iypt + npt + 1], &c__1, &w[iypt + npt + 1], &c__1);    i__1 = *n;    for (i__ = 1; i__ <= i__1; ++i__) {      /* L90: */      diag[i__] = ys / yy;    }  } else {    *iflag = 2;    return 0;  } L100:  if (*diagco != 0) {    i__1 = *n;    for (i__ = 1; i__ <= i__1; ++i__) {      /* L110: */      if (diag[i__] <= zero) {        goto L195;      }    }  }  /*     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;  if (point == 0) {    cp = *m;  }  w[*n + cp] = one / ys;  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    /* L112: */    w[i__] = -g[i__];  }  cp = point;  i__1 = bound;  for (i__ = 1; i__ <= i__1; ++i__) {    --cp;    if (cp == -1) {      cp = *m - 1;    }    sq = ddot_(n, &w[ispt + cp * *n + 1], &c__1, &w[1], &c__1);    inmc = *n + *m + cp + 1;    iycn = iypt + cp * *n;    w[inmc] = w[*n + cp + 1] * sq;    d__1 = -w[inmc];    daxpy_(n, &d__1, &w[iycn + 1], &c__1, &w[1], &c__1);    /* L125: */  }  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    /* L130: */    w[i__] = diag[i__] * w[i__];  }  i__1 = bound;  for (i__ = 1; i__ <= i__1; ++i__) {    yr = ddot_(n, &w[iypt + cp * *n + 1], &c__1, &w[1], &c__1);    beta = w[*n + cp + 1] * yr;    inmc = *n + *m + cp + 1;    beta = w[inmc] - beta;    iscn = ispt + cp * *n;    daxpy_(n, &beta, &w[iscn + 1], &c__1, &w[1], &c__1);    ++cp;    if (cp == *m) {      cp = 0;    }    /* L145: */  }  /*     STORE THE NEW SEARCH DIRECTION */  /*     ------------------------------ */  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    /* L160: */    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;  stp = one;  if (iter == 1) {    stp = stp1;  }  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    /* L170: */    w[i__] = g[i__];  } L172:  mcsrch_(n, &x[1], f, &g[1], &w[ispt + point * *n + 1], &stp, &ftol, xtol,          &maxfev, &info, &nfev, &diag[1]);  if (info == -1) {    *iflag = 1;    return 0;  }  if (info != 1) {    goto L190;  }  nfun += nfev;  /*     COMPUTE THE NEW STEP AND GRADIENT CHANGE */  /*     ----------------------------------------- */  npt = point * *n;  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    w[ispt + npt + i__] = stp * w[ispt + npt + i__];    /* L175: */    w[iypt + npt + i__] = g[i__] - w[i__];  }  ++point;  if (point == *m) {    point = 0;  }  /*     TERMINATION TEST */  /*     ---------------- */  gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));  xnorm = sqrt(ddot_(n, &x[1], &c__1, &x[1], &c__1));  xnorm = max(1.,xnorm);  if (gnorm / xnorm <= *eps) {    finish = TRUE_;  }  /*    if (iprint[1] >= 0) {        lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &        finish);        }*/  if (finish) {    *iflag = 0;    return 0;  }  goto L80;  /*     ------------------------------------------------------------ */  /*     END OF MAIN ITERATION LOOP. ERROR EXITS. */  /*     ------------------------------------------------------------ */ L190:  *iflag = -1;  return 0; L195:  *iflag = -2;  return 0; L196:  *iflag = -3;  return 0;} /* lbfgs_ *//*   ---------------------------------------------------------- *//* Subroutine */ static int daxpy_(n, da, dx, incx, dy, incy)     integer *n;     doublereal *da, *dx;     integer *incx;     doublereal *dy;     integer *incy;{  /* System generated locals */  integer i__1;  /* Local variables */  static integer i__, m, ix, iy, mp1;  /*     constant times a vector plus a vector. */  /*     uses unrolled loops for increments equal to one. */  /*     jack dongarra, linpack, 3/11/78. */  /* Parameter adjustments */  --dy;  --dx;  /* Function Body */  if (*n <= 0) {    return 0;  }  if (*da == 0.) {    return 0;  }  if (*incx == 1 && *incy == 1) {    goto L20;  }  /*        code for unequal increments or equal increments */  /*          not equal to 1 */  ix = 1;  iy = 1;  if (*incx < 0) {    ix =(-(*n) + 1) * *incx + 1;  }  if (*incy < 0) {    iy =(-(*n) + 1) * *incy + 1;  }  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {    dy[iy] += *da * dx[ix];    ix += *incx;    iy += *incy;    /* L10: */  }  return 0;  /*        code for both increments equal to 1 */  /*        clean-up loop */ L20:  m = *n % 4;  if (m == 0) {    goto L40;  }  i__1 = m;  for (i__ = 1; i__ <= i__1; ++i__) {    dy[i__] += *da * dx[i__];    /* L30: */  }  if (*n < 4) {    return 0;  } L40:  mp1 = m + 1;  i__1 = *n;  for (i__ = mp1; i__ <= i__1; i__ += 4) {    dy[i__] += *da * dx[i__];    dy[i__ + 1] += *da * dx[i__ + 1];    dy[i__ + 2] += *da * dx[i__ + 2];    dy[i__ + 3] += *da * dx[i__ + 3];    /* L50: */  }  return 0;} /* daxpy_ *//*   ---------------------------------------------------------- */static doublereal ddot_(n, dx, incx, dy, incy)     integer *n;     doublereal *dx;     integer *incx;     doublereal *dy;     integer *incy;{  /* System generated locals */  integer i__1;  doublereal ret_val;  /* Local variables */  static integer i__, m;  static doublereal dtemp;  static integer ix, iy, mp1;  /*     forms the dot product of two vectors. */  /*     uses unrolled loops for increments equal to one. */  /*     jack dongarra, linpack, 3/11/78. */  /* Parameter adjustments */  --dy;  --dx;  /* Function Body */  ret_val = 0.;  dtemp = 0.;  if (*n <= 0) {    return ret_val;  }  if (*incx == 1 && *incy == 1) {    goto L20;  }  /*        code for unequal increments or equal increments */  /*          not equal to 1 */  ix = 1;  iy = 1;  if (*incx < 0) {    ix =(-(*n) + 1) * *incx + 1;  }  if (*incy < 0) {    iy =(-(*n) + 1) * *incy + 1;  }  i__1 = *n;  for (i__ = 1; i__ <= i__1; ++i__) {