6 #include <gsl/gsl_errno.h>
7 #include <gsl/gsl_sf.h>
8 #include <gsl/gsl_multiroots.h>
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]);
16 static inline _Bool IS_SMALL(double v) {
20 static inline double sinc(double x) {
21 return gsl_sf_sinc(x / M_PI);
24 static int cb_f(const gsl_vector *x, void *params, gsl_vector *f) {
25 double *target = params;
31 for (i=0; i<N; i++) gsl_vector_set(f,i, F[i] - target[i]);
35 static int cb_df(const gsl_vector *x, void *params, gsl_matrix *J) {
43 static int cb_fdf(const gsl_vector *x, void *params,
44 gsl_vector *f, gsl_matrix *J) {
50 static double scan1double(void) {
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); }
62 double target[N], epsilon;
65 gsl_multiroot_function_fdf function;
68 function.fdf = cb_fdf;
70 function.params = target;
72 gsl_vector *startpoint = gsl_vector_alloc(N);
74 gsl_multiroot_fdfsolver *solver
75 = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_hybridsj, N);
78 /* 2N+1 doubles: target, initial guess, epsilon for residual */
80 target[i] = scan1double();
82 gsl_vector_set(startpoint,i,scan1double());
83 epsilon = scan1double();
85 gsl_multiroot_fdfsolver_set(solver, &function, startpoint);
89 r = gsl_multiroot_fdfsolver_iterate(solver);
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));
98 r = gsl_multiroot_test_residual(solver->f, epsilon);
99 if (r != GSL_CONTINUE) break;
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]);
109 fprintf(stderr," %10.5f,", gsl_matrix_get(j_dump, i, j));
110 fprintf(stderr,"\n");
112 fprintf(stderr,"ERROR %s\n",gsl_strerror(r));