Hi,I am using the root finding algorithms (mainly brent), and I realized that it's not possible (at least not straightforwardly intended)
1) to provide function values when initializing the root solver or 2) to retrieve the function values at x_lower and x_upper after an iteration directly from the methods.
I attach a patch (as a diff to commit e38f2146) which implements the first of those two things I mentioned. Basically, besides the existing function gsl_root_fsolver_set(...) it provides another such function gsl_root_fsolver_set2(..., f_lower, f_upper) which allows to also provide the function values at x_lower and x_upper and in return skips the function evaluation in the initialization routine.
I don't think there is another way to implement this because my main intention is to get rid of the two (in my case costly) function evaluations.
Explanation why I want this:Evaluating my function f(x), of which I want to find a root, takes a few minutes time (here, f(x) is a characteristic value of the solution of an elliptic PDE). Furthermore, I don't know a priori what a good choice of x_lower and x_upper is. I have to iterate first through some generic interval before I can start the Brent search; this means that I have already calculated f(x_lower) and f(x_upper). The Brent initialization would calculate those values again and waste a few minutes time.
In the same way, I will write a patch that allows me to access the function value at the end of the Brent iteration, rather than evaluating f(x) yet another time. In the end, those changes would allow me to save 20% - 30% computational time.
If this patch is any interesting, I'm happy for it to be implemented in GSL and I'm also happy to make changes (e.g. a nicer name as *set2() is not very inventive or do the same for the fdfsolver methods) beforehand.
Christian
diff --git a/roots/bisection.c b/roots/bisection.c
index 9665b736..6905b221 100644
--- a/roots/bisection.c
+++ b/roots/bisection.c
@@ -41,6 +41,7 @@ typedef struct
bisection_state_t;
static int bisection_init (void * vstate, gsl_function * f, double * root,
double x_lower, double x_upper);
+static int bisection_init2 (void * vstate, gsl_function * f, double * root,
double x_lower, double f_lower, double x_upper, double f_upper);
static int bisection_iterate (void * vstate, gsl_function * f, double * root,
double * x_lower, double * x_upper);
static int
@@ -67,6 +68,25 @@ bisection_init (void * vstate, gsl_function * f, double *
root, double x_lower,
}
+static int
+bisection_init2 (void * vstate, gsl_function * f, double * root, double
x_lower, double f_lower, double x_upper, double f_upper)
+{
+ bisection_state_t * state = (bisection_state_t *) vstate;
+
+ *root = 0.5 * (x_lower + x_upper) ;
+
+ state->f_lower = f_lower;
+ state->f_upper = f_upper;
+
+ if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+ {
+ GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+ }
+
+ return GSL_SUCCESS;
+
+}
+
static int
bisection_iterate (void * vstate, gsl_function * f, double * root, double *
x_lower, double * x_upper)
{
@@ -129,6 +149,7 @@ static const gsl_root_fsolver_type bisection_type =
{"bisection", /* name */
sizeof (bisection_state_t),
&bisection_init,
+ &bisection_init2,
&bisection_iterate};
const gsl_root_fsolver_type * gsl_root_fsolver_bisection = &bisection_type;
diff --git a/roots/brent.c b/roots/brent.c
index 4f0c9cd3..c2a81b09 100644
--- a/roots/brent.c
+++ b/roots/brent.c
@@ -42,6 +42,7 @@ typedef struct
brent_state_t;
static int brent_init (void * vstate, gsl_function * f, double * root, double
x_lower, double x_upper);
+static int brent_init2 (void * vstate, gsl_function * f, double * root, double
x_lower, double f_lower, double x_upper, double f_upper);
static int brent_iterate (void * vstate, gsl_function * f, double * root,
double * x_lower, double * x_upper);
@@ -69,6 +70,38 @@ brent_init (void * vstate, gsl_function * f, double * root,
double x_lower, doub
state->d = x_upper - x_lower ;
state->e = x_upper - x_lower ;
+ printf("# BRENT INIT with %e %e -- %e %e\n", x_lower, f_lower, x_upper,
f_upper);
+
+ if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+ {
+ GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+ }
+
+ return GSL_SUCCESS;
+
+}
+
+static int
+brent_init2 (void * vstate, gsl_function * f, double * root, double x_lower,
double f_lower, double x_upper, double f_upper)
+{
+ brent_state_t * state = (brent_state_t *) vstate;
+
+ *root = 0.5 * (x_lower + x_upper) ;
+
+ state->a = x_lower;
+ state->fa = f_lower;
+
+ state->b = x_upper;
+ state->fb = f_upper;
+
+ state->c = x_upper;
+ state->fc = f_upper;
+
+ state->d = x_upper - x_lower ;
+ state->e = x_upper - x_lower ;
+
+ printf("# BRENT INIT with %e %e -- %e %e\n", x_lower, f_lower, x_upper,
f_upper);
+
if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
{
GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
@@ -239,6 +272,7 @@ static const gsl_root_fsolver_type brent_type =
{"brent", /* name */
sizeof (brent_state_t),
&brent_init,
+ &brent_init2,
&brent_iterate};
const gsl_root_fsolver_type * gsl_root_fsolver_brent = &brent_type;
diff --git a/roots/falsepos.c b/roots/falsepos.c
index f6cda5c8..95080173 100644
--- a/roots/falsepos.c
+++ b/roots/falsepos.c
@@ -52,6 +52,7 @@ typedef struct
falsepos_state_t;
static int falsepos_init (void * vstate, gsl_function * f, double * root,
double x_lower, double x_upper);
+static int falsepos_init2 (void * vstate, gsl_function * f, double * root,
double x_lower, double f_lower, double x_upper, double f_upper);
static int falsepos_iterate (void * vstate, gsl_function * f, double * root,
double * x_lower, double * x_upper);
static int
@@ -78,6 +79,25 @@ falsepos_init (void * vstate, gsl_function * f, double *
root, double x_lower, d
}
+static int
+falsepos_init2 (void * vstate, gsl_function * f, double * root, double
x_lower, double f_lower, double x_upper, double f_upper)
+{
+ falsepos_state_t * state = (falsepos_state_t *) vstate;
+
+ *root = 0.5 * (x_lower + x_upper);
+
+ state->f_lower = f_lower;
+ state->f_upper = f_upper;
+
+ if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
+ {
+ GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
+ }
+
+ return GSL_SUCCESS;
+
+}
+
static int
falsepos_iterate (void * vstate, gsl_function * f, double * root, double *
x_lower, double * x_upper)
{
@@ -173,6 +193,7 @@ static const gsl_root_fsolver_type falsepos_type =
{"falsepos", /* name */
sizeof (falsepos_state_t),
&falsepos_init,
+ &falsepos_init2,
&falsepos_iterate};
const gsl_root_fsolver_type * gsl_root_fsolver_falsepos = &falsepos_type;
diff --git a/roots/fsolver.c b/roots/fsolver.c
index 8f586bb6..8d4b7ad4 100644
--- a/roots/fsolver.c
+++ b/roots/fsolver.c
@@ -66,6 +66,22 @@ gsl_root_fsolver_set (gsl_root_fsolver * s, gsl_function *
f, double x_lower, do
return (s->type->set) (s->state, s->function, &(s->root), x_lower, x_upper);
}
+int
+gsl_root_fsolver_set2 (gsl_root_fsolver * s, gsl_function * f, double x_lower,
double f_lower, double x_upper, double f_upper)
+{
+ if (x_lower > x_upper)
+ {
+ GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
+ }
+
+ s->function = f;
+ s->root = 0.5 * (x_lower + x_upper); /* initial estimate */
+ s->x_lower = x_lower;
+ s->x_upper = x_upper;
+
+ return (s->type->set2) (s->state, s->function, &(s->root), x_lower, f_lower,
x_upper, f_upper);
+}
+
int
gsl_root_fsolver_iterate (gsl_root_fsolver * s)
{
diff --git a/roots/gsl_roots.h b/roots/gsl_roots.h
index 46e45870..fef4e9e8 100644
--- a/roots/gsl_roots.h
+++ b/roots/gsl_roots.h
@@ -41,6 +41,7 @@ typedef struct
const char *name;
size_t size;
int (*set) (void *state, gsl_function * f, double * root, double x_lower,
double x_upper);
+ int (*set2) (void *state, gsl_function * f, double * root, double x_lower,
double f_lower, double x_upper, double f_upper);
int (*iterate) (void *state, gsl_function * f, double * root, double *
x_lower, double * x_upper);
}
gsl_root_fsolver_type;
@@ -82,6 +83,11 @@ int gsl_root_fsolver_set (gsl_root_fsolver * s,
gsl_function * f,
double x_lower, double x_upper);
+int gsl_root_fsolver_set2 (gsl_root_fsolver * s,
+ gsl_function * f,
+ double x_lower, double f_lower,
+ double x_upper, double f_upper);
+
int gsl_root_fsolver_iterate (gsl_root_fsolver * s);
const char * gsl_root_fsolver_name (const gsl_root_fsolver * s);
OpenPGP_0x5D5587E0B7C49590.asc
Description: application/pgp-keys
OpenPGP_signature
Description: OpenPGP digital signature
