void
vmmin(int n0, double *b, double *Fmin, optimfn fminfn, optimgr fmingr,
      int maxit, int trace, int *mask,
      double abstol, double reltol, int nREPORT, void *ex,
      int *fncount, int *grcount, int *fail)
{
  Rboolean accpoint, enough;
  //double *g, *t, *X, *c, **B;
  double *g, *t, *X, *c;
  int   count, funcount, gradcount;
  double f, gradproj;
  int   i, j, ilast, iter = 0;
  double s, steplength;
  double D1, D2;
  int   n, *l;
  int INCX=1, INCY=1;
  double alphaOne=1.0, betaZero=0.0;
  double alphaMOne=-1.0;
  double *B0, *B1, *B2, *B3;
  char *transa = "N", *transb = "T";
  char *UPLO = "U";
  int nn;
  double D1inv;
  
  //Rprintf("Adjusted.\n", *Fmin);
  if (maxit <= 0) {
    *fail = 0;
    *Fmin = fminfn(n0, b, ex);
    *fncount = *grcount = 0;
    return;
  }
  
  if (nREPORT <= 0)
    error(_("REPORT must be > 0 (method = \"BFGS\")"));
  l = (int *) R_alloc(n0, sizeof(int));
  n = 0;
  for (i = 0; i < n0; i++) if (mask[i]) l[n++] = i;
  g = vect(n0);
  t = vect(n);
  X = vect(n);
  c = vect(n);
  //B = Lmatrix(n);
  B0 = vect(n*n);
  //B1 = vect(n*n);
  //B2 = vect(n*n);
  //B3 = vect(n*n);
  f = fminfn(n0, b, ex);
  if (!R_FINITE(f))
    error(_("initial value in 'vmmin' is not finite"));
  if (trace) Rprintf("initial  value %f \n", f);
  *Fmin = f;
  funcount = gradcount = 1;
  fmingr(n0, b, g, ex);
  iter++;
  ilast = gradcount;
  
  do {
    if (ilast == gradcount) {
      // *** Ｂマトリックスを初期化
      // *** 単位行列に
      for (i = 0; i < n; i++) {
        //for (j = 0; j < i; j++) B[i][j] = 0.0;
        //B[i][i] = 1.0;
        for (j = 0; j < n; j++) B0[j+n*i] = 0.0;
        B0[i+n*i] = 1.0;
      }
    }
    for (i = 0; i < n; i++) {
      X[i] = b[l[i]]; /* 元のｂを保存 */
      c[i] = g[l[i]];
    }
    //for (i = 0; i < n; i++) {
    //  F77_CALL(dcopy)(&n, B[i], &INCX, &B0[i*n], &INCY);
    //}
    //for (i = 0; i < n; i++) {
    //  for (j = 0; j <= i; j++){
    //    B0[j+n*i] = B[i][j];
    //    B0[i+n*j] = B[i][j];
    //  }
    //}
    // *** -B0 %*% c -> t
    // *** ｔは降下方向ベクトル
    F77_CALL(dsymv)(UPLO, &n, &alphaMOne, B0, &n, c, &INCX, &betaZero, t, &INCY);
    // *** -c %*% t = -c %*% B0 %*% c -> gradproj
    // gradprojが負なら降下方向
    gradproj = F77_CALL(ddot)(&n, t, &INCX, c, &INCY);
    //Rprintf("gradproj : %f\n", gradproj);
    //gradproj = 0.0;
    //for (i = 0; i < n; i++) {
    //  s = 0.0;
    //  //for (j = 0; j <= i; j++) s -= B[i][j] * g[l[j]];
    //  //for (j = i + 1; j < n; j++) s -= B[j][i] * g[l[j]];
    //  for (j = 0; j <= i; j++) s -= B0[j+n*i] * g[l[j]];
    //  for (j = i + 1; j < n; j++) s -= B0[i+n*j] * g[l[j]];
    //  t[i] = s;
    //  gradproj += s * g[l[i]];
    //}
    //Rprintf("gradproj : %f\n", gradproj);
    
    if (gradproj < 0.0) {	/* search direction is downhill */
      steplength = 1.0;
      accpoint = FALSE;
      do {
        count = 0;
        for (i = 0; i < n; i++) {
          b[l[i]] = X[i] + steplength * t[i];
          if (reltest + X[i] == reltest + b[l[i]]) /* no change */
            count++;
        }
        if (count < n) {
          f = fminfn(n0, b, ex);
          funcount++;
          accpoint = R_FINITE(f) &&
            (f <= *Fmin + gradproj * steplength * acctol);
          if (!accpoint) {
            steplength *= stepredn;
          }
        }
      } while (!(count == n || accpoint));
      enough = (f > abstol) &&
        fabs(f - *Fmin) > reltol * (fabs(*Fmin) + reltol);
      /* stop if value if small or if relative change is low */
      if (!enough) {
        count = n;
        *Fmin = f;
      }
      if (count < n) {/* making progress */
        *Fmin = f;
        fmingr(n0, b, g, ex);
        gradcount++;
        iter++;
        //D1 = 0.0;
        for (i = 0; i < n; i++) {
          //t[i] = steplength * t[i];
          c[i] = g[l[i]] - c[i];
          //D1 += t[i] * c[i];
        }
        F77_CALL(dscal)(&n, &steplength, t, &INCX);
        D1 = F77_CALL(ddot)(&n, t, &INCX, c, &INCY);
        //Rprintf("D1 : %f\n", D1);
        if (D1 > 0) {
          F77_CALL(dsymv)(UPLO, &n, &alphaOne, B0, &n, c, &INCX, &betaZero, X, &INCY);
          D2 = F77_CALL(ddot)(&n, X, &INCX, c, &INCY);
          //Rprintf("D2 : %f\n", D2);
          //D2 = 0.0;
          //for (i = 0; i < n; i++) {
          //  s = 0.0;
          //  for (j = 0; j <= i; j++)
          //    //s += B[i][j] * c[j];
          //    s += B0[j+n*i] * c[j];
          //  for (j = i + 1; j < n; j++)
          //    //s += B[j][i] * c[j];
          //    s += B0[i+n*j] * c[j];
          //  X[i] = s;
          //  D2 += s * c[i];
          //}
          //Rprintf("D2 : %f\n", D2);
          D2 = 1.0 + D2 / D1;
          //for (i = 0; i < n; i++) {
          //  for (j = 0; j <= i; j++)
          //    B[i][j] += (D2 * t[i] * t[j]
          //                - X[i] * t[j] - t[i] * X[j]) / D1;
          //}
          //F77_CALL(dgemm)(transa, transb, &n, &n, &INCX, &D1inv,        t, &n, t, &n, &betaZero, B1, &n);
          //F77_CALL(dgemm)(transa, transb, &n, &n, &INCX, &alphaMOne, X, &n, t, &n, &alphaOne, B1, &n);
          //F77_CALL(dgemm)(transa, transb, &n, &n, &INCX, &alphaMOne, t, &n, X, &n, &alphaOne, B1, &n);
          D1inv = D2/D1;
          F77_CALL(dsyr)(UPLO, &n, &D1inv, t, &INCX, B0, &n);
          D1inv = -1/D1;
          //F77_CALL(dger)(&n, &n, &D1inv, X, &INCX, t, &INCY, B0, &n);
          //F77_CALL(dger)(&n, &n, &D1inv, t, &INCX, X, &INCY, B0, &n);
          F77_CALL(dsyr2)(UPLO, &n, &D1inv, t, &INCX, X, &INCY, B0, &n);
          //nn = n * n;
          //F77_CALL(daxpy)(&nn, &alphaMOne, B2, &INCX, B1, &INCY);
          //F77_CALL(daxpy)(&nn, &alphaMOne, B3, &INCX, B1, &INCY);
          //F77_CALL(dcopy)(&nn, B1, &INCX, B0, &INCY);
          //F77_CALL(dscal)(&nn, &D1, B0, &INCX);
          //D1inv = 1/D1;
          //F77_CALL(daxpy)(&nn, &D1inv, B1, &INCX, B0, &INCY);
          //for (i = 0; i < nn; i++) {
          //  //B0[j+n*i] += B1[j+n*i] / D1;
          //  B0[i] += B1[i] / D1;
          //}
        } else {	/* D1 < 0 */
          ilast = gradcount;
        }
      } else {	/* no progress */
        if (ilast < gradcount) {
          count = 0;
          ilast = gradcount;
        }
      }
    } else {		/* uphill search */
      count = 0;
      if (ilast == gradcount) count = n;
      else ilast = gradcount;
      /* Resets unless has just been reset */
    }
    if (trace && (iter % nREPORT == 0))
      Rprintf("iter%4d value %f\n", iter, f);
    if (iter >= maxit) break;
    if (gradcount - ilast > 2 * n)
      ilast = gradcount;	/* periodic restart */
  } while (count != n || ilast != gradcount);
  if (trace) {
    Rprintf("final  value %f \n", *Fmin);
    if (iter < maxit) Rprintf("converged\n");
    else Rprintf("stopped after %i iterations\n", iter);
  }
  *fail = (iter < maxit) ? 0 : 1;
  *fncount = funcount;
  *grcount = gradcount;
}