docs: public static StatusCode lbfgs(
          double[] x,
          MutableDouble ptr_fx,
          LBFGSCallback cd,
          LBFGSPrm _param
          )     
      {
          final int n = x.length;
          StatusCode ret, ls_ret;
          MutableInt ls = new MutableInt(0);
          int i, j, k, end, bound;        
          double step;
      
          /* Constant parameters and their default values. */
          LBFGSPrm param = (_param != null) ? _param : new LBFGSPrm();
          final int m = param.m;
  
          double[] xp;
          double[] g, gp, pg = null;
          double[] d, w, pf = null;
          IterationData[] lm = null;
          IterationData it = null;
          double ys, yy;
          double xnorm, gnorm, beta;
          double fx = 0.;
          double rate = 0.;
          LineSearchAlgInternal linesearch_choice = LineSearchAlgInternal.morethuente;
  
          /* Check the input parameters for errors. */
          if (n <= 0) {
              return LBFGSERR_INVALID_N;
          }
          if (param.epsilon < 0.) {
              return LBFGSERR_INVALID_EPSILON;
          }
          if (param.past < 0) {
              return LBFGSERR_INVALID_TESTPERIOD;
          }
          if (param.delta < 0.) {
              return LBFGSERR_INVALID_DELTA;
          }
          if (param.min_step < 0.) {
              return LBFGSERR_INVALID_MINSTEP;
          }
          if (param.max_step < param.min_step) {
              return LBFGSERR_INVALID_MAXSTEP;
          }
          if (param.ftol < 0.) {
              return LBFGSERR_INVALID_FTOL;
          }
          if (param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_WOLFE ||
              param.linesearch == LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE) {
              if (param.wolfe <= param.ftol || 1. <= param.wolfe) {
                  return LBFGSERR_INVALID_WOLFE;
              }
          }
          if (param.gtol < 0.) {
              return LBFGSERR_INVALID_GTOL;
          }
          if (param.xtol < 0.) {
              return LBFGSERR_INVALID_XTOL;
          }
          if (param.max_linesearch <= 0) {
              return LBFGSERR_INVALID_MAXLINESEARCH;
          }
          if (param.orthantwise_c < 0.) {
              return LBFGSERR_INVALID_ORTHANTWISE;
          }
          if (param.orthantwise_start < 0 || n < param.orthantwise_start) {
              return LBFGSERR_INVALID_ORTHANTWISE_START;
          }
          if (param.orthantwise_end < 0) {
              param.orthantwise_end = n;
          }
          if (n < param.orthantwise_end) {
              return LBFGSERR_INVALID_ORTHANTWISE_END;
          }
          if (param.orthantwise_c != 0.) {
              switch (param.linesearch) {
              case LBFGS_LINESEARCH_BACKTRACKING_WOLFE:
                  linesearch_choice = LineSearchAlgInternal.backtracking_owlqn;
                  break;
              default:
                  /* Only the backtracking method is available. */
                  return LBFGSERR_INVALID_LINESEARCH;
              }
          } else {
              switch (param.linesearch) {
              case LBFGS_LINESEARCH_MORETHUENTE:
                  linesearch_choice = LineSearchAlgInternal.morethuente;
                  break;
              case LBFGS_LINESEARCH_BACKTRACKING_ARMIJO:
              case LBFGS_LINESEARCH_BACKTRACKING_WOLFE:
              case LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE:
                  linesearch_choice = LineSearchAlgInternal.backtracking;
                  break;
              default:
                  return LBFGSERR_INVALID_LINESEARCH;
              }
          }
  
          /* Allocate working space. */
          xp = new double[n];
          g = new double[n];
          gp = new double[n];
          d = new double[n];
          w = new double[n];
  
          if (param.orthantwise_c != 0.) {
              /* Allocate working space for OW-LQN. */
              pg = new double[n];
          }
  
          /* Allocate limited memory storage. */
          lm = new IterationData[m];
  
          /* Initialize the limited memory. */
          for (i = 0;i < m;++i) {
              lm[i] = new IterationData();
              it = lm[i];
              it.alpha = 0;
              it.ys = 0;
              it.s = new double[n];
              it.y = new double[n];
          }
  
          /* Allocate an array for storing previous values of the objective function. */
          if (0 < param.past) {
              pf = new double[param.past];
          }
  
          /* Evaluate the function value and its gradient. */
          fx = cd.proc_evaluate(x, g, 0);
          if (0. != param.orthantwise_c) {
              /* Compute the L1 norm of the variable and add it to the object value. */
              xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end);
              fx += xnorm * param.orthantwise_c;
              owlqn_pseudo_gradient(
                  pg, x, g, n,
                  param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
                  );
          }
  
          /* Store the initial value of the objective function. */
          if (pf != null) {
              pf[0] = fx;
          }
  
          /*
              Compute the direction;
              we assume the initial hessian matrix H_0 as the identity matrix.
           */
          if (param.orthantwise_c == 0.) {
              vecncpy(d, g, n);
          } else {
              vecncpy(d, pg, n);
          }
  
          /*
             Make sure that the initial variables are not a minimizer.
           */
          xnorm = vec2norm(x, n);
          if (param.orthantwise_c == 0.) {
              gnorm = vec2norm(g, n);
          } else {
              gnorm = vec2norm(pg, n);
          }
          if (xnorm < 1.0) xnorm = 1.0;
          if (gnorm / xnorm <= param.epsilon) {
              ptr_fx.v = fx;
              return LBFGS_ALREADY_MINIMIZED;
          }
  
          /* Compute the initial step:
              step = 1.0 / sqrt(vecdot(d, d, n))
           */
          step = vec2norminv(d, n);
  
          k = 1;
          end = 0;
          for (;;) {
              /* Store the current position and gradient vectors. */
              veccpy(xp, x, n);
              veccpy(gp, g, n);
  
              /* Search for an optimal step. */
              MutableDouble fxRef = new MutableDouble(fx);
              MutableDouble stepRef = new MutableDouble(step);
              if (param.orthantwise_c == 0.) {
                  ls_ret = linesearch(n, x, fxRef, g, d, stepRef, xp, gp, w, cd, param, ls, linesearch_choice);
              } else {
                  ls_ret = linesearch(n, x, fxRef, g, d, stepRef, xp, pg, w, cd, param, ls, linesearch_choice);
                  owlqn_pseudo_gradient(
                      pg, x, g, n,
                      param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
                      );
              }
              fx = fxRef.v;
              step = stepRef.v;
              
              if (ls_ret.ret < 0) {
                  /* Revert to the previous point. */
                  veccpy(x, xp, n);
                  veccpy(g, gp, n);
                  ptr_fx.v = fx;
                  return ls_ret;
              }
  
              /* Compute x and g norms. */
              xnorm = vec2norm(x, n);
              if (param.orthantwise_c == 0.) {
                  gnorm = vec2norm(g, n);
              } else {
                  gnorm = vec2norm(pg, n);
              }
  
              /* Report the progress. */
              ret = cd.proc_progress(x, g, fx, xnorm, gnorm, step, k, ls.v);
              if (ret.ret != 0) {
                  ptr_fx.v = fx;
                  return ret;
              }
  
              /*
                  Convergence test.
                  The criterion is given by the following formula:
                      |g(x)| / \max(1, |x|) < \epsilon
               */
              if (xnorm < 1.0) xnorm = 1.0;
              if (gnorm / xnorm <= param.epsilon) {
                  /* Convergence. */
                  ret = LBFGS_SUCCESS;
                  break;
              }
  
              /*
                  Test for stopping criterion.
                  The criterion is given by the following formula:
                      (f(past_x) - f(x)) / f(x) < \delta
               */
              if (pf != null) {
                  /* We don't test the stopping criterion while k < past. */
                  if (param.past <= k) {
                      /* Compute the relative improvement from the past. */
                      rate = (pf[k % param.past] - fx) / fx;
  
                      /* The stopping criterion. */
                      if (rate < param.delta) {
                          ret = LBFGS_STOP;
                          break;
                      }
                  }
  
                  /* Store the current value of the objective function. */
                  pf[k % param.past] = fx;
              }
  
              if (param.max_iterations != 0 && param.max_iterations < k+1) {
                  /* Maximum number of iterations. */
                  ret = LBFGSERR_MAXIMUMITERATION;
                  break;
              }
  
              /*
                  Update vectors s and y:
                      s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}.
                      y_{k+1} = g_{k+1} - g_{k}.
               */
              it = lm[end];
              vecdiff(it.s, x, xp, n);
              vecdiff(it.y, g, gp, n);
  
              /*
                  Compute scalars ys and yy:
                      ys = y^t \cdot s = 1 / \rho.
                      yy = y^t \cdot y.
                  Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor).
               */
              ys = vecdot(it.y, it.s, n);
              yy = vecdot(it.y, it.y, n);
              it.ys = ys;
  
              /*
                  Recursive formula to compute dir = -(H \cdot g).
                      This is described in page 779 of:
                      Jorge Nocedal.
                      Updating Quasi-Newton Matrices with Limited Storage.
                      Mathematics of Computation, Vol. 35, No. 151,
                      pp. 773--782, 1980.
               */
              bound = (m <= k) ? m : k;
              ++k;
              end = (end + 1) % m;
  
              /* Compute the steepest direction. */
              if (param.orthantwise_c == 0.) {
                  /* Compute the negative of gradients. */
                  vecncpy(d, g, n);
              } else {
                  vecncpy(d, pg, n);
              }
  
              j = end;
              for (i = 0;i < bound;++i) {
                  j = (j + m - 1) % m;    /* if (--j == -1) j = m-1; */
                  it = lm[j];
                  /* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */
                  it.alpha = vecdot(it.s, d, n);
                  it.alpha /= it.ys;
                  /* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */
                  vecadd(d, it.y, -it.alpha, n);
              }
  
              vecscale(d, ys / yy, n);
  
              for (i = 0;i < bound;++i) {
                  it = lm[j];
                  /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */
                  beta = vecdot(it.y, d, n);
                  beta /= it.ys;
                  /* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */
                  vecadd(d, it.s, it.alpha - beta, n);
                  j = (j + 1) % m;        /* if (++j == m) j = 0; */
              }
  
              /*
                  Constrain the search direction for orthant-wise updates.
               */
              if (param.orthantwise_c != 0.) {
                  for (i = param.orthantwise_start;i < param.orthantwise_end;++i) {
                      if (d[i] * pg[i] >= 0) {
                          d[i] = 0;
                      }
                  }
              }
  
              /*
                  Now the search direction d is ready. We try step = 1 first.
               */
              step = 1.0;
          }
  
          /* Return the final value of the objective function. */
          ptr_fx.v = fx;
          
          return ret;
      }