findcurve.c

   1 /*
   2  */
   3
   4 #include <gsl/gsl_multiroots.h>
   5
   6 #include "symbolic.c"
   7
   8 #define X(i) gsl_vector_get(x,i)
   9 #define J_END_COL(i) \
  10   for (j=0; j<N; j++) gsl_matrix_set(f,i,j,J_COL[j]);
  11
  12 static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
  13   double *target = params;
  14   double F[6];
  15
  16   int i;
  17   X_EXTRACT;
  18   F_POPULATE;
  19   for (i=0; i<N; i++) gsl_vector_set(f,i, F(i) - target[i]);
  20   return 0;
  21 }
  22
  23 static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
  24   double J_COL[N];
  25   int j;
  26   X_EXTRACT;
  27   J_POPULATE;
  28 }
  29
  30 static int cb_fdf(const gsl_vector *x, void *params,
  31                   gsl_vector *f, gsl_matrix *J) {
  32   cb_f(x,params,f);
  33   cb_df(x,params,df);
  34 }
  35
  36 static double *scan1double(void) {
  37   double v;
  38   int r;
  39
  40   r = scanf("%lf",&v);
  41   if (ferror(stdin)) { perror("stdin"); exit(-1); }
  42   if (feof(stdin)) exit(0);
  43   if (r!=1) { fputs("bad input\n",stderr); exit(-1); }
  44   return v;
  45 }
  46
  47 static int main(void) {
  48   double target[N], epsilon;
  49
  50   gsl_multiroot_function_fdf function;
  51   function.f = cb_f;
  52   function.df = cb_df;
  53   function.fdf = cb_fdf;
  54   function.n = N;
  55   function.params = target;
  56
  57   gsl_vector *startpoint = gsl_vector_alloc(N);
  58
  59   gsl_multiroot_fdfsolver *solver
  60     = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybrid, N);
  61
  62   for (;;) {
  63     /* 2N+1 doubles: target, initial guess, epsilon for residual */
  64     for (i=0; i<N; i++)
  65       target[i] = scan1double();
  66     for (i=0; i<N; i++)
  67       gsl_vector_set(startpoint,i,scan1double());
  68     epsilon = scan1double();
  69
  70     gsl_multiroot_fdfsolver_set(solver, function, startpoint);
  71
  72     for (;;) {
  73       r = gsl_multiroot_fdfsolver_iterate(solver);
  74       if (r) break;
  75
  76       r = gsl_multiroot_test_residual(solver->f, epsilon);
  77       if (r != GSL_CONTINUE) break;
  78     }
  79
  80     if (r) {
  81       fprintf(stderr,"ERROR %s\n",gsl_strerror(gsl_errno));
  82       exit(-1);
  83     }
  84   }
  85 }