"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Numeric/NumericUtils.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

NumericUtils.cpp  (getdp-3.4.0-source.tgz)	:	NumericUtils.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege		// GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
//		//
// See the LICENSE.txt file for license information. Please report all		// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues.		// issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include "GetDPConfig.h"		#include "GetDPConfig.h"
#include "Message.h"		#include "Message.h"
#if !defined(HAVE_GSL) && !defined(HAVE_NR)		#if !defined(HAVE_GSL) && !defined(HAVE_NR)
double brent(double ax, double bx, double cx,		double brent(double ax, double bx, double cx, double (*f)(double), double tol,
double (*f) (double), double tol, double *xmin)		double *xmin)
{		{
Message::Error("Minimization routines require GSL or NR");		Message::Error("Minimization routines require GSL or NR");
return 0;		return 0;
}		}
void mnbrak(double *ax, double *bx, double *cx,		void mnbrak(double *ax, double *bx, double *cx, double *fa_dummy,
double *fa_dummy, double *fb_dummy, double *fc_dummy,		double *fb_dummy, double *fc_dummy, double (*func)(double))
double (*func) (double))
{		{
Message::Error("Minimization routines require GSL or NR");		Message::Error("Minimization routines require GSL or NR");
}		}
#endif		#endif
#if defined(HAVE_GSL)		#if defined(HAVE_GSL)
#include <gsl/gsl_errno.h>		#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>		#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>		#include <gsl/gsl_min.h>
static double (*nrfunc) (double);		static double (*nrfunc)(double);
double fn1(double x, void *params)		double fn1(double x, void *params)
{		{
double val = nrfunc(x);		double val = nrfunc(x);
return val;		return val;
}		}
#define MAXITER 100		#define MAXITER 100
// Returns the minimum betwen ax and cx to a given tolerance tol using		// Returns the minimum betwen ax and cx to a given tolerance tol using
// brent's method.		// brent's method.
double brent(double ax, double bx, double cx,		double brent(double ax, double bx, double cx, double (*f)(double), double tol,
double (*f) (double), double tol, double *xmin)		double *xmin)
{		{
int status;		int status;
int iter = 0;		int iter = 0;
double a, b, c; // a < b < c		double a, b, c; // a < b < c
const gsl_min_fminimizer_type *T;		const gsl_min_fminimizer_type *T;
gsl_min_fminimizer *s;		gsl_min_fminimizer *s;
gsl_function F;		gsl_function F;
// gsl wants a<b		// gsl wants a<b
b = bx;		b = bx;
if(ax < cx) {		if(ax < cx) {
a = ax;		a = ax;
c = cx;		c = cx;
}		}
skipping to change at line 85		skipping to change at line 84
F.function = &fn1;		F.function = &fn1;
F.params = 0;		F.params = 0;
T = gsl_min_fminimizer_brent;		T = gsl_min_fminimizer_brent;
s = gsl_min_fminimizer_alloc(T);		s = gsl_min_fminimizer_alloc(T);
gsl_min_fminimizer_set(s, &F, b, a, c);		gsl_min_fminimizer_set(s, &F, b, a, c);
do {		do {
iter++;		iter++;
status = gsl_min_fminimizer_iterate(s);		status = gsl_min_fminimizer_iterate(s);
if(status)		if(status) break; // solver problem
break; // solver problem
b = gsl_min_fminimizer_minimum(s);		b = gsl_min_fminimizer_minimum(s);
// this is deprecated: we should use gsl_min_fminimizer_x_minimum(s) instead		// this is deprecated: we should use gsl_min_fminimizer_x_minimum(s) instead
a = gsl_min_fminimizer_x_lower(s);		a = gsl_min_fminimizer_x_lower(s);
c = gsl_min_fminimizer_x_upper(s);		c = gsl_min_fminimizer_x_upper(s);
// we should think about the tolerance more carefully...		// we should think about the tolerance more carefully...
status = gsl_min_test_interval(a, c, tol, tol);		status = gsl_min_test_interval(a, c, tol, tol);
}		} while(status == GSL_CONTINUE && iter < MAXITER);
while(status == GSL_CONTINUE && iter < MAXITER);
if(status != GSL_SUCCESS)		if(status != GSL_SUCCESS)
Message::Error("MIN1D not converged: f(%g) = %g", b, fn1(b, NULL));		Message::Error("MIN1D not converged: f(%g) = %g", b, fn1(b, NULL));
*xmin = b;		*xmin = b;
gsl_min_fminimizer_free(s);		gsl_min_fminimizer_free(s);
return fn1(b, NULL);		return fn1(b, NULL);
}		}
// Find an initial bracketting of the minimum of func. Given 2 initial		// Find an initial bracketting of the minimum of func. Given 2 initial
// points ax and bx, mnbrak checks in which direction func decreases,		// points ax and bx, mnbrak checks in which direction func decreases,
// and takes some steps in that direction, until the function		// and takes some steps in that direction, until the function
// increases--at cx. mnbrak returns ax and cx (possibly switched),		// increases--at cx. mnbrak returns ax and cx (possibly switched),
// bracketting a minimum.		// bracketting a minimum.
#define MYGOLD_ 1.618034		#define MYGOLD_ 1.618034
#define MYLIMIT_ 100.0		#define MYLIMIT_ 100.0
#define MYTINY_ 1.0e-20		#define MYTINY_ 1.0e-20
#define SIGN(a,b)((b) >= 0.0 ? fabs(a) : -fabs(a))		#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void mnbrak(double *ax, double *bx, double *cx,		void mnbrak(double *ax, double *bx, double *cx, double *fa_dummy,
double *fa_dummy, double *fb_dummy, double *fc_dummy,		double *fb_dummy, double *fc_dummy, double (*func)(double))
double (*func) (double))
{		{
double ulim, u, r, q;		double ulim, u, r, q;
volatile double f_a, f_b, f_c, f_u;		volatile double f_a, f_b, f_c, f_u;
f_a = (*func) (*ax);		f_a = (*func)(*ax);
f_b = (*func) (*bx);		f_b = (*func)(*bx);
if(f_b > f_a) {		if(f_b > f_a) {
double tmp;		double tmp;
tmp = *ax;		tmp = *ax;
*ax = *bx;		*ax = *bx;
*bx = tmp;		*bx = tmp;
tmp = f_b;		tmp = f_b;
f_b = f_a;		f_b = f_a;
f_a = tmp;		f_a = tmp;
}		}
*cx = *bx + MYGOLD_ * (*bx - *ax);		*cx = *bx + MYGOLD_ * (*bx - *ax);
f_c = (*func) (*cx);		f_c = (*func)(*cx);
while(f_b > f_c) {		while(f_b > f_c) {
r = (*bx - *ax) * (f_b - f_c);		r = (*bx - *ax) * (f_b - f_c);
q = (*bx - *cx) * (f_b - f_a);		q = (*bx - *cx) * (f_b - f_a);
u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /		u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
(2.0 * SIGN(std::max(fabs(q - r), MYTINY_), q - r));		(2.0 * SIGN(std::max(fabs(q - r), MYTINY_), q - r));
ulim = *bx + MYLIMIT_ * (*cx - *bx);		ulim = *bx + MYLIMIT_ * (*cx - *bx);
if((*bx - u) * (u - *cx) > 0.0) {		if((*bx - u) * (u - *cx) > 0.0) {
f_u = (*func) (u);		f_u = (*func)(u);
if(f_u < f_c) {		if(f_u < f_c) {
*ax = *bx;		*ax = *bx;
*bx = u;		*bx = u;
return;		return;
}		}
else if(f_u > f_b) {		else if(f_u > f_b) {
*cx = u;		*cx = u;
return;		return;
}		}
u = *cx + MYGOLD_ * (*cx - *bx);		u = *cx + MYGOLD_ * (*cx - *bx);
f_u = (*func) (u);		f_u = (*func)(u);
}		}
else if((*cx - u) * (u - ulim) > 0.0) {		else if((*cx - u) * (u - ulim) > 0.0) {
f_u = (*func) (u);		f_u = (*func)(u);
if(f_u < f_c) {		if(f_u < f_c) {
*bx = *cx;		*bx = *cx;
*cx = u;		*cx = u;
u = *cx + MYGOLD_ * (*cx - *bx);		u = *cx + MYGOLD_ * (*cx - *bx);
f_b = f_c;		f_b = f_c;
f_c = f_u;		f_c = f_u;
f_u = (*func) (u);		f_u = (*func)(u);
}		}
}		}
else if((u - ulim) * (ulim - *cx) >= 0.0) {		else if((u - ulim) * (ulim - *cx) >= 0.0) {
u = ulim;		u = ulim;
f_u = (*func) (u);		f_u = (*func)(u);
}		}
else {		else {
u = *cx + MYGOLD_ * (*cx - *bx);		u = *cx + MYGOLD_ * (*cx - *bx);
f_u = (*func) (u);		f_u = (*func)(u);
}		}
*ax = *bx;		*ax = *bx;
*bx = *cx;		*bx = *cx;
*cx = u;		*cx = u;
f_a = f_b;		f_a = f_b;
f_b = f_c;		f_b = f_c;
f_c = f_u;		f_c = f_u;
}		}
}		}
End of changes. 20 change blocks.
30 lines changed or deleted		26 lines changed or added

---