I am a very bad programmer, so beware of the high probability that I have done a very stupid mistake. I am writing C and using gsl_multiroot_fsolver_hybrid wrong, any help will be much apprteciated!
My code contains the following two functions:
/* --- --- */
int quasilognormal_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
const double x0 = gsl_vector_get (x, 0);
const double x1 = gsl_vector_get (x, 1);
const double x2 = gsl_vector_get (x, 2);
qparams->N =x0;
qparams->A =x1;
qparams->S =x2;
printf("INNER POINT 1 -- Attention!\n");
printf("N,A,S = %e %e %e \n",
x0, x1, x2);
const double y0 = norm_quasilognormal(qparams) - 1.0;
//const double y0 = 0;
printf("INNER POINT 2\n");
const double y1 = mean_quasilognormal(qparams) - 0.0;
//const double y1 =0;
printf("INNER POINT 3\n");
const double y2 = var_quasilognormal(qparams) - kappa_var(gcosmo-> theta);
printf("INNER POINT 4\n");
gsl_vector_set (f, 0, y0);
gsl_vector_set (f, 1, y1);
gsl_vector_set (f, 2, y2);
return GSL_SUCCESS;
}
int COMPUTE_QUASILOGNORMAL_PARAMETERS()
{
gsl_vector * root;
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
int status;
size_t i, iter = 0;
const size_t n = 3;
gsl_multiroot_function F;
F.f = &quasilognormal_f;
F.params = NULL;
F.n = 3;
double x_init[3] = {1.0, 0.0,gcosmo->ausgelagertes_lognormal_sigma};
gsl_vector *x = gsl_vector_alloc (n);
gsl_vector_set (x, 0, x_init[0]);
gsl_vector_set (x, 1, x_init[1]);
gsl_vector_set (x, 2, x_init[2]);
T = gsl_multiroot_fsolver_hybrid;
s = gsl_multiroot_fsolver_alloc (T, 3);
gsl_multiroot_fsolver_set (s, &F, x);
printf("THE SOLVER WAS SET!\n");
do
{
iter++;
printf("ABOUT TO ITERATE SOLVER\n");
printf("iter = %i\n",
iter);
printf("x_init = %e %e %e \n",
x_init[0], x_init[1], x_init[2]);
printf("N,A,S = %e %e %e \n",
root[0], root[1], root[2]);
status = gsl_multiroot_fsolver_iterate (s);
printf("THE SOLVER WAS ITERATED\n");
printf("iter = %i\n",
iter);
printf("x_init = %e %e %e \n",
x_init[0], x_init[1], x_init[2]);
printf("N,A,S = %e %e %e \n",
root[0], root[1], root[2]);
if (status) /* check if solver is stuck */
break;
status =
gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
root = gsl_multiroot_fsolver_root(s);
printf("x = % .3e \n",
root[0], root[1], root[2]);
gsl_multiroot_fsolver_free (s);
gsl_vector_free (x);
return 0;
}
/* --- --- */
The function int COMPUTE_QUASILOGNORMAL_PARAMETERS should solve a system of three equations that are defined in int quasilognormal_f. This goes wrong, because as soon as the second (iter ==1) iteration starts, the parameters N, A and S are NAN inside int quasilognormal_f. This I concluded from the programme output:
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 1.490116e-08 2.071317e-01
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071318e-01
INNER POINT 2
INNER POINT 3
INNER POINT 4
THE SOLVER WAS SET!
ABOUT TO ITERATE SOLVER
iter = 1
x_init = 1.000000e+00 0.000000e+00 2.071317e-01
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01
INNER POINT 1 -- Attention!
N,A,S = nan nan nan
I don't know what to do because the solver is a blackbox to me and I tried following the GSL manual.
