diff --git a/Numeric/gsl_newt.cpp b/Numeric/gsl_newt.cpp index bba2b52eca1542b4313a6a25f25931d17bc5bc00..34a5167f368af915208b1b20681827a980545c99 100644 --- a/Numeric/gsl_newt.cpp +++ b/Numeric/gsl_newt.cpp @@ -1,4 +1,4 @@ -// $Id: gsl_newt.cpp,v 1.5 2003-02-20 16:49:21 geuzaine Exp $ +// $Id: gsl_newt.cpp,v 1.6 2003-02-20 17:02:00 geuzaine Exp $ // // Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle // @@ -37,10 +37,10 @@ #define MAXITER 10000 #define PREC 1.e-8 -int gsl_dim; -double gsl_u[MAX_DIM_NEWT], gsl_v[MAX_DIM_NEWT]; -static void (*nrfunc)(int n, double x[],double y[]); -struct gsl_dummy{;}; +int nrdim; +double nru[MAX_DIM_NEWT], nrv[MAX_DIM_NEWT]; +static void (*nrfunc)(int n, double x[],double y[]); +struct gsl_dummy{;}; void convert_vector_from_gsl(const gsl_vector *gv, double * v){ int i, m; @@ -58,15 +58,15 @@ void convert_vector_to_gsl(double *v, int n, gsl_vector *gv){ } int gslfunc(const gsl_vector *xx, void *params, gsl_vector *f){ - convert_vector_from_gsl(xx,gsl_u); - (*nrfunc)(gsl_dim,gsl_u,gsl_v); - // Msg(INFO, "f(%lf,%lf) = %lf %lf\n",gsl_u[1],gsl_u[2],gsl_v[1],gsl_v[2]); - convert_vector_to_gsl(gsl_v,gsl_dim,f); + convert_vector_from_gsl(xx,nru); + (*nrfunc)(nrdim,nru,nrv); + // Msg(INFO, "f(%lf,%lf) = %lf %lf\n",nru[1],nru[2],nrv[1],nrv[2]); + convert_vector_to_gsl(nrv,nrdim,f); return GSL_SUCCESS; } -// Warning: for compatibility with the old newt from NR, x[] is -// indexed from 1 to N! +// Warning: for compatibility with the old newt from NR, x[] and the +// arguments of func() are indexed from 1 to n... void newt(double x[], int n, int *check, void (*func)(int, double [],double [])){ const gsl_multiroot_fsolver_type *T; @@ -78,7 +78,7 @@ void newt(double x[], int n, int *check, void (*func)(int, double [],double [])) gsl_vector *xx = gsl_vector_alloc (n); if(n > MAX_DIM_NEWT-1) Msg(FATAL, "Maximum Newton dimension exceeded\n"); - gsl_dim = n; + nrdim = n; nrfunc = func; convert_vector_to_gsl(x,n,xx); @@ -90,22 +90,20 @@ void newt(double x[], int n, int *check, void (*func)(int, double [],double [])) do{ iter++; status = gsl_multiroot_fsolver_iterate(s); - // Msg(INFO, "status %d %d %d %lf %lf\n",status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); + // Msg(INFO, "status %d %d %d %lf %lf\n",status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); if(status) break; // solver problem status = gsl_multiroot_test_residual(s->f, n*PREC); } while(status == GSL_CONTINUE && iter < MAXITER); - if (status == GSL_CONTINUE)// problem !!! - { - *check=1; - } - else// converged - { - // Msg(INFO, "status %d %d %d %lf %lf\n",status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); - convert_vector_from_gsl(s->x,x); - *check=0; - } + if (status == GSL_CONTINUE){ + *check=1; // problem !!! + } + else{ + // Msg(INFO, "status %d %d %d %lf %lf\n",status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1)); + convert_vector_from_gsl(s->x,x); + *check=0; // converged + } gsl_multiroot_fsolver_free(s); gsl_vector_free(xx);