Dear listmembers I have struggeled with the "gsl_multifit_fdfsolver_lmsder", and it works fine *once*. I have tried to look over the source of "gsl_multifit_", but haven't found anything. I still think that it's probably my own routine, but it follows the example ing gsl-ref closely... The "lmsder" converges correctly the firste time, but when the routine is called again [with the same parameters] it fails, and sends X= ['nan' 'nan' 'nan'] to the solverfunctions f and df. I do allocate space for the solver, initialize the solver, and I do free the solver afterwards. And I cant find any reason for the behaivour of the solver. 3 calls as the following fails... [And I dont get it :-)] solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); /* OK */ solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); /* FAILS */ solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); /* FAILS */ ANy insight will be much appreciated. The problem solved is n>=3 equations: 0 = P0 - X0*sin(X1*P1+X2) Regards Nicolai Hanssing Denmark I have included a compilable listing of the messy source: ------------ fit_Asterix.c: -- /* * fit_Asterix.h for use in Asterix_main.c * * Nicolai Hanssing, nuh@kampsax.dtu.dk * * Masters Thesis 01, Foss-Electric * */ #include #include #include #include #include #include "fit_Asterix.h" /* "Private" prototyper */ int ulineval_f(const gsl_vector *x, void *params, gsl_vector *f); int ulineval_df (const gsl_vector * x, void *params,gsl_matrix *df); int ulineval_fdf (const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *df); int print_state(size_t iter, gsl_multifit_fdfsolver * s); struct data { double L[5]; gsl_vector *v_K; }; #define DEBUG #define DEBUG_DELRIN #define DEBUG_SNS int solve_nonlin_sys_multifit_nlin(gsl_vector *v_topp,gsl_vector *v_alpha){ /* Løser det overbestemte ulineære system * med: 0 = L - a3sin(a4*k+a5) */ const gsl_multifit_fdfsolver_type * T; gsl_multifit_fdfsolver *s; gsl_multifit_function_fdf FDF; struct data params; int status, i; size_t iter = 0; const size_t n = 3; /* number of equations >= 3 */ const size_t p = 3; /* number of parameters (3) */ double x_init[3] = { 0.1600, 0.0, 0.0 }; gsl_vector_view x = gsl_vector_view_array (x_init, p); params.v_K = gsl_vector_alloc(v_topp->size); gsl_vector_memcpy(params.v_K,v_topp); gsl_vector_scale (params.v_K, 0.01); /* To avoid badly conditioned Jacobian */ for (i=0;i<5;i++) params.L[i] = Delrin_lambda[i]/10000.0; /* To avoid badly conditioned Jacobian */ FDF.n = n; /* 3 equations */ FDF.p = p; /* 3 variables */ FDF.params = (void *) ¶ms; FDF.f = &ulineval_f; FDF.df = &ulineval_df; FDF.fdf = &ulineval_fdf; T = gsl_multifit_fdfsolver_lmsder; s = gsl_multifit_fdfsolver_alloc (T, n, p); status = gsl_multifit_fdfsolver_set (s, &FDF, &x.vector); #ifdef DEBUG_SNS printf ("\n status = %s\n", gsl_strerror (status)); print_state (iter, s); #endif do { iter++; status = gsl_multifit_fdfsolver_iterate (s); #ifdef DEBUG_SNS if (!(int)fmod(iter,1)) { printf ("\n status = %s\n", gsl_strerror (status)); print_state (iter, s); } #endif if (status) break; status = gsl_multifit_test_delta (s->dx, s->x, 0.000001, 0.000001); } while (status == GSL_CONTINUE && iter < 500); /* Scaling back */ gsl_vector_set(v_alpha,2,gsl_vector_get(s->x,0)*10000); gsl_vector_set(v_alpha,3,gsl_vector_get(s->x,1)/100); gsl_vector_set(v_alpha,4,gsl_vector_get(s->x,2)); gsl_vector_free(params.v_K); gsl_multifit_fdfsolver_free(s); return GSL_SUCCESS; } int ulineval_f(const gsl_vector * x, void *params, gsl_vector * f){ size_t i; double f_val; long double Alpha3=gsl_vector_get(x,0); long double Alpha4=gsl_vector_get(x,1); long double Alpha5=gsl_vector_get(x,2); long double Ki,L; #ifdef DEBUG_SNS printf(" X = %3.6e %3.6e %3.6e\n ", gsl_vector_get(x,0), gsl_vector_get(x,1), gsl_vector_get(x,2)); #endif for (i=0;i<3;i++){ Ki = gsl_vector_get(((struct data *)params)->v_K,i); L = ((struct data *)params)->L[i]; f_val= L - Alpha3 * sin(Alpha4 * Ki + Alpha5); gsl_vector_set(f,i,f_val); #ifdef DEBUG_SNS printf("%e ",f_val); #endif } #ifdef DEBUG_SNS printf("\n"); #endif return GSL_SUCCESS; } int ulineval_df (const gsl_vector * x, void *params,gsl_matrix *df){ size_t i,h; double Alpha3=gsl_vector_get(x,0); double Alpha4=gsl_vector_get(x,1); double Alpha5=gsl_vector_get(x,2); double j[3]; double Ki; for (i=0;i<3;i++){ Ki=gsl_vector_get(((struct data *)params)->v_K,i); j[0] = - sin(Alpha4 * Ki + Alpha5); j[1] = - Alpha3 * cos(Alpha4 * Ki + Alpha5)*Ki; j[2] = - Alpha3 * cos(Alpha4 * Ki + Alpha5); for (h=0;h<3;h++) gsl_matrix_set (df, i, h, j[h]); #ifdef DEBUG_SNS printf("\n | %e %e %e | ",j[0],j[1],j[2]); #endif } #ifdef DEBUG_SNS printf("\n "); #endif return GSL_SUCCESS; } int ulineval_fdf (const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *df) { ulineval_f(x, params, f); ulineval_df(x, params, df); return GSL_SUCCESS; } int print_state(size_t iter, gsl_multifit_fdfsolver * s){ printf ("\niter: %3u x = % 15.8f % 15.8f % 15.8f |f(x)| = %g\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->x, 2), gsl_blas_dnrm2 (s->f)); } -- main.c -- /* Testprogram for debugging lmsder */ #include "fit_Asterix.h" #include int main(void){ double topp[] = {525.80226278891245783598, 620.25431012228921190399, 725.58581627818944070896, 854.57557708963383902301, 960.14100685986863936705}; gsl_vector_view v_topp; gsl_vector *v_alpha; v_alpha=gsl_vector_alloc(5); v_topp = gsl_vector_view_array(topp, 5); printf("\n\n Calling solver\n"); solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); /* The program fails on the following calls */ printf("\n\n\n Calling solver, again...\n"); solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); printf("\n\n\n Calling solver, again...\n"); solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha); gsl_vector_free(v_alpha); return 0; } -- fit_Asterix.h: -- #define Topp_samples 7 /* Antal mÃ¥linger der skal * bruges ved fit af toppunkt * OBS: Skal være ulige */ #define Delrin_MinMax 5 /* Antal kendte lokale * min/max punkter: */ static const double Delrin_lambda[]={903.408416839855, 960.556012816831, 1022.49708053353, 0, 0}; /* i nm */ static const double Delrin_Arb[]={ 1.61981074275959, 1.5127518464565, 1.59433851771899, 0, 0}; /* OBS - Fundet via splinefit.... [wb04.m] */ /* 1.61981074275959 1.5127518464565 1.59433851771899 Arb */ /* 903.408416839855 960.556012816831 1022.49708053353 lambda */ #include int solve_nonlin_sys_multifit_nlin(gsl_vector *v_topp,gsl_vector *v_alpha); -- fit_Asterix.c