Error in multifit_lmsder (?)

[Nicolai Hanssing <nuh at kampsax dot dtu dot dk>]

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 <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <math.h>
#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 *) &params;	
	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 <gsl/gsl_vector.h>


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 <gsl/gsl_vector.h>

int solve_nonlin_sys_multifit_nlin(gsl_vector *v_topp,gsl_vector
*v_alpha);
--

fit_Asterix.c