static double linesearch(MultivariateFunction func, double[] xold, double fold, double[] g, double[] p, double[] x, double stpmax) {
        if (stpmax <= 0) {
            throw new IllegalArgumentException("Invalid upper bound of linear search step: " + stpmax);
        }

        // Termination occurs when the relative width of the interval
        // of uncertainty is at most xtol.
        final double xtol = EPSILON;
        // Tolerance for the sufficient decrease condition.
        final double ftol = 1.0E-4;

        int n = xold.length;

        // Scale if attempted step is too big
        double pnorm = norm(p);
        if (pnorm > stpmax) {
            double r = stpmax / pnorm;
            for (int i = 0; i < n; i++) {
                p[i] *= r;
            }
        }

        // Check if s is a descent direction.
        double slope = 0.0;
        for (int i = 0; i < n; i++) {
            slope += g[i] * p[i];
        }

        if (slope >= 0) {
            throw new IllegalArgumentException("Line Search: the search direction is not a descent direction, which may be caused by roundoff problem.");
        }

        // Calculate minimum step.
        double test = 0.0;
        for (int i = 0; i < n; i++) {
            double temp = abs(p[i]) / max(xold[i], 1.0);
            if (temp > test) {
                test = temp;
            }
        }

        double alammin = xtol / test;
        double alam = 1.0;

        double alam2 = 0.0, f2 = 0.0;
        double a, b, disc, rhs1, rhs2, tmpalam;
        while (true) {
            // Evaluate the function and gradient at stp
            // and compute the directional derivative.
            for (int i = 0; i < n; i++) {
                x[i] = xold[i] + alam * p[i];
            }

            double f = func.f(x);

            // Convergence on &Delta; x.
            if (alam < alammin) {
                System.arraycopy(xold, 0, x, 0, n);
                return f;
            } else if (f <= fold + ftol * alam * slope) {
                // Sufficient function decrease.
                return f;
            } else {
                // Backtrack
                if (alam == 1.0) {
                    // First time
                    tmpalam = -slope / (2.0 * (f - fold - slope));
                } else {
                    // Subsequent backtracks.
                    rhs1 = f - fold - alam * slope;
                    rhs2 = f2 - fold - alam2 * slope;
                    a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
                    b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
                    if (a == 0.0) {
                        tmpalam = -slope / (2.0 * b);
                    } else {
                        disc = b * b - 3.0 * a * slope;
                        if (disc < 0.0) {
                            tmpalam = 0.5 * alam;
                        } else if (b <= 0.0) {
                            tmpalam = (-b + sqrt(disc)) / (3.0 * a);
                        } else {
                            tmpalam = -slope / (b + sqrt(disc));
                        }
                    }
                    if (tmpalam > 0.5 * alam) {
                        tmpalam = 0.5 * alam;
                    }
                }
            }
            alam2 = alam;
            f2 = f;
            alam = max(tmpalam, 0.1 * alam);
        }
    }
static double linesearch ( MultivariateFunction func , double [ ] xold , double fold , double [ ] g , double [ ] p , double [ ] x , double stpmax ) { if ( stpmax <= 0 ) { throw new IllegalArgumentException ( "Invalid upper bound of linear search step: " + stpmax ) ; } final double xtol = EPSILON ; final double ftol = 1.0E-4 ; int n = xold . length ; double pnorm = norm ( p ) ; if ( pnorm > stpmax ) { double r = stpmax / pnorm ; for ( int i = 0 ; i < n ; i ++ ) { p [ i ] *= r ; } } double slope = 0.0 ; for ( int i = 0 ; i < n ; i ++ ) { slope += g [ i ] * p [ i ] ; } if ( slope >= 0 ) { throw new IllegalArgumentException ( "Line Search: the search direction is not a descent direction, which may be caused by roundoff problem." ) ; } double test = 0.0 ; for ( int i = 0 ; i < n ; i ++ ) { double temp = abs ( p [ i ] ) / max ( xold [ i ] , 1.0 ) ; if ( temp > test ) { test = temp ; } } double alammin = xtol / test ; double alam = 1.0 ; double alam2 = 0.0 , f2 = 0.0 ; double a , b , disc , rhs1 , rhs2 , tmpalam ; while ( true ) { for ( int i = 0 ; i < n ; i ++ ) { x [ i ] = xold [ i ] + alam * p [ i ] ; } double f = func . f ( x ) ; if ( alam < alammin ) { System . arraycopy ( xold , 0 , x , 0 , n ) ; return f ; } else if ( f <= fold + ftol * alam * slope ) { return f ; } else { if ( alam == 1.0 ) { tmpalam = - slope / ( 2.0 * ( f - fold - slope ) ) ; } else { rhs1 = f - fold - alam * slope ; rhs2 = f2 - fold - alam2 * slope ; a = ( rhs1 / ( alam * alam ) - rhs2 / ( alam2 * alam2 ) ) / ( alam - alam2 ) ; b = ( - alam2 * rhs1 / ( alam * alam ) + alam * rhs2 / ( alam2 * alam2 ) ) / ( alam - alam2 ) ; if ( a == 0.0 ) { tmpalam = - slope / ( 2.0 * b ) ; } else { disc = b * b - 3.0 * a * slope ; if ( disc < 0.0 ) { tmpalam = 0.5 * alam ; } else if ( b <= 0.0 ) { tmpalam = ( - b + sqrt ( disc ) ) / ( 3.0 * a ) ; } else { tmpalam = - slope / ( b + sqrt ( disc ) ) ; } } if ( tmpalam > 0.5 * alam ) { tmpalam = 0.5 * alam ; } } } alam2 = alam ; f2 = f ; alam = max ( tmpalam , 0.1 * alam ) ; } }
Minimize a function along a search direction by find a step which satisfies
a sufficient decrease condition and a curvature condition.
<p>
At each stage this function updates an interval of uncertainty with
endpoints <code>stx</code> and <code>sty</code>. The interval of
uncertainty is initially chosen so that it contains a
minimizer of the modified function
<pre>
f(x+stp*s) - f(x) - ftol*stp*(gradf(x)'s).
</pre>
If a step is obtained for which the modified function
has a nonpositive function value and nonnegative derivative,
then the interval of uncertainty is chosen so that it
contains a minimizer of <code>f(x+stp*s)</code>.
<p>
The algorithm is designed to find a step which satisfies
the sufficient decrease condition
<pre>
f(x+stp*s) &lt;= f(X) + ftol*stp*(gradf(x)'s),
</pre>
and the curvature condition
<pre>
abs(gradf(x+stp*s)'s)) &lt;= gtol*abs(gradf(x)'s).
</pre>
If <code>ftol</code> is less than <code>gtol</code> and if, for example,
the function is bounded below, then there is always a step which
satisfies both conditions. If no step can be found which satisfies both
conditions, then the algorithm usually stops when rounding
errors prevent further progress. In this case <code>stp</code> only
satisfies the sufficient decrease condition.
<p>

@param xold on input this contains the base point for the line search.

@param fold on input this contains the value of the objective function
at <code>x</code>.

@param g on input this contains the gradient of the objective function
at <code>x</code>.

@param p the search direction.

@param x on output, it contains <code>xold + alam*p</code>.

@param stpmax specify upper bound for the step in the line search so that
we do not try to evaluate the function in regions where it is
undefined or subject to overflow.

@return the new function value.