diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp index fe2832789c81cdc52e32b9cb4bb7ac712a29fe41..d8ca75020602196a93346b4363176966b8a7d768 100644 --- a/Geo/Geo.cpp +++ b/Geo/Geo.cpp @@ -2874,11 +2874,9 @@ static Vertex *VERTEX; static double min1d(double (*funct) (double), double *xmin) { - // we should think about the tolerance more carefully... - double ax = 1.e-15, bx = 1.e-12, cx = 1.e-11, fa, fb, fx, tol = 1.e-4; - mnbrak(&ax, &bx, &cx, &fa, &fx, &fb, funct); - //Msg::Info("--MIN1D : ax %12.5E bx %12.5E cx %12.5E",ax,bx,cx); - return (brent(ax, bx, cx, funct, tol, xmin)); + // 0. for tolerance allows for maximum as code in gsl_brent + return (brent(CURVE->ubeg, 0.5*(CURVE->ubeg + CURVE->uend), CURVE->uend, + funct, 0., xmin)); } static void projectPS(int N, double x[], double res[]) @@ -2900,10 +2898,6 @@ static void projectPS(int N, double x[], double res[]) static double projectPC(double u) { - if(u < CURVE->ubeg) - u = CURVE->ubeg; - if(u < CURVE->ubeg) - u = CURVE->ubeg; Vertex c = InterpolateCurve(CURVE, u, 0); return sqrt(SQU(c.Pos.X - VERTEX->Pos.X) + SQU(c.Pos.Y - VERTEX->Pos.Y) + diff --git a/Numeric/gsl_brent.cpp b/Numeric/gsl_brent.cpp index 3d18a68377d1d44bcfaa7dfe2b5869846e9b69b4..3cabe77250589a3695d7db8ed5c3cc5f4c80f6ec 100644 --- a/Numeric/gsl_brent.cpp +++ b/Numeric/gsl_brent.cpp @@ -20,6 +20,29 @@ double fn1(double x, void *params) return val; } +int gsl_min_fminimizer_set_with_values_MOD(gsl_min_fminimizer *s, + gsl_function *f, + double x_minimum, double f_minimum, + double x_lower, double f_lower, + double x_upper, double f_upper) +{ + s->function = f; + s->x_minimum = x_minimum; + s->x_lower = x_lower; + s->x_upper = x_upper; + if (x_lower > x_upper) + { + GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL); + } + s->f_lower = f_lower; + s->f_upper = f_upper; + s->f_minimum = f_minimum; + return (s->type->set) (s->state, s->function, + x_minimum, f_minimum, + x_lower, f_lower, + x_upper, f_upper); +} + #define MAXITER 100 // Returns the minimum betwen ax and cx to a given tolerance tol using @@ -28,9 +51,14 @@ double fn1(double x, void *params) double brent(double ax, double bx, double cx, double (*f) (double), double tol, double *xmin) { + // The solver can stall at the following internal tolerance in brent - so + // this is about the lowest possible tolerance. + if(tol == 0.) tol = 10.*GSL_SQRT_DBL_EPSILON; + const double tolsq = tol*tol; int status; int iter = 0; double a, b, c; // a < b < c + double fa, fb, fc; const gsl_min_fminimizer_type *T; gsl_min_fminimizer *s; gsl_function F; @@ -42,38 +70,56 @@ double brent(double ax, double bx, double cx, c = cx; } else { - a = ax; - c = cx; - } - - // if a-b < tol, return func(a) - if(fabs(c - a) < tol) { - *xmin = ax; - return (f(*xmin)); + a = cx; + c = ax; } nrfunc = f; - F.function = &fn1; F.params = 0; + fa = f(a); + fb = f(b); + fc = f(c); + if(fb > fa || fb > fc) { + if(a < c) { + b = a; + fb = fa; + } + else { + b = c; + fb = fc; + } + } + + // if a-c < tol, return func(b) + if(gsl_min_test_interval(a, c, tolsq, tol) == GSL_SUCCESS) { + *xmin = b; + return (fb); + } + T = gsl_min_fminimizer_brent; s = gsl_min_fminimizer_alloc(T); - gsl_min_fminimizer_set(s, &F, b, a, c); +// gsl_min_fminimizer_set(s, &F, b, a, c); // Standard + // This mod avoids some error checks so we can have a minimum at an extent. + gsl_min_fminimizer_set_with_values_MOD(s, &F, b, fb, a, fa, c, fc); do { iter++; status = gsl_min_fminimizer_iterate(s); if(status) break; // solver problem - b = gsl_min_fminimizer_minimum(s); - // this is deprecated: we should use gsl_min_fminimizer_x_minimum(s) instead a = gsl_min_fminimizer_x_lower(s); c = gsl_min_fminimizer_x_upper(s); - // we should think about the tolerance more carefully... - status = gsl_min_test_interval(a, c, tol, tol); +// fb = gsl_min_fminimizer_f_minimum(s); +// fa = gsl_min_fminimizer_f_lower(s); +// fc = gsl_min_fminimizer_f_upper(s); + status = gsl_min_test_interval(a, c, tolsq, tol); + // ... in case curve becomes too flat - not required? +// std::fabs(fa - fc) < (std::fabs(fb) + tol)*tol } while(status == GSL_CONTINUE && iter < MAXITER); + b = gsl_min_fminimizer_x_minimum(s); if(status != GSL_SUCCESS) Msg::Error("MIN1D not converged: f(%g) = %g", b, fn1(b, NULL));