findcurve.c

   1 /*
   2  */
   3
   4 #include <math.h>
   5
   6 #include <gsl/gsl_errno.h>
   7 #include <gsl/gsl_sf.h>
   8 #include <gsl/gsl_multiroots.h>
   9
  10 #include "symbolic.c"
  11
  12 #define X(i) gsl_vector_get(x,i)
  13 #define J_END_COL(i) \
  14   for (j=0; j<N; j++) gsl_matrix_set(J,i,j,J_COL[j]);
  15
  16 static inline _Bool IS_SMALL(double v) {
  17   return v < 1E6;
  18 }
  19
  20 static inline double sinc(double x) {
  21   return gsl_sf_sinc(x / M_PI);
  22 }
  23
  24 static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
  25   double *target = params;
  26   double F[6];
  27
  28   int i;
  29   X_EXTRACT;
  30   F_POPULATE;
  31   for (i=0; i<N; i++) gsl_vector_set(f,i, F[i] - target[i]);
  32   return 0;
  33 }
  34
  35 static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
  36   double J_COL[N];
  37   int j;
  38   X_EXTRACT;
  39   J_POPULATE;
  40   return 0;
  41 }
  42
  43 static int cb_fdf(const gsl_vector *x, void *params,
  44                   gsl_vector *f, gsl_matrix *J) {
  45   cb_f(x,params,f);
  46   cb_df(x,params,J);
  47   return 0;
  48 }
  49
  50 static double scan1double(void) {
  51   double v;
  52   int r;
  53
  54   r = scanf("%lf",&v);
  55   if (ferror(stdin)) { perror("stdin"); exit(-1); }
  56   if (feof(stdin)) exit(0);
  57   if (r!=1) { fputs("bad input\n",stderr); exit(-1); }
  58   return v;
  59 }
  60
  61 int main(void) {
  62   double target[N], epsilon;
  63   int i, j;
  64
  65   gsl_multiroot_function_fdf function;
  66   function.f = cb_f;
  67   function.df = cb_df;
  68   function.fdf = cb_fdf;
  69   function.n = N;
  70   function.params = target;
  71
  72   gsl_vector *startpoint = gsl_vector_alloc(N);
  73
  74   gsl_multiroot_fdfsolver *solver
  75     = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybridsj, N);
  76
  77   for (;;) {
  78     /* 2N+1 doubles: target, initial guess, epsilon for residual */
  79     for (i=0; i<N; i++)
  80       target[i] = scan1double();
  81     for (i=0; i<N; i++)
  82       gsl_vector_set(startpoint,i,scan1double());
  83     epsilon = scan1double();
  84
  85     gsl_multiroot_fdfsolver_set(solver, &function, startpoint);
  86
  87     int r;
  88     for (;;) {
  89       r = gsl_multiroot_fdfsolver_iterate(solver);
  90       if (r) break;
  91
  92       printf("[");
  93       for (i=0; i<6; i++) printf(" %.18g,", gsl_vector_get(solver->x, i));
  94       for (i=0; i<6; i++) printf(" %.18g,",
  95                                  target[i] + gsl_vector_get(solver->f, i));
  96       printf(" ]\n");
  97
  98       r = gsl_multiroot_test_residual(solver->f, epsilon);
  99       if (r != GSL_CONTINUE) break;
 100     }
 101
 102     if (r) {
 103       gsl_matrix *j_dump = gsl_matrix_alloc(N,N);
 104       cb_df(solver->x, target, j_dump);
 105       fprintf(stderr,"// FC J %35s %35s", "d(Q)/", "d(dQ)/\n");
 106       for (i=0; i<6; i++) {
 107         fprintf(stderr,"// FC J d./d%-6s ", PARAM_NAMES[i]);
 108         for (j=0; j<6; j++)
 109           fprintf(stderr," %10.5f,", gsl_matrix_get(j_dump, i, j));
 110         fprintf(stderr,"\n");
 111       }
 112       fprintf(stderr,"ERROR %s\n",gsl_strerror(r));
 113       exit(-1);
 114     }
 115
 116     printf("[]\n");
 117   }
 118 }