Skip to content
Snippets Groups Projects
Commit 1f56bf3b authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

rename some variables to avoid confusion
parent fec8d9bc
No related branches found
No related tags found
No related merge requests found
// $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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment