/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:46 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sfmin.h"
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
/* PARAMETER translations */
#define FIVE 5.0e0
#define HALF 0.5e0
#define PT95 0.95e0
#define THREE 3.0e0
#define TWO 2.0e0
#define ZERO 0.0e0
/* end of PARAMETER translations */
void /*FUNCTION*/ sfmin(
float *x,
float *xorf,
long *mode,
float tol)
{
static long int _l0, next;
static float a, asave, b, bsave, c, c3, d, e, fu, fv, fw, fx,
p, q, r, tol1, tol2, toli, u, v, w, xi, xm;
static float eps = ZERO;
static long np = 0;
static long ic = 0;
/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
* ALL RIGHTS RESERVED.
* Based on Government Sponsored Research NAS7-03001.
*>> 2004-11-11 SFMIN Krogh Removed if (..) C3 = statement doing nil.
*>> 2000-12-01 SFMIN Krogh Removed unused parameter ONE.
*>> 1996-03-30 SFMIN Krogh Added external statement.
*>> 1994-10-20 SFMIN Krogh Changes to use M77CON
*>> 1994-04-20 SFMIN CLL Edited to make DP & SP files similar.
*>> 1987-12-09 SFMIN Lawson Initial code.
*
* Finds a local minimum of a function f(x) in the closed interval
* between A and B.
* The function f(x) is defined by code in the calling program.
* The user must initially provide A, B, and a tolerance, TOL.
* The function will not be evaluated at A or B unless the
* minimization search leads to one of these endpoints.
*
* This subroutine uses reverse communication, i.e., it returns to
* the calling program each time it needs to have f() evaluated at a
* new value of x.
* ------------------------------------------------------------------
* Subroutine Arguments
*
* X, XORF, MODE [all are inout]
* On the initial call to this subroutine to solve a
* new problem, the user must set MODE = 0, and must set
* X = A
* XORF = B
* TOL = desired tolerance
* A and B denote endpoints defining a closed
* interval in which a local minimum is to be found.
* Permit A < B or A > B or A = B. If A = B the solution, x = A,
* will be found immediately.
* On any call the user can set MODE negative to start or stop
* detail printing. On such an entry we set NP = -1 - MODE, and
* the normal algorithm will not be executed. On the next NP calls
* with MODE = 0 or 1, a detail line will be printed.
*
* TOL [in] is an absolute tolerance on the uncertainty in the final
* estimate of the minimizing abcissa, X1. Recommend TOL .ge. 0.
* This subr will set TOLI = TOL if TOL > 0., and otherwise sets
* TOLI = 0.. The operational tolerance at any trial abcissa, X,
* will be
* TOL2 = 2 * abs(X) * sqrt(R1MACH(4)) + (2/3) * C3
* where C3 = TOLI/3.E0 + R1MACH(4)**2 where TOLI = max(TOL,0.E0).
*
* On each return this subr will set MODE to a value in the range
* [1:3] to indicate the action needed from the calling program or
* the status on termination.
*
* = 1 means the calling program must evaluate
* f(X), store the value in XORF, and then call this
* subr again.
*
* = 2 means normal termination. X contains the point of the
* minimum and XORF contains f(X).
*
* = 3 as for 2 but the requested accuracy could not be obtained.
*
* = 4 means termination on an error condition:
* MODE on entry not in [0:1].
* ------------------------------------------------------------------
* This algorithm is due to Richard. P. Brent. It is presented as
* an ALGOL procedure, LOCALMIN, in his 1973 book,
* "Algorithms for minimization without derivatives", Prentice-Hall.
* Published as subroutine FMIN in Forsythe, Malcolm, and Moler,
* "Computer Methods for Mathematical Computations", Prentice-Hall,
* 1977.
* The current subroutine adapted from F.,M.& M. by C. L. Lawson, and
* F. T. Krogh, JPL, Oct 1987, for use in Fortran 77 in the JPL
* MATH77 library. The changes improve performance when a minimum is
* at an endpoint, and change the user interface to reverse
* communication. Unlike the original, the changed algorithm may
* evaluate the function at one of the end points.
* ------------------------------------------------------------------
* Subprograms referenced: R1MACH, IERM1
* ------------------------------------------------------------------
* THE METHOD USED IS A COMBINATION OF GOLDEN SECTION SEARCH AND
* SUCCESSIVE PARABOLIC INTERPOLATION. CONVERGENCE IS NEVER MUCH SLOWER
* THAN THAT FOR A FIBONACCI SEARCH. IF F HAS A CONTINUOUS SECOND
* DERIVATIVE WHICH IS POSITIVE AT THE MINIMUM (WHICH IS NOT AT AX OR
* BX), THEN CONVERGENCE IS SUPERLINEAR, AND USUALLY OF THE ORDER OF
* ABOUT 1.324....
* THE FUNCTION F IS NEVER EVALUATED AT TWO POINTS CLOSER TOGETHER
* THAN EPS*ABS(FMIN) + (TOLI/3), WHERE EPS IS APPROXIMATELY THE SQUARE
* ROOT OF THE RELATIVE MACHINE PRECISION. IF F IS A UNIMODAL
* FUNCTION AND THE COMPUTED VALUES OF F ARE ALWAYS UNIMODAL WHEN
* SEPARATED BY AT LEAST EPS*ABS(X) + (TOLI/3), THEN FMIN APPROXIMATES
* THE ABCISSA OF THE GLOBAL MINIMUM OF F ON THE INTERVAL [AX,BX] WITH
* AN ERROR LESS THAN 3*EPS*ABS(FMIN) + TOLI. IF F IS NOT UNIMODAL,
* THEN FMIN MAY APPROXIMATE A LOCAL, BUT PERHAPS NON-GLOBAL, MINIMUM TO
* THE SAME ACCURACY. (Comments from F., M. & M.)
* ------------------------------------------------------------------
*--S replaces "?": ?FMIN
* Both versions use IERM1
* ------------------------------------------------------------------ */
/* ------------------------------------------------------------------ */
if (*mode < 0)
{
np = -1 - *mode;
return;
}
if (np > 0)
{
np -= 1;
printf(" In SFMIN: X= %g , XORF= %g , IC= %ld \n", *x, *xorf, ic);
}
if (*mode == 1)
{
switch (next)
{
case 1: goto L_301;
case 2: goto L_302;
}
}
if (*mode != 0)
{
/* Error: MODE > 1 on entry. */
ierm1( "SFMIN", *mode, 0, "Input value of MODE exceeds 1."
, "MODE", *mode, '.' );
*mode = 4;
return;
}
/* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* C is the squared inverse of the "Golden Ratio", 1.618...
* C = 0.381966
* */
c = HALF*(THREE - sqrtf( FIVE ));
/* EPS IS THE SQUARE ROOT OF THE RELATIVE MACHINE PRECISION.
* */
if (eps == ZERO)
eps = sqrtf( FLT_EPSILON );
/* INITIALIZATION
* */
a = *x;
b = *xorf;
if (a > b)
{
xi = a;
a = b;
b = xi;
}
/* Now we have A .le. B */
asave = a;
bsave = b;
toli = tol;
if (toli <= ZERO)
toli = ZERO;
c3 = toli/THREE + powif(eps,4);
v = a + c*(b - a);
w = v;
xi = v;
ic = 0;
e = ZERO;
/* Return to calling prog for computation of FX = f(XI) */
*x = xi;
*mode = 1;
next = 1;
return;
/* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
L_301:
;
next = 2;
fx = *xorf;
fv = fx;
fw = fx;
/* MAIN LOOP STARTS HERE */
L_20:
xm = HALF*(a + b);
tol1 = eps*fabsf( xi ) + c3;
tol2 = TWO*tol1;
/* CHECK STOPPING CRITERION
* */
if (fabsf( xi - xm ) <= (tol2 - HALF*(b - a)))
goto L_90;
/* IS GOLDEN-SECTION NECESSARY ?
* */
if (fabsf( e ) <= tol1)
goto L_40;
/* FIT PARABOLA
* */
r = (xi - w)*(fx - fv);
q = (xi - v)*(fx - fw);
p = (xi - v)*q - (xi - w)*r;
q = TWO*(q - r);
if (q > ZERO)
p = -p;
q = fabsf( q );
r = e;
e = d;
/* IS PARABOLA ACCEPTABLE
* */
if (fabsf( p ) >= fabsf( HALF*q*r ))
goto L_40;
if (p <= q*(a - xi))
goto L_40;
if (p >= q*(b - xi))
goto L_40;
/* A PARABOLIC INTERPOLATION STEP
* */
d = p/q;
u = xi + d;
/* F MUST NOT BE EVALUATED TOO CLOSE TO A OR B
* */
if ((u - a) < tol2 || (b - u) < tol2)
d = signf( tol1, xm - xi );
goto L_50;
/* A GOLDEN-SECTION STEP
* */
L_40:
ic += 1;
if (ic > 3)
{
if (a == asave)
{
if (ic == 4)
{
e = a + (eps*fabsf( a ) + c3) - xi;
}
else
{
if (b != w)
goto L_45;
e = a - xi;
}
}
else if (b == bsave)
{
if (ic == 4)
{
e = b - (eps*fabsf( b ) + c3) - xi;
}
else
{
if (a != w)
goto L_45;
e = b - xi;
}
}
else
{
ic = -99;
goto L_45;
}
d = e;
goto L_50;
}
L_45:
if (xi >= xm)
{
e = a - xi;
}
else
{
e = b - xi;
}
d = c*e;
/* F MUST NOT BE EVALUATED TOO CLOSE TO XI
* */
L_50:
if (fabsf( d ) >= tol1)
{
u = xi + d;
}
else
{
if (b - a > THREE*tol1)
tol1 = tol2;
u = fminf( b, fmaxf( a, xi + signf( PT95*tol1, d ) ) );
}
/* Return to calling prog for computation of FU = f(U)
* Returning with MODE = 1 and NEXT = 2 */
*x = u;
return;
/* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
L_302:
;
fu = *xorf;
/* UPDATE A, B, V, W, AND XI
* */
if (fu <= fx)
{
if (fu == fx)
{
if (b - fminf( xi, u ) > fmaxf( xi, u ) - a)
{
b = fmaxf( xi, u );
}
else
{
a = fminf( xi, u );
}
}
else
{
if (u >= xi)
{
a = xi;
}
else
{
b = xi;
}
}
v = w;
fv = fw;
w = xi;
fw = fx;
xi = u;
fx = fu;
goto L_20;
}
if (u < xi)
{
a = u;
}
else
{
b = u;
}
if (fu <= fw || w == xi)
{
v = w;
fv = fw;
w = u;
fw = fu;
}
else if ((fu <= fv || v == xi) || v == w)
{
v = u;
fv = fu;
}
goto L_20;
/* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
L_90:
;
*x = xi;
*xorf = fx;
*mode = 2;
if ((b - a > THREE*toli) && (toli != ZERO))
*mode = 3;
return;
} /* end of function */