gt:       fplus = minclass.apply(x);
context: /*
      java copyright claim:
  
      This software is based on the public domain UNCMIN routines.
      It was translated from FORTRAN to Java by a US government employee
      on official time.  Thus this software is also in the public domain.
  
  
      The translator's mail address is:
  
      Steve Verrill
      USDA Forest Products Laboratory
      1 Gifford Pinchot Drive
      Madison, Wisconsin
      53705
  
  
      The translator's e-mail address is:
  
      steve@www1.fpl.fs.fed.us
  
  
  ***********************************************************************
  
  DISCLAIMER OF WARRANTIES:
  
  THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND.
  THE TRANSLATOR DOES NOT WARRANT, GUARANTEE OR MAKE ANY REPRESENTATIONS
  REGARDING THE SOFTWARE OR DOCUMENTATION IN TERMS OF THEIR CORRECTNESS,
  RELIABILITY, CURRENCY, OR OTHERWISE. THE ENTIRE RISK AS TO
  THE RESULTS AND PERFORMANCE OF THE SOFTWARE IS ASSUMED BY YOU.
  IN NO CASE WILL ANY PARTY INVOLVED WITH THE CREATION OR DISTRIBUTION
  OF THE SOFTWARE BE LIABLE FOR ANY DAMAGE THAT MAY RESULT FROM THE USE
  OF THIS SOFTWARE.
  
  Sorry about that.
  
  ***********************************************************************
  
  
  History:
  
  Date        Translator        Changes
  
  4/14/98     Steve Verrill     Translated
  
  */
  
  
  package org.renjin.stats.internals.optimize;
  
  
  /**
   * <p/>
   * This class contains Java translations of the UNCMIN unconstrained
   * optimization routines.  See R.B. Schnabel, J.E. Koontz, and B.E.
   * Weiss, <i>A Modular System
   * of Algorithms for Unconstrained Minimization</i>, Report CU-CS-240-82,
   * Comp. Sci. Dept., University of Colorado at Boulder, 1982.
   * <p/>
   * <p/>
   * <b>IMPORTANT:</b>  The "_f77" suffixes indicate that these routines use
   * FORTRAN style indexing.  For example, you will see
   * <pre>
   *   for (i = 1; i <= n; i++)
   * </pre>
   * rather than
   * <pre>
   *   for (i = 0; i < n; i++)
   * </pre>
   * To use the "_f77" routines you will have to declare your vectors
   * and matrices to be one element larger (e.g., v[101] rather than
   * v[100], and a[101][101] rather than a[100][100]), and you will have
   * to fill elements 1 through n rather than elements 0 through n - 1.
   * Versions of these programs that use C/Java style indexing will
   * eventually be available.  They will end with the suffix "_j".
   * <p/>
   * <p/>
   * This class was translated by a statistician from a FORTRAN version of
   * UNCMIN.  It is NOT an official translation.  It wastes
   * memory by failing to use the first elements of vectors.  When
   * public domain Java optimization routines become available
   * from the people who produced UNCMIN, then <b>THE CODE PRODUCED
   * BY THE NUMERICAL ANALYSTS SHOULD BE USED</b>.
   * <p/>
   * <p/>
   * Meanwhile, if you have suggestions for improving this
   * code, please contact Steve Verrill at steve@www1.fpl.fs.fed.us.
   *
   * @author (translator)Steve Verrill
   * @version .5 --- April 14, 1998
   */
  class Uncmin {
  
    public enum Method {
      LINE_SEARCH,
      DOUBLE_DOGLEG,
      MORE_HEBDON
    }
  
    /**
     * Minimizes a smooth nonlinear function of n variables.
     * A method that computes the function value at any point
     * must be supplied.  (See Uncmin_methods.java and UncminTest.java.)
     * Derivative values are not required.
     * The optif9 method provides complete user access to the UNCMIN
     * minimization routines.  The user has full control over options.
     * For details, see the Schnabel et al reference and the comments in the code.
     * <p/>
     * Translated by Steve Verrill, August 4, 1998.
     *
     * @param n        The number of arguments of the function to minimize
     * @param x        The initial estimate of the minimum point
     * @param minclass A class that implements the Uncmin_methods
     *                 interface (see the definition in
     *                 Uncmin_methods.java).  See UncminTest_f77.java for an
     *                 example of such a class.  The class must define:
     *                 1.) a method, f_to_minimize, to minimize.
     *                 f_to_minimize must have the form
     *                 <p/>
     *                 public static double f_to_minimize(double x[])
     *                 <p/>
     *                 where x is the vector of arguments to the function
     *                 and the return value is the value of the function
     *                 evaluated at x.
     *                 2.) a method, gradient, that has the form
     *                 <p/>
     *                 public static void gradient(double x[],
     *                 double g[])
     *                 <p/>
     *                 where g is the gradient of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the gradient.
     *                 3.) a method, hessian, that has the form
     *                 <p/>
     *                 public static void hessian(double x[],
     *                 double h[][])
     *                 where h is the Hessian of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the Hessian.  If the user wants Uncmin
     *                 to check the Hessian, then the hessian method
     *                 should only fill the lower triangle (and diagonal)
     *                 of h.
     * @param typsiz   Typical size for each component of x
     * @param fscale   Estimate of the scale of the objective function
     * @param method   Algorithm to use to solve the minimization problem
     *                 = 1  line search
     *                 = 2  double dogleg
     *                 = 3  More-Hebdon
     * @param iexp     = 1 if the optimization function f_to_minimize
     *                 is expensive to evaluate, = 0 otherwise.  If iexp = 1,
     *                 then the Hessian will be evaluated by secant update
     *                 rather than analytically or by finite differences.
     * @param msg      Message to inhibit certain automatic checks and output
     * @param ndigit   Number of good digits in the minimization function
     * @param itnlim   Maximum number of allowable iterations
     * @param iagflg   = 0 if an analytic gradient is not supplied
     * @param iahflg   = 0 if an analytic Hessian is not supplied
     * @param dlt      Trust region radius
     * @param gradtl   Tolerance at which the gradient is considered close enough
     *                 to zero to terminate the algorithm
     * @param stepmx   Maximum allowable step size
     * @param steptl   Relative step size at which successive iterates
     *                 are considered close enough to terminate the algorithm
     * @param xpls     The final estimate of the minimum point
     * @param fpls     The value of f_to_minimize at xpls
     * @param gpls     The gradient at the local minimum xpls
     * @param itrmcd   Termination code
     *                 ITRMCD =  0:  Optimal solution found
     *                 ITRMCD =  1:  Terminated with gradient small,
     *                 X is probably optimal
     *                 ITRMCD =  2:  Terminated with stepsize small,
     *                 X is probably optimal
     *                 ITRMCD =  3:  Lower point cannot be found,
     *                 X is probably optimal
     *                 ITRMCD =  4:  Iteration limit (150) exceeded
     *                 ITRMCD =  5:  Too many large steps,
     *                 function may be unbounded
     * @param a        Workspace for the Hessian (or its estimate)
     *                 and its Cholesky decomposition
     * @param udiag    Workspace for the diagonal of the Hessian
     */
  
    public static void optif9_f77(int n, double x[], UncminFunction minclass,
                                  double typsiz[], double fscale, Method method,
                                  boolean iexp, int msg[], int ndigit[], int itnlim[],
                                  int iagflg[], int iahflg[], double dlt[],
                                  double gradtl[], double stepmx[], double steptl[],
                                  double xpls[], double fpls[], double gpls[],
                                  int itrmcd[], double a[][], double udiag[]) {
  
  
  
  // CHECK FUNCTION SCALE
  
      if (fscale == 0) {
        fscale = 1.0;
      }
      if (fscale < 0.0) {
        fscale = -fscale;
      }
  
  
      double g[] = new double[n + 1];
      double p[] = new double[n + 1];
      double sx[] = new double[n + 1];
      double wrk0[] = new double[n + 1];
      double wrk1[] = new double[n + 1];
      double wrk2[] = new double[n + 1];
      double wrk3[] = new double[n + 1];
  
  // MINIMIZE FUNCTION
  
      optdrv_f77(n, x, minclass, typsiz, fscale, method,
          iexp, msg, ndigit, itnlim, iagflg, iahflg,
          dlt, gradtl, stepmx, steptl, xpls,
          fpls, gpls, itrmcd, a, udiag, g, p, sx,
          wrk0, wrk1, wrk2, wrk3);
  
      if (itrmcd[1] == 1) {
  
  //         System.out.print("\nUncmin WARNING --- itrmcd = 1, " +
  //                          "probably converged, gradient small\n\n");
  
      } else if (itrmcd[1] == 2) {
  
  //         System.out.print("\nUncmin WARNING --- itrmcd = 2, " +
  //                          "probably converged, stepsize small\n\n");
  
      } else if (itrmcd[1] == 3) {
  
  //         System.out.print("\nUncmin WARNING --- itrmcd = 3, " +
  //                          "cannot find lower point\n\n");
  
      } else if (itrmcd[1] == 4) {
  
  //         System.out.print("\nUncmin WARNING --- itrmcd = 4, " +
  //                          "too many iterations\n\n");
  
      } else if (itrmcd[1] == 5) {
  
  //         System.out.print("\nUncmin WARNING --- itrmcd = 5, " +
  //                          "too many large steps, " +
  //                          "possibly unbounded\n\n");
  
      }
    }
  
  
    /**
     * Solves Ax = b where A is an upper triangular
     * matrix.  Note that A is input as a lower triangular matrix and
     * this method takes its transpose implicitly.
     * <p/>
     * Translated by Steve Verrill, April 14, 1998.
     *
     * @param n Dimension of the problem
     * @param a n by n lower triangular matrix (preserved)
     * @param x The solution vector
     * @param b The right-hand side vector
     */
  
    public static void bakslv_f77(int n, double a[][], double x[], double b[]) {
  
      int i, ip1, j;
      double sum;
  
  // SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE)
  
      i = n;
      x[i] = b[i] / a[i][i];
  
      while (i > 1) {
  
        ip1 = i;
        i--;
  
        sum = 0.0;
  
        for (j = ip1; j <= n; j++) {
          sum += a[j][i] * x[j];
        }
  
        x[i] = (b[i] - sum) / a[i][i];
      }
    }
  
  
    /**
     * <p/>
     * The chlhsn_f77 method finds
     * "THE L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION OF THE PERTURBED
     * MODEL HESSIAN MATRIX A+MU*I(WHERE MU\0 AND I IS THE IDENTITY MATRIX)
     * WHICH IS SAFELY POSITIVE DEFINITE.  IF A IS SAFELY POSITIVE DEFINITE
     * UPON ENTRY, THEN MU=0."
     * <p/>
     * Translated by Steve Verrill, April 14, 1998.
     *
     * @param n     Dimension of the problem
     * @param a     On entry: A is the model Hessian (only the lower
     *              triangle and diagonal stored)
     *              On exit: A contains L of the LL+ decomposition of
     *              the perturbed model Hessian in the lower triangle
     *              and diagonal, and contains the Hessian in the upper
     *              triangle and udiag
     * @param epsm  Machine epsilon
     * @param sx    Scaling vector for x
     * @param udiag On exit: Contains the diagonal of the Hessian
     */
  
    public static void chlhsn_f77(int n, double a[][], double epsm,
                                  double sx[], double udiag[]) {
  
  /*C INTERNAL VARIABLES
  C ------------------
  C TOL              TOLERANCE
  C DIAGMN           MINIMUM ELEMENT ON DIAGONAL OF A
  C DIAGMX           MAXIMUM ELEMENT ON DIAGONAL OF A
  C OFFMAX           MAXIMUM OFF-DIAGONAL ELEMENT OF A
  C OFFROW           SUM OF OFF-DIAGONAL ELEMENTS IN A ROW OF A
  C EVMIN            MINIMUM EIGENVALUE OF A
  C EVMAX            MAXIMUM EIGENVALUE OF A
  C
  C DESCRIPTION
  C -----------
  C 1. IF "A" HAS ANY NEGATIVE DIAGONAL ELEMENTS, THEN CHOOSE MU>0
  C SUCH THAT THE DIAGONAL OF A:=A+MU*I IS ALL POSITIVE
  C WITH THE RATIO OF ITS SMALLEST TO LARGEST ELEMENT ON THE
  C ORDER OF SQRT(EPSM).
  C
  C 2. "A" UNDERGOES A PERTURBED CHOLESKY DECOMPOSITION WHICH
  C RESULTS IN AN LL+ DECOMPOSITION OF A+D, WHERE D IS A
  C NON-NEGATIVE DIAGONAL MATRIX WHICH IS IMPLICITLY ADDED TO
  C "A" DURING THE DECOMPOSITION IF "A" IS NOT POSITIVE DEFINITE.
  C "A" IS RETAINED AND NOT CHANGED DURING THIS PROCESS BY
  C COPYING L INTO THE UPPER TRIANGULAR PART OF "A" AND THE
  C DIAGONAL INTO UDIAG.  THEN THE CHOLESKY DECOMPOSITION ROUTINE
  C IS CALLED.  ON RETURN, ADDMAX CONTAINS MAXIMUM ELEMENT OF D.
  C
  C 3. IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2
  C AND RETURN IS MADE TO CALLING PROGRAM.  OTHERWISE,
  C THE MINIMUM NUMBER SDD WHICH MUST BE ADDED TO THE
  C DIAGONAL OF A TO MAKE IT SAFELY STRICTLY DIAGONALLY DOMINANT
  C IS CALCULATED.  SINCE A+ADDMAX*I AND A+SDD*I ARE SAFELY
  C POSITIVE DEFINITE, CHOOSE MU=MIN(ADDMAX,SDD) AND DECOMPOSE
  C A+MU*I TO OBTAIN L.
  C
  
  */
  
      int i, j, im1, jm1;
  
      /**
       * tolerance
       */
      double tol;
  
      double diagmx, diagmn, posmax, amu, offmax,
          evmin, evmax, offrow, sdd;
  
      double addmax[] = new double[2];
  
  // SCALE HESSIAN
  // PRE- AND POST- MULTIPLY "A" BY INV(SX)
  
      for (j = 1; j <= n; j++) {
  
        for (i = j; i <= n; i++) {
  
          a[i][j] /= (sx[i] * sx[j]);
  
        }
  
      }
  
  // STEP1
  // -----
  // NOTE:  IF A DIFFERENT TOLERANCE IS DESIRED THROUGHOUT THIS
  // ALGORITHM, CHANGE TOLERANCE HERE:
  
      tol = Math.sqrt(epsm);
  
      diagmx = a[1][1];
      diagmn = a[1][1];
  
      for (i = 2; i <= n; i++) {
  
        if (a[i][i] < diagmn) diagmn = a[i][i];
        if (a[i][i] > diagmx) diagmx = a[i][i];
  
      }
  
      posmax = Math.max(diagmx, 0.0);
  
  // DIAGMN .LE. 0
  
      if (diagmn <= posmax * tol) {
  
        amu = tol * (posmax - diagmn) - diagmn;
  
        if (amu == 0.0) {
  
  // FIND LARGEST OFF-DIAGONAL ELEMENT OF A
  
          offmax = 0.0;
  
          for (i = 2; i <= n; i++) {
  
            im1 = i - 1;
  
            for (j = 1; j <= im1; j++) {
  
              if (Math.abs(a[i][j]) > offmax) offmax =
                  Math.abs(a[i][j]);
  
            }
  
          }
  
          amu = offmax;
  
          if (amu == 0.0) {
  
            amu = 1.0;
  
          } else {
  
            amu *= 1.0 + tol;
  
          }
  
        }
  
  // A = A + MU*I
  
        for (i = 1; i <= n; i++) {
  
          a[i][i] += amu;
  
        }
  
        diagmx += amu;
  
      }
  
  // STEP2
  // -----
  // COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
  // AND DIAGONAL OF "A" TO UDIAG
  
  
      for (j = 1; j <= n; j++) {
  
        udiag[j] = a[j][j];
  
        for (i = j + 1; i <= n; i++) {
  
          a[j][i] = a[i][j];
  
        }
  
      }
  
      choldc_f77(n, a, diagmx, tol, addmax);
  
  // STEP3
  // -----
  // IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2,
  // THE LL+ DECOMPOSITION HAS BEEN DONE, AND WE RETURN.
  // OTHERWISE, ADDMAX>0.  PERTURB "A" SO THAT IT IS SAFELY
  // DIAGONALLY DOMINANT AND FIND LL+ DECOMPOSITION
  
      if (addmax[1] > 0.0) {
  
  // RESTORE ORIGINAL "A" (LOWER TRIANGULAR PART AND DIAGONAL)
  
        for (j = 1; j <= n; j++) {
  
          a[j][j] = udiag[j];
  
          for (i = j + 1; i <= n; i++) {
  
            a[i][j] = a[j][i];
  
          }
  
        }
  
  // FIND SDD SUCH THAT A+SDD*I IS SAFELY POSITIVE DEFINITE
  // NOTE:  EVMIN<0 SINCE A IS NOT POSITIVE DEFINITE;
  
        evmin = 0.0;
        evmax = a[1][1];
  
        for (i = 1; i <= n; i++) {
  
          offrow = 0.0;
          im1 = i - 1;
  
          for (j = 1; j <= im1; j++) {
  
            offrow += Math.abs(a[i][j]);
  
          }
  
          for (j = i + 1; j <= n; j++) {
  
            offrow += Math.abs(a[j][i]);
  
          }
  
          evmin = Math.min(evmin, a[i][i] - offrow);
          evmax = Math.max(evmax, a[i][i] + offrow);
  
        }
  
        sdd = tol * (evmax - evmin) - evmin;
  
  // PERTURB "A" AND DECOMPOSE AGAIN
  
        amu = Math.min(sdd, addmax[1]);
  
        for (i = 1; i <= n; i++) {
  
          a[i][i] += amu;
          udiag[i] = a[i][i];
  
        }
  
  // "A" NOW GUARANTEED SAFELY POSITIVE DEFINITE
  
        choldc_f77(n, a, 0.0, tol, addmax);
  
      }
  
  // UNSCALE HESSIAN AND CHOLESKY DECOMPOSITION MATRIX
  
      for (j = 1; j <= n; j++) {
  
        for (i = j; i <= n; i++) {
  
          a[i][j] *= sx[i];
  
        }
  
        jm1 = j - 1;
  
        for (i = 1; i <= jm1; i++) {
  
          a[i][j] *= sx[i] * sx[j];
  
        }
  
        udiag[j] *= sx[j] * sx[j];
  
      }
    }
  
  
    /**
     * <p/>
     * The choldc_f77 method finds
     * "THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION
     * OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO A IF
     * NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE."
     * <p/>
     * Translated by Steve Verrill, April 15, 1998.
     *
     * @param n      Dimension of the problem
     * @param a      On entry: matrix for which to find the perturbed
     *               Cholesky decomposition
     *               On exit: contains L of the LL+ decomposition
     *               in lower triangle
     * @param diagmx Maximum diagonal element of "A"
     * @param tol    Tolerance
     * @param addmax Maximum amount implicitly added to diagonal
     *               of "A" in forming the Cholesky decomposition
     *               of A+D
     */
  
    public static void choldc_f77(int n, double a[][], double diagmx,
                                  double tol, double addmax[]) {
  
  /*
  
  Here is a copy of the choldc FORTRAN documentation:
  
        SUBROUTINE CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)
  
  C
  C PURPOSE
  C -------
  C FIND THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION
  C OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO A IF
  C NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE.
  C
  C PARAMETERS
  C ----------
  C NR           --> ROW DIMENSION OF MATRIX
  C N            --> DIMENSION OF PROBLEM
  C A(N,N)      <--> ON ENTRY: MATRIX FOR WHICH TO FIND PERTURBED
  C                       CHOLESKY DECOMPOSITION
  C                  ON EXIT:  CONTAINS L OF LL+ DECOMPOSITION
  C                  IN LOWER TRIANGULAR PART AND DIAGONAL OF "A"
  C DIAGMX       --> MAXIMUM DIAGONAL ELEMENT OF "A"
  C TOL          --> TOLERANCE
  C ADDMAX      <--  MAXIMUM AMOUNT IMPLICITLY ADDED TO DIAGONAL OF "A"
  C                  IN FORMING THE CHOLESKY DECOMPOSITION OF A+D
  C INTERNAL VARIABLES
  C ------------------
  C AMINL    SMALLEST ELEMENT ALLOWED ON DIAGONAL OF L
  C AMNLSQ   =AMINL**2
  C OFFMAX   MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN OF A
  C
  C
  C DESCRIPTION
  C -----------
  C THE NORMAL CHOLESKY DECOMPOSITION IS PERFORMED.  HOWEVER, IF AT ANY
  C POINT THE ALGORITHM WOULD ATTEMPT TO SET L(I,I)=SQRT(TEMP)
  C WITH TEMP < TOL*DIAGMX, THEN L(I,I) IS SET TO SQRT(TOL*DIAGMX)
  C INSTEAD.  THIS IS EQUIVALENT TO ADDING TOL*DIAGMX-TEMP TO A(I,I)
  C
  C
  
  */
  
      int i, j, jm1, jp1, k;
      double aminl, amnlsq, offmax, sum, temp;
  
      addmax[1] = 0.0;
      aminl = Math.sqrt(diagmx * tol);
      amnlsq = aminl * aminl;
  
  // FORM COLUMN J OF L
  
      for (j = 1; j <= n; j++) {
  
  // FIND DIAGONAL ELEMENTS OF L
  
        sum = 0.0;
        jm1 = j - 1;
        jp1 = j + 1;
  
        for (k = 1; k <= jm1; k++) {
  
          sum += a[j][k] * a[j][k];
  
        }
  
        temp = a[j][j] - sum;
  
        if (temp >= amnlsq) {
          a[j][j] = Math.sqrt(temp);
        } else {
  
  // FIND MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN
  
          offmax = 0.0;
  
          for (i = jp1; i <= n; i++) {
            if (Math.abs(a[i][j]) > offmax) {
              offmax = Math.abs(a[i][j]);
            }
          }
  
          if (offmax <= amnlsq) offmax = amnlsq;
  
  // ADD TO DIAGONAL ELEMENT  TO ALLOW CHOLESKY DECOMPOSITION TO CONTINUE
  
          a[j][j] = Math.sqrt(offmax);
          addmax[1] = Math.max(addmax[1], offmax - temp);
  
        }
  
  // FIND I,J ELEMENT OF LOWER TRIANGULAR MATRIX
  
        for (i = jp1; i <= n; i++) {
          sum = 0.0;
  
          for (k = 1; k <= jm1; k++) {
            sum += a[i][k] * a[j][k];
          }
  
          a[i][j] = (a[i][j] - sum) / a[j][j];
        }
      }
    }
  
  
    /**
     * <p/>
     * The dogdrv_f77 method finds the next Newton iterate (xpls) by the double dogleg
     * method.  It drives dogstp_f77.
     * <p/>
     * Translated by Steve Verrill, April 15, 1998.
     *
     * @param n        Dimension of the problem
     * @param x        The old iterate
     * @param f        Function value at the old iterate
     * @param g        Gradient or approximation at the old iterate
     * @param a        Cholesky decomposition of Hessian
     *                 in lower triangular part and diagonal
     * @param p        Newton step
     * @param xpls     The new iterate
     * @param fpls     Function value at the new iterate
     * @param minclass A class that implements the Uncmin_methods
     *                 interface (see the definition in
     *                 Uncmin_methods.java).  See UncminTest_f77.java for an
     *                 example of such a class.  The class must define:
     *                 1.) a method, f_to_minimize, to minimize.
     *                 f_to_minimize must have the form
     *                 <p/>
     *                 public static double f_to_minimize(double x[])
     *                 <p/>
     *                 where x is the vector of arguments to the function
     *                 and the return value is the value of the function
     *                 evaluated at x.
     *                 2.) a method, gradient, that has the form
     *                 <p/>
     *                 public static void gradient(double x[],
     *                 double g[])
     *                 <p/>
     *                 where g is the gradient of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the gradient.
     *                 3.) a method, hessian, that has the form
     *                 <p/>
     *                 public static void hessian(double x[],
     *                 double h[][])
     *                 where h is the Hessian of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the Hessian.  If the user wants Uncmin
     *                 to check the Hessian, then the hessian method
     *                 should only fill the lower triangle (and diagonal)
     *                 of h.
     * @param sx       Scaling vector for x
     * @param stepmx   Maximum allowable step size
     * @param steptl   Relative step size at which successive iterates
     *                 are considered close enough to terminate the
     *                 algorithm
     * @param dlt      Trust region radius (value needs to be retained
     *                 between successive calls)
     * @param iretcd   Return code:
     *                 0 --- satisfactory xpls found
     *                 1 --- failed to find satisfactory xpls
     *                 sufficently distinct from x
     * @param mxtake   Boolean flag indicating that a step of maximum
     *                 length was used length
     * @param sc       Workspace (current step)
     * @param wrk1     Workspace (and place holding argument to tregup)
     * @param wrk2     Workspace
     * @param wrk3     Workspace
     */
  
    public static void dogdrv_f77(int n, double x[], double f[], double g[],
                                  double a[][], double p[], double xpls[],
                                  double fpls[], UncminFunction minclass,
                                  double sx[], double stepmx[], double steptl[],
                                  double dlt[], int iretcd[], boolean mxtake[],
                                  double sc[], double wrk1[], double wrk2[],
                                  double wrk3[]) {
  
  /*
  
  Here is a copy of the dogdrv FORTRAN documentation:
  
        SUBROUTINE DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX,
       +     STEPTL,DLT,IRETCD,MXTAKE,SC,WRK1,WRK2,WRK3,IPR)
  
  C
  C PURPOSE
  C -------
  C FIND A NEXT NEWTON ITERATE (XPLS) BY THE DOUBLE DOGLEG METHOD
  C
  C PARAMETERS
  C ----------
  C NR           --> ROW DIMENSION OF MATRIX
  C N            --> DIMENSION OF PROBLEM
  C X(N)         --> OLD ITERATE X[K-1]
  C F            --> FUNCTION VALUE AT OLD ITERATE, F(X)
  C G(N)         --> GRADIENT  AT OLD ITERATE, G(X), OR APPROXIMATE
  C A(N,N)       --> CHOLESKY DECOMPOSITION OF HESSIAN
  C                  IN LOWER TRIANGULAR PART AND DIAGONAL
  C P(N)         --> NEWTON STEP
  C XPLS(N)     <--  NEW ITERATE X[K]
  C FPLS        <--  FUNCTION VALUE AT NEW ITERATE, F(XPLS)
  C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
  C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
  C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
  C STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
  C                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
  C DLT         <--> TRUST REGION RADIUS
  C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
  C IRETCD      <--  RETURN CODE
  C                    =0 SATISFACTORY XPLS FOUND
  C                    =1 FAILED TO FIND SATISFACTORY XPLS SUFFICIENTLY
  C                       DISTINCT FROM X
  C MXTAKE      <--  BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
  C SC(N)        --> WORKSPACE [CURRENT STEP]
  C WRK1(N)      --> WORKSPACE (AND PLACE HOLDING ARGUMENT TO TREGUP)
  C WRK2(N)      --> WORKSPACE
  C WRK3(N)      --> WORKSPACE
  C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
  C
  
  */
  
      int i;
      double tmp, rnwtln;
      double fplsp[] = new double[2];
      double cln[] = new double[2];
      double eta[] = new double[2];
      boolean fstdog[] = new boolean[2];
      boolean nwtake[] = new boolean[2];
  
  
      iretcd[1] = 4;
      fstdog[1] = true;
      tmp = 0.0;
  
      for (i = 1; i <= n; i++) {
  
        tmp += sx[i] * sx[i] * p[i] * p[i];
  
      }
  
      rnwtln = Math.sqrt(tmp);
  
      while (iretcd[1] > 1) {
  
  // FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
  
        dogstp_f77(n, g, a, p, sx, rnwtln, dlt, nwtake, fstdog,
            wrk1, wrk2, cln, eta, sc, stepmx);
  
  // CHECK NEW POINT AND UPDATE TRUST REGION
  
        tregup_f77(n, x, f, g, a, minclass, sc, sx, nwtake, stepmx,
            steptl, dlt, iretcd, wrk3, fplsp, xpls, fpls,
            mxtake, 2, wrk1);
      }
    }
  
  
    /**
     * <p/>
     * The dogstp_f77 method finds the new step by the double dogleg
     * appproach.
     * <p/>
     * Translated by Steve Verrill, April 21, 1998.
     *
     * @param n      DIMENSION OF PROBLEM
     * @param g      GRADIENT AT CURRENT ITERATE, G(X)
     * @param a      CHOLESKY DECOMPOSITION OF HESSIAN IN
     *               LOWER PART AND DIAGONAL
     * @param p      NEWTON STEP
     * @param sx     Scaling vector for x
     * @param rnwtln NEWTON STEP LENGTH
     * @param dlt    TRUST REGION RADIUS
     * @param nwtake BOOLEAN, = true IF NEWTON STEP TAKEN
     * @param fstdog BOOLEAN, = true IF ON FIRST LEG OF DOGLEG
     * @param ssd    WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE
     *               QUADRATIC MODEL IN THE SCALED STEEPEST DESCENT
     *               DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
     * @param v      WORKSPACE  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
     * @param cln    CAUCHY LENGTH
     *               [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
     * @param eta    [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
     * @param sc     CURRENT STEP
     * @param stepmx MAXIMUM ALLOWABLE STEP SIZE
     */
  
  
    public static void dogstp_f77(int n, double g[], double a[][], double p[],
                                  double sx[], double rnwtln, double dlt[],
                                  boolean nwtake[], boolean fstdog[], double ssd[],
                                  double v[], double cln[], double eta[],
                                  double sc[], double stepmx[]) {
  
  /*
  
  Here is a copy of the dogstp FORTRAN documentation:
  
  C
  C PURPOSE
  C -------
  C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
  C
  C PARAMETERS
  C ----------
  C NR           --> ROW DIMENSION OF MATRIX
  C N            --> DIMENSION OF PROBLEM
  C G(N)         --> GRADIENT AT CURRENT ITERATE, G(X)
  C A(N,N)       --> CHOLESKY DECOMPOSITION OF HESSIAN IN
  C                  LOWER PART AND DIAGONAL
  C P(N)         --> NEWTON STEP
  C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
  C RNWTLN       --> NEWTON STEP LENGTH
  C DLT         <--> TRUST REGION RADIUS
  C NWTAKE      <--> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
  C FSTDOG      <--> BOOLEAN, =.TRUE. IF ON FIRST LEG OF DOGLEG
  C SSD(N)      <--> WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE
  C                  QUADRATIC MODEL IN THE SCALED STEEPEST DESCENT
  C                  DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
  C V(N)        <--> WORKSPACE  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
  C CLN         <--> CAUCHY LENGTH
  C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
  C ETA              [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
  C SC(N)       <--  CURRENT STEP
  C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
  C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
  C
  C INTERNAL VARIABLES
  C ------------------
  C CLN              LENGTH OF CAUCHY STEP
  C
  
  */
  
      double alpha, beta, tmp, dot1, dot2, alam;
      int i, j;
  
  // CAN WE TAKE NEWTON STEP
  
      if (rnwtln <= dlt[1]) {
  
        nwtake[1] = true;
  
        for (i = 1; i <= n; i++) {
  
          sc[i] = p[i];
  
        }
  
        dlt[1] = rnwtln;
  
      } else {
  
  // NEWTON STEP TOO LONG
  // CAUCHY STEP IS ON DOUBLE DOGLEG CURVE
  
        nwtake[1] = false;
  
        if (fstdog[1]) {
  
  // CALCULATE DOUBLE DOGLEG CURVE (SSD)
  
          fstdog[1] = false;
          alpha = 0.0;
  
          for (i = 1; i <= n; i++) {
  
            alpha += (g[i] * g[i]) / (sx[i] * sx[i]);
  
          }
  
          beta = 0.0;
  
          for (i = 1; i <= n; i++) {
  
            tmp = 0.0;
  
            for (j = i; j <= n; j++) {
  
              tmp += (a[j][i] * g[j]) / (sx[j] * sx[j]);
  
            }
  
            beta += tmp * tmp;
  
          }
  
          for (i = 1; i <= n; i++) {
  
            ssd[i] = -(alpha / beta) * g[i] / sx[i];
  
          }
  
          cln[1] = alpha * Math.sqrt(alpha) / beta;
  
          eta[1] = .2 +
              (.8 * alpha * alpha) / (-beta * Blas_f77.ddot_f77(n, g, 1, p, 1));
  
          for (i = 1; i <= n; i++) {
  
            v[i] = eta[1] * sx[i] * p[i] - ssd[i];
  
          }
  
          if (dlt[1] == -1.0) dlt[1] = Math.min(cln[1], stepmx[1]);
  
        }
  
        if (eta[1] * rnwtln <= dlt[1]) {
  
  // TAKE PARTIAL STEP IN NEWTON DIRECTION
  
          for (i = 1; i <= n; i++) {
  
            sc[i] = (dlt[1] / rnwtln) * p[i];
  
          }
  
        } else {
  
          if (cln[1] >= dlt[1]) {
  
  // TAKE STEP IN STEEPEST DESCENT DIRECTION
  
            for (i = 1; i <= n; i++) {
  
              sc[i] = (dlt[1] / cln[1]) * ssd[i] / sx[i];
  
            }
  
          } else {
  
  // CALCULATE CONVEX COMBINATION OF SSD AND ETA*P
  // WHICH HAS SCALED LENGTH DLT
  
            dot1 = Blas_f77.ddot_f77(n, v, 1, ssd, 1);
            dot2 = Blas_f77.ddot_f77(n, v, 1, v, 1);
  
            alam = (-dot1 + Math.sqrt((dot1 * dot1) -
                dot2 * (cln[1] * cln[1] - dlt[1] * dlt[1]))) / dot2;
  
            for (i = 1; i <= n; i++) {
  
              sc[i] = (ssd[i] + alam * v[i]) / sx[i];
  
            }
          }
        }
      }
    }
  
  
    /**
     * <p/>
     * The forslv_f77 method solves Ax = b where A is a lower triangular matrix.
     * <p/>
     * Translated by Steve Verrill, April 21, 1998.
     *
     * @param n The dimension of the problem
     * @param a The lower triangular matrix (preserved)
     * @param x The solution vector
     * @param b The right-hand side vector
     */
  
    public static void forslv_f77(int n, double a[][], double x[],
                                  double b[]) {
  
  /*
  
  Here is a copy of the forslv FORTRAN documentation:
  
        SUBROUTINE FORSLV(NR,N,A,X,B)
  
  C
  C PURPOSE
  C -------
  C SOLVE  AX=B  WHERE A IS LOWER TRIANGULAR MATRIX
  C
  C PARAMETERS
  C ----------
  C NR           --> ROW DIMENSION OF MATRIX
  C N            --> DIMENSION OF PROBLEM
  C A(N,N)       --> LOWER TRIANGULAR MATRIX (PRESERVED)
  C X(N)        <--  SOLUTION VECTOR
  C B(N)         --> RIGHT-HAND SIDE VECTOR
  C
  C NOTE
  C ----
  C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE,
  C THEN VECTORS B AND X MAY SHARE THE SAME STORAGE.
  C
  
  */
  
      int i, im1, j;
      double sum;
  
  // SOLVE LX=B. (FORWARD SOLVE)
  
      x[1] = b[1] / a[1][1];
  
      for (i = 2; i <= n; i++) {
  
        sum = 0.0;
        im1 = i - 1;
  
        for (j = 1; j <= im1; j++) {
  
          sum += a[i][j] * x[j];
  
        }
  
        x[i] = (b[i] - sum) / a[i][i];
  
      }
  
      return;
  
    }
  
  
    /**
     * <p/>
     * The fstocd_f77 method finds a central difference approximation to the
     * gradient of the function to be minimized.
     * <p/>
     * Translated by Steve Verrill, April 21, 1998.
     *
     * @param n        The dimension of the problem
     * @param x        The point at which the gradient is to be approximated
     * @param minclass A class that implements the Uncmin_methods
     *                 interface (see the definition in
     *                 Uncmin_methods.java).  See UncminTest_f77.java for an
     *                 example of such a class.  The class must define:
     *                 1.) a method, f_to_minimize, to minimize.
     *                 f_to_minimize must have the form
     *                 <p/>
     *                 public static double f_to_minimize(double x[])
     *                 <p/>
     *                 where x is the vector of arguments to the function
     *                 and the return value is the value of the function
     *                 evaluated at x.
     *                 2.) a method, gradient, that has the form
     *                 <p/>
     *                 public static void gradient(double x[],
     *                 double g[])
     *                 <p/>
     *                 where g is the gradient of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the gradient.
     *                 3.) a method, hessian, that has the form
     *                 <p/>
     *                 public static void hessian(double x[],
     *                 double h[][])
     *                 where h is the Hessian of f evaluated at x.  This
     *                 method will have an empty body if the user
     *                 does not wish to provide an analytic estimate
     *                 of the Hessian.  If the user wants Uncmin
     *                 to check the Hessian, then the hessian method
     *                 should only fill the lower triangle (and diagonal)
     *                 of h.
     * @param sx       Scaling vector for x
     * @param rnoise   Relative noise in the function to be minimized
     * @param g        A central difference approximation to the gradient
     */
  
  
    public static void fstocd_f77(int n, double x[], UncminFunction minclass,
                                  double sx[], double rnoise, double g[]) {
  
  /*
  
  Here is a copy of the fstocd FORTRAN documentation:
  
        SUBROUTINE FSTOCD (N, X, FCN, SX, RNOISE, G)
  
  C PURPOSE
  C -------
  C FIND CENTRAL DIFFERENCE APPROXIMATION G TO THE FIRST DERIVATIVE
  C (GRADIENT) OF THE FUNCTION DEFINED BY FCN AT THE POINT X.
  C
  C PARAMETERS
  C ----------
  C N            --> DIMENSION OF PROBLEM
  C X            --> POINT AT WHICH GRADIENT IS TO BE APPROXIMATED.
  C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION.
  C SX           --> DIAGONAL SCALING MATRIX FOR X.
  C RNOISE       --> RELATIVE NOISE IN FCN [F(X)].
  C G           <--  CENTRAL DIFFERENCE APPROXIMATION TO GRADIENT.
  C
  C
  
  */
  
      double stepi, xtempi, fplus, fminus, xmult;
      int i;
  
  // FIND I-TH  STEPSIZE, EVALUATE TWO NEIGHBORS IN DIRECTION OF I-TH
  // UNIT VECTOR, AND EVALUATE I-TH  COMPONENT OF GRADIENT.
  
      xmult = Math.pow(rnoise, 1.0 / 3.0);
  
      for (i = 1; i <= n; i++) {
  
        stepi = xmult * Math.max(Math.abs(x[i]), 1.0 / sx[i]);
        xtempi = x[i];
  
        x[i] = xtempi + stepi;
  