diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-08 08:48:12 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-08 08:48:12 +0200 |
commit | 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 (patch) | |
tree | fad79f909d9c3be53d68e6ebd67202650536d387 /packages/hmatrix/src/Numeric/GSL/gsl-aux.c | |
parent | eb3f702d065a4a967bb754977233e6eec408fd1f (diff) |
empty hmatrix-base
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/gsl-aux.c')
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/gsl-aux.c | 1541 |
1 files changed, 1541 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c new file mode 100644 index 0000000..410d157 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c | |||
@@ -0,0 +1,1541 @@ | |||
1 | #include <gsl/gsl_complex.h> | ||
2 | |||
3 | #define RVEC(A) int A##n, double*A##p | ||
4 | #define RMAT(A) int A##r, int A##c, double* A##p | ||
5 | #define KRVEC(A) int A##n, const double*A##p | ||
6 | #define KRMAT(A) int A##r, int A##c, const double* A##p | ||
7 | |||
8 | #define CVEC(A) int A##n, gsl_complex*A##p | ||
9 | #define CMAT(A) int A##r, int A##c, gsl_complex* A##p | ||
10 | #define KCVEC(A) int A##n, const gsl_complex*A##p | ||
11 | #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p | ||
12 | |||
13 | #define FVEC(A) int A##n, float*A##p | ||
14 | #define FMAT(A) int A##r, int A##c, float* A##p | ||
15 | #define KFVEC(A) int A##n, const float*A##p | ||
16 | #define KFMAT(A) int A##r, int A##c, const float* A##p | ||
17 | |||
18 | #define QVEC(A) int A##n, gsl_complex_float*A##p | ||
19 | #define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p | ||
20 | #define KQVEC(A) int A##n, const gsl_complex_float*A##p | ||
21 | #define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p | ||
22 | |||
23 | #include <gsl/gsl_blas.h> | ||
24 | #include <gsl/gsl_math.h> | ||
25 | #include <gsl/gsl_errno.h> | ||
26 | #include <gsl/gsl_fft_complex.h> | ||
27 | #include <gsl/gsl_integration.h> | ||
28 | #include <gsl/gsl_deriv.h> | ||
29 | #include <gsl/gsl_poly.h> | ||
30 | #include <gsl/gsl_multimin.h> | ||
31 | #include <gsl/gsl_multiroots.h> | ||
32 | #include <gsl/gsl_min.h> | ||
33 | #include <gsl/gsl_complex_math.h> | ||
34 | #include <gsl/gsl_rng.h> | ||
35 | #include <gsl/gsl_randist.h> | ||
36 | #include <gsl/gsl_roots.h> | ||
37 | #include <gsl/gsl_multifit_nlin.h> | ||
38 | #include <string.h> | ||
39 | #include <stdio.h> | ||
40 | |||
41 | #define MACRO(B) do {B} while (0) | ||
42 | #define ERROR(CODE) MACRO(return CODE;) | ||
43 | #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) | ||
44 | #define OK return 0; | ||
45 | |||
46 | #define MIN(A,B) ((A)<(B)?(A):(B)) | ||
47 | #define MAX(A,B) ((A)>(B)?(A):(B)) | ||
48 | |||
49 | #ifdef DBG | ||
50 | #define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); | ||
51 | #else | ||
52 | #define DEBUGMSG(M) | ||
53 | #endif | ||
54 | |||
55 | #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) | ||
56 | |||
57 | #ifdef DBG | ||
58 | #define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); | ||
59 | #else | ||
60 | #define DEBUGMAT(MSG,X) | ||
61 | #endif | ||
62 | |||
63 | #ifdef DBG | ||
64 | #define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); | ||
65 | #else | ||
66 | #define DEBUGVEC(MSG,X) | ||
67 | #endif | ||
68 | |||
69 | #define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) | ||
70 | #define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) | ||
71 | #define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) | ||
72 | #define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) | ||
73 | #define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) | ||
74 | #define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) | ||
75 | #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) | ||
76 | #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) | ||
77 | |||
78 | #define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) | ||
79 | #define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) | ||
80 | #define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) | ||
81 | #define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) | ||
82 | #define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) | ||
83 | #define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) | ||
84 | #define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) | ||
85 | #define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) | ||
86 | |||
87 | #define V(a) (&a.vector) | ||
88 | #define M(a) (&a.matrix) | ||
89 | |||
90 | #define GCVEC(A) int A##n, gsl_complex*A##p | ||
91 | #define KGCVEC(A) int A##n, const gsl_complex*A##p | ||
92 | |||
93 | #define GQVEC(A) int A##n, gsl_complex_float*A##p | ||
94 | #define KGQVEC(A) int A##n, const gsl_complex_float*A##p | ||
95 | |||
96 | #define BAD_SIZE 2000 | ||
97 | #define BAD_CODE 2001 | ||
98 | #define MEM 2002 | ||
99 | #define BAD_FILE 2003 | ||
100 | |||
101 | |||
102 | void no_abort_on_error() { | ||
103 | gsl_set_error_handler_off(); | ||
104 | } | ||
105 | |||
106 | |||
107 | int sumF(KFVEC(x),FVEC(r)) { | ||
108 | DEBUGMSG("sumF"); | ||
109 | REQUIRES(rn==1,BAD_SIZE); | ||
110 | int i; | ||
111 | float res = 0; | ||
112 | for (i = 0; i < xn; i++) res += xp[i]; | ||
113 | rp[0] = res; | ||
114 | OK | ||
115 | } | ||
116 | |||
117 | int sumR(KRVEC(x),RVEC(r)) { | ||
118 | DEBUGMSG("sumR"); | ||
119 | REQUIRES(rn==1,BAD_SIZE); | ||
120 | int i; | ||
121 | double res = 0; | ||
122 | for (i = 0; i < xn; i++) res += xp[i]; | ||
123 | rp[0] = res; | ||
124 | OK | ||
125 | } | ||
126 | |||
127 | int sumQ(KQVEC(x),QVEC(r)) { | ||
128 | DEBUGMSG("sumQ"); | ||
129 | REQUIRES(rn==1,BAD_SIZE); | ||
130 | int i; | ||
131 | gsl_complex_float res; | ||
132 | res.dat[0] = 0; | ||
133 | res.dat[1] = 0; | ||
134 | for (i = 0; i < xn; i++) { | ||
135 | res.dat[0] += xp[i].dat[0]; | ||
136 | res.dat[1] += xp[i].dat[1]; | ||
137 | } | ||
138 | rp[0] = res; | ||
139 | OK | ||
140 | } | ||
141 | |||
142 | int sumC(KCVEC(x),CVEC(r)) { | ||
143 | DEBUGMSG("sumC"); | ||
144 | REQUIRES(rn==1,BAD_SIZE); | ||
145 | int i; | ||
146 | gsl_complex res; | ||
147 | res.dat[0] = 0; | ||
148 | res.dat[1] = 0; | ||
149 | for (i = 0; i < xn; i++) { | ||
150 | res.dat[0] += xp[i].dat[0]; | ||
151 | res.dat[1] += xp[i].dat[1]; | ||
152 | } | ||
153 | rp[0] = res; | ||
154 | OK | ||
155 | } | ||
156 | |||
157 | int prodF(KFVEC(x),FVEC(r)) { | ||
158 | DEBUGMSG("prodF"); | ||
159 | REQUIRES(rn==1,BAD_SIZE); | ||
160 | int i; | ||
161 | float res = 1; | ||
162 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
163 | rp[0] = res; | ||
164 | OK | ||
165 | } | ||
166 | |||
167 | int prodR(KRVEC(x),RVEC(r)) { | ||
168 | DEBUGMSG("prodR"); | ||
169 | REQUIRES(rn==1,BAD_SIZE); | ||
170 | int i; | ||
171 | double res = 1; | ||
172 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
173 | rp[0] = res; | ||
174 | OK | ||
175 | } | ||
176 | |||
177 | int prodQ(KQVEC(x),QVEC(r)) { | ||
178 | DEBUGMSG("prodQ"); | ||
179 | REQUIRES(rn==1,BAD_SIZE); | ||
180 | int i; | ||
181 | gsl_complex_float res; | ||
182 | float temp; | ||
183 | res.dat[0] = 1; | ||
184 | res.dat[1] = 0; | ||
185 | for (i = 0; i < xn; i++) { | ||
186 | temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; | ||
187 | res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; | ||
188 | res.dat[0] = temp; | ||
189 | } | ||
190 | rp[0] = res; | ||
191 | OK | ||
192 | } | ||
193 | |||
194 | int prodC(KCVEC(x),CVEC(r)) { | ||
195 | DEBUGMSG("prodC"); | ||
196 | REQUIRES(rn==1,BAD_SIZE); | ||
197 | int i; | ||
198 | gsl_complex res; | ||
199 | double temp; | ||
200 | res.dat[0] = 1; | ||
201 | res.dat[1] = 0; | ||
202 | for (i = 0; i < xn; i++) { | ||
203 | temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; | ||
204 | res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; | ||
205 | res.dat[0] = temp; | ||
206 | } | ||
207 | rp[0] = res; | ||
208 | OK | ||
209 | } | ||
210 | |||
211 | int dotF(KFVEC(x), KFVEC(y), FVEC(r)) { | ||
212 | DEBUGMSG("dotF"); | ||
213 | REQUIRES(xn==yn,BAD_SIZE); | ||
214 | REQUIRES(rn==1,BAD_SIZE); | ||
215 | DEBUGMSG("dotF"); | ||
216 | KFVVIEW(x); | ||
217 | KFVVIEW(y); | ||
218 | gsl_blas_sdot(V(x),V(y),rp); | ||
219 | OK | ||
220 | } | ||
221 | |||
222 | int dotR(KRVEC(x), KRVEC(y), RVEC(r)) { | ||
223 | DEBUGMSG("dotR"); | ||
224 | REQUIRES(xn==yn,BAD_SIZE); | ||
225 | REQUIRES(rn==1,BAD_SIZE); | ||
226 | DEBUGMSG("dotR"); | ||
227 | KDVVIEW(x); | ||
228 | KDVVIEW(y); | ||
229 | gsl_blas_ddot(V(x),V(y),rp); | ||
230 | OK | ||
231 | } | ||
232 | |||
233 | int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) { | ||
234 | DEBUGMSG("dotQ"); | ||
235 | REQUIRES(xn==yn,BAD_SIZE); | ||
236 | REQUIRES(rn==1,BAD_SIZE); | ||
237 | DEBUGMSG("dotQ"); | ||
238 | KQVVIEW(x); | ||
239 | KQVVIEW(y); | ||
240 | gsl_blas_cdotu(V(x),V(y),rp); | ||
241 | OK | ||
242 | } | ||
243 | |||
244 | int dotC(KCVEC(x), KCVEC(y), CVEC(r)) { | ||
245 | DEBUGMSG("dotC"); | ||
246 | REQUIRES(xn==yn,BAD_SIZE); | ||
247 | REQUIRES(rn==1,BAD_SIZE); | ||
248 | DEBUGMSG("dotC"); | ||
249 | KCVVIEW(x); | ||
250 | KCVVIEW(y); | ||
251 | gsl_blas_zdotu(V(x),V(y),rp); | ||
252 | OK | ||
253 | } | ||
254 | |||
255 | int toScalarR(int code, KRVEC(x), RVEC(r)) { | ||
256 | REQUIRES(rn==1,BAD_SIZE); | ||
257 | DEBUGMSG("toScalarR"); | ||
258 | KDVVIEW(x); | ||
259 | double res; | ||
260 | switch(code) { | ||
261 | case 0: { res = gsl_blas_dnrm2(V(x)); break; } | ||
262 | case 1: { res = gsl_blas_dasum(V(x)); break; } | ||
263 | case 2: { res = gsl_vector_max_index(V(x)); break; } | ||
264 | case 3: { res = gsl_vector_max(V(x)); break; } | ||
265 | case 4: { res = gsl_vector_min_index(V(x)); break; } | ||
266 | case 5: { res = gsl_vector_min(V(x)); break; } | ||
267 | default: ERROR(BAD_CODE); | ||
268 | } | ||
269 | rp[0] = res; | ||
270 | OK | ||
271 | } | ||
272 | |||
273 | int toScalarF(int code, KFVEC(x), FVEC(r)) { | ||
274 | REQUIRES(rn==1,BAD_SIZE); | ||
275 | DEBUGMSG("toScalarF"); | ||
276 | KFVVIEW(x); | ||
277 | float res; | ||
278 | switch(code) { | ||
279 | case 0: { res = gsl_blas_snrm2(V(x)); break; } | ||
280 | case 1: { res = gsl_blas_sasum(V(x)); break; } | ||
281 | case 2: { res = gsl_vector_float_max_index(V(x)); break; } | ||
282 | case 3: { res = gsl_vector_float_max(V(x)); break; } | ||
283 | case 4: { res = gsl_vector_float_min_index(V(x)); break; } | ||
284 | case 5: { res = gsl_vector_float_min(V(x)); break; } | ||
285 | default: ERROR(BAD_CODE); | ||
286 | } | ||
287 | rp[0] = res; | ||
288 | OK | ||
289 | } | ||
290 | |||
291 | |||
292 | int toScalarC(int code, KCVEC(x), RVEC(r)) { | ||
293 | REQUIRES(rn==1,BAD_SIZE); | ||
294 | DEBUGMSG("toScalarC"); | ||
295 | KCVVIEW(x); | ||
296 | double res; | ||
297 | switch(code) { | ||
298 | case 0: { res = gsl_blas_dznrm2(V(x)); break; } | ||
299 | case 1: { res = gsl_blas_dzasum(V(x)); break; } | ||
300 | default: ERROR(BAD_CODE); | ||
301 | } | ||
302 | rp[0] = res; | ||
303 | OK | ||
304 | } | ||
305 | |||
306 | int toScalarQ(int code, KQVEC(x), FVEC(r)) { | ||
307 | REQUIRES(rn==1,BAD_SIZE); | ||
308 | DEBUGMSG("toScalarQ"); | ||
309 | KQVVIEW(x); | ||
310 | float res; | ||
311 | switch(code) { | ||
312 | case 0: { res = gsl_blas_scnrm2(V(x)); break; } | ||
313 | case 1: { res = gsl_blas_scasum(V(x)); break; } | ||
314 | default: ERROR(BAD_CODE); | ||
315 | } | ||
316 | rp[0] = res; | ||
317 | OK | ||
318 | } | ||
319 | |||
320 | |||
321 | inline double sign(double x) { | ||
322 | if(x>0) { | ||
323 | return +1.0; | ||
324 | } else if (x<0) { | ||
325 | return -1.0; | ||
326 | } else { | ||
327 | return 0.0; | ||
328 | } | ||
329 | } | ||
330 | |||
331 | inline float float_sign(float x) { | ||
332 | if(x>0) { | ||
333 | return +1.0; | ||
334 | } else if (x<0) { | ||
335 | return -1.0; | ||
336 | } else { | ||
337 | return 0.0; | ||
338 | } | ||
339 | } | ||
340 | |||
341 | inline gsl_complex complex_abs(gsl_complex z) { | ||
342 | gsl_complex r; | ||
343 | r.dat[0] = gsl_complex_abs(z); | ||
344 | r.dat[1] = 0; | ||
345 | return r; | ||
346 | } | ||
347 | |||
348 | inline gsl_complex complex_signum(gsl_complex z) { | ||
349 | gsl_complex r; | ||
350 | double mag; | ||
351 | if (z.dat[0] == 0 && z.dat[1] == 0) { | ||
352 | r.dat[0] = 0; | ||
353 | r.dat[1] = 0; | ||
354 | } else { | ||
355 | mag = gsl_complex_abs(z); | ||
356 | r.dat[0] = z.dat[0]/mag; | ||
357 | r.dat[1] = z.dat[1]/mag; | ||
358 | } | ||
359 | return r; | ||
360 | } | ||
361 | |||
362 | #define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK } | ||
363 | #define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK } | ||
364 | int mapR(int code, KRVEC(x), RVEC(r)) { | ||
365 | int k; | ||
366 | REQUIRES(xn == rn,BAD_SIZE); | ||
367 | DEBUGMSG("mapR"); | ||
368 | switch (code) { | ||
369 | OP(0,sin) | ||
370 | OP(1,cos) | ||
371 | OP(2,tan) | ||
372 | OP(3,fabs) | ||
373 | OP(4,asin) | ||
374 | OP(5,acos) | ||
375 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
376 | OP(7,sinh) | ||
377 | OP(8,cosh) | ||
378 | OP(9,tanh) | ||
379 | OP(10,gsl_asinh) | ||
380 | OP(11,gsl_acosh) | ||
381 | OP(12,gsl_atanh) | ||
382 | OP(13,exp) | ||
383 | OP(14,log) | ||
384 | OP(15,sign) | ||
385 | OP(16,sqrt) | ||
386 | default: ERROR(BAD_CODE); | ||
387 | } | ||
388 | } | ||
389 | |||
390 | int mapF(int code, KFVEC(x), FVEC(r)) { | ||
391 | int k; | ||
392 | REQUIRES(xn == rn,BAD_SIZE); | ||
393 | DEBUGMSG("mapF"); | ||
394 | switch (code) { | ||
395 | OP(0,sin) | ||
396 | OP(1,cos) | ||
397 | OP(2,tan) | ||
398 | OP(3,fabs) | ||
399 | OP(4,asin) | ||
400 | OP(5,acos) | ||
401 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
402 | OP(7,sinh) | ||
403 | OP(8,cosh) | ||
404 | OP(9,tanh) | ||
405 | OP(10,gsl_asinh) | ||
406 | OP(11,gsl_acosh) | ||
407 | OP(12,gsl_atanh) | ||
408 | OP(13,exp) | ||
409 | OP(14,log) | ||
410 | OP(15,sign) | ||
411 | OP(16,sqrt) | ||
412 | default: ERROR(BAD_CODE); | ||
413 | } | ||
414 | } | ||
415 | |||
416 | |||
417 | int mapCAux(int code, KGCVEC(x), GCVEC(r)) { | ||
418 | int k; | ||
419 | REQUIRES(xn == rn,BAD_SIZE); | ||
420 | DEBUGMSG("mapC"); | ||
421 | switch (code) { | ||
422 | OP(0,gsl_complex_sin) | ||
423 | OP(1,gsl_complex_cos) | ||
424 | OP(2,gsl_complex_tan) | ||
425 | OP(3,complex_abs) | ||
426 | OP(4,gsl_complex_arcsin) | ||
427 | OP(5,gsl_complex_arccos) | ||
428 | OP(6,gsl_complex_arctan) | ||
429 | OP(7,gsl_complex_sinh) | ||
430 | OP(8,gsl_complex_cosh) | ||
431 | OP(9,gsl_complex_tanh) | ||
432 | OP(10,gsl_complex_arcsinh) | ||
433 | OP(11,gsl_complex_arccosh) | ||
434 | OP(12,gsl_complex_arctanh) | ||
435 | OP(13,gsl_complex_exp) | ||
436 | OP(14,gsl_complex_log) | ||
437 | OP(15,complex_signum) | ||
438 | OP(16,gsl_complex_sqrt) | ||
439 | |||
440 | // gsl_complex_arg | ||
441 | // gsl_complex_abs | ||
442 | default: ERROR(BAD_CODE); | ||
443 | } | ||
444 | } | ||
445 | |||
446 | int mapC(int code, KCVEC(x), CVEC(r)) { | ||
447 | return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp); | ||
448 | } | ||
449 | |||
450 | |||
451 | gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a) | ||
452 | { | ||
453 | gsl_complex c; | ||
454 | gsl_complex r; | ||
455 | |||
456 | gsl_complex_float float_r; | ||
457 | |||
458 | c.dat[0] = a.dat[0]; | ||
459 | c.dat[1] = a.dat[1]; | ||
460 | |||
461 | r = (*cf)(c); | ||
462 | |||
463 | float_r.dat[0] = r.dat[0]; | ||
464 | float_r.dat[1] = r.dat[1]; | ||
465 | |||
466 | return float_r; | ||
467 | } | ||
468 | |||
469 | gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), | ||
470 | gsl_complex_float a,gsl_complex_float b) | ||
471 | { | ||
472 | gsl_complex c1; | ||
473 | gsl_complex c2; | ||
474 | gsl_complex r; | ||
475 | |||
476 | gsl_complex_float float_r; | ||
477 | |||
478 | c1.dat[0] = a.dat[0]; | ||
479 | c1.dat[1] = a.dat[1]; | ||
480 | |||
481 | c2.dat[0] = b.dat[0]; | ||
482 | c2.dat[1] = b.dat[1]; | ||
483 | |||
484 | r = (*cf)(c1,c2); | ||
485 | |||
486 | float_r.dat[0] = r.dat[0]; | ||
487 | float_r.dat[1] = r.dat[1]; | ||
488 | |||
489 | return float_r; | ||
490 | } | ||
491 | |||
492 | #define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK } | ||
493 | #define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK } | ||
494 | int mapQAux(int code, KGQVEC(x), GQVEC(r)) { | ||
495 | int k; | ||
496 | REQUIRES(xn == rn,BAD_SIZE); | ||
497 | DEBUGMSG("mapQ"); | ||
498 | switch (code) { | ||
499 | OPC(0,gsl_complex_sin) | ||
500 | OPC(1,gsl_complex_cos) | ||
501 | OPC(2,gsl_complex_tan) | ||
502 | OPC(3,complex_abs) | ||
503 | OPC(4,gsl_complex_arcsin) | ||
504 | OPC(5,gsl_complex_arccos) | ||
505 | OPC(6,gsl_complex_arctan) | ||
506 | OPC(7,gsl_complex_sinh) | ||
507 | OPC(8,gsl_complex_cosh) | ||
508 | OPC(9,gsl_complex_tanh) | ||
509 | OPC(10,gsl_complex_arcsinh) | ||
510 | OPC(11,gsl_complex_arccosh) | ||
511 | OPC(12,gsl_complex_arctanh) | ||
512 | OPC(13,gsl_complex_exp) | ||
513 | OPC(14,gsl_complex_log) | ||
514 | OPC(15,complex_signum) | ||
515 | OPC(16,gsl_complex_sqrt) | ||
516 | |||
517 | // gsl_complex_arg | ||
518 | // gsl_complex_abs | ||
519 | default: ERROR(BAD_CODE); | ||
520 | } | ||
521 | } | ||
522 | |||
523 | int mapQ(int code, KQVEC(x), QVEC(r)) { | ||
524 | return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp); | ||
525 | } | ||
526 | |||
527 | |||
528 | int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { | ||
529 | int k; | ||
530 | double val = *pval; | ||
531 | REQUIRES(xn == rn,BAD_SIZE); | ||
532 | DEBUGMSG("mapValR"); | ||
533 | switch (code) { | ||
534 | OPV(0,val*xp[k]) | ||
535 | OPV(1,val/xp[k]) | ||
536 | OPV(2,val+xp[k]) | ||
537 | OPV(3,val-xp[k]) | ||
538 | OPV(4,pow(val,xp[k])) | ||
539 | OPV(5,pow(xp[k],val)) | ||
540 | default: ERROR(BAD_CODE); | ||
541 | } | ||
542 | } | ||
543 | |||
544 | int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) { | ||
545 | int k; | ||
546 | float val = *pval; | ||
547 | REQUIRES(xn == rn,BAD_SIZE); | ||
548 | DEBUGMSG("mapValF"); | ||
549 | switch (code) { | ||
550 | OPV(0,val*xp[k]) | ||
551 | OPV(1,val/xp[k]) | ||
552 | OPV(2,val+xp[k]) | ||
553 | OPV(3,val-xp[k]) | ||
554 | OPV(4,pow(val,xp[k])) | ||
555 | OPV(5,pow(xp[k],val)) | ||
556 | default: ERROR(BAD_CODE); | ||
557 | } | ||
558 | } | ||
559 | |||
560 | int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) { | ||
561 | int k; | ||
562 | gsl_complex val = *pval; | ||
563 | REQUIRES(xn == rn,BAD_SIZE); | ||
564 | DEBUGMSG("mapValC"); | ||
565 | switch (code) { | ||
566 | OPV(0,gsl_complex_mul(val,xp[k])) | ||
567 | OPV(1,gsl_complex_div(val,xp[k])) | ||
568 | OPV(2,gsl_complex_add(val,xp[k])) | ||
569 | OPV(3,gsl_complex_sub(val,xp[k])) | ||
570 | OPV(4,gsl_complex_pow(val,xp[k])) | ||
571 | OPV(5,gsl_complex_pow(xp[k],val)) | ||
572 | default: ERROR(BAD_CODE); | ||
573 | } | ||
574 | } | ||
575 | |||
576 | int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) { | ||
577 | return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp); | ||
578 | } | ||
579 | |||
580 | |||
581 | int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) { | ||
582 | int k; | ||
583 | gsl_complex_float val = *pval; | ||
584 | REQUIRES(xn == rn,BAD_SIZE); | ||
585 | DEBUGMSG("mapValQ"); | ||
586 | switch (code) { | ||
587 | OPCA(0,gsl_complex_mul,val,xp[k]) | ||
588 | OPCA(1,gsl_complex_div,val,xp[k]) | ||
589 | OPCA(2,gsl_complex_add,val,xp[k]) | ||
590 | OPCA(3,gsl_complex_sub,val,xp[k]) | ||
591 | OPCA(4,gsl_complex_pow,val,xp[k]) | ||
592 | OPCA(5,gsl_complex_pow,xp[k],val) | ||
593 | default: ERROR(BAD_CODE); | ||
594 | } | ||
595 | } | ||
596 | |||
597 | int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) { | ||
598 | return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp); | ||
599 | } | ||
600 | |||
601 | |||
602 | #define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK } | ||
603 | #define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK } | ||
604 | int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) { | ||
605 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
606 | int k; | ||
607 | switch(code) { | ||
608 | OPZE(4,"zipR Pow",pow) | ||
609 | OPZE(5,"zipR ATan2",atan2) | ||
610 | } | ||
611 | KDVVIEW(a); | ||
612 | KDVVIEW(b); | ||
613 | DVVIEW(r); | ||
614 | gsl_vector_memcpy(V(r),V(a)); | ||
615 | int res; | ||
616 | switch(code) { | ||
617 | OPZV(0,"zipR Add",gsl_vector_add) | ||
618 | OPZV(1,"zipR Sub",gsl_vector_sub) | ||
619 | OPZV(2,"zipR Mul",gsl_vector_mul) | ||
620 | OPZV(3,"zipR Div",gsl_vector_div) | ||
621 | default: ERROR(BAD_CODE); | ||
622 | } | ||
623 | } | ||
624 | |||
625 | |||
626 | int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) { | ||
627 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
628 | int k; | ||
629 | switch(code) { | ||
630 | OPZE(4,"zipF Pow",pow) | ||
631 | OPZE(5,"zipF ATan2",atan2) | ||
632 | } | ||
633 | KFVVIEW(a); | ||
634 | KFVVIEW(b); | ||
635 | FVVIEW(r); | ||
636 | gsl_vector_float_memcpy(V(r),V(a)); | ||
637 | int res; | ||
638 | switch(code) { | ||
639 | OPZV(0,"zipF Add",gsl_vector_float_add) | ||
640 | OPZV(1,"zipF Sub",gsl_vector_float_sub) | ||
641 | OPZV(2,"zipF Mul",gsl_vector_float_mul) | ||
642 | OPZV(3,"zipF Div",gsl_vector_float_div) | ||
643 | default: ERROR(BAD_CODE); | ||
644 | } | ||
645 | } | ||
646 | |||
647 | |||
648 | int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) { | ||
649 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
650 | int k; | ||
651 | switch(code) { | ||
652 | OPZE(0,"zipC Add",gsl_complex_add) | ||
653 | OPZE(1,"zipC Sub",gsl_complex_sub) | ||
654 | OPZE(2,"zipC Mul",gsl_complex_mul) | ||
655 | OPZE(3,"zipC Div",gsl_complex_div) | ||
656 | OPZE(4,"zipC Pow",gsl_complex_pow) | ||
657 | //OPZE(5,"zipR ATan2",atan2) | ||
658 | } | ||
659 | //KCVVIEW(a); | ||
660 | //KCVVIEW(b); | ||
661 | //CVVIEW(r); | ||
662 | //gsl_vector_memcpy(V(r),V(a)); | ||
663 | //int res; | ||
664 | switch(code) { | ||
665 | default: ERROR(BAD_CODE); | ||
666 | } | ||
667 | } | ||
668 | |||
669 | |||
670 | int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { | ||
671 | return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp); | ||
672 | } | ||
673 | |||
674 | |||
675 | #define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK } | ||
676 | int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) { | ||
677 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
678 | int k; | ||
679 | switch(code) { | ||
680 | OPCZE(0,"zipQ Add",gsl_complex_add) | ||
681 | OPCZE(1,"zipQ Sub",gsl_complex_sub) | ||
682 | OPCZE(2,"zipQ Mul",gsl_complex_mul) | ||
683 | OPCZE(3,"zipQ Div",gsl_complex_div) | ||
684 | OPCZE(4,"zipQ Pow",gsl_complex_pow) | ||
685 | //OPZE(5,"zipR ATan2",atan2) | ||
686 | } | ||
687 | //KCVVIEW(a); | ||
688 | //KCVVIEW(b); | ||
689 | //CVVIEW(r); | ||
690 | //gsl_vector_memcpy(V(r),V(a)); | ||
691 | //int res; | ||
692 | switch(code) { | ||
693 | default: ERROR(BAD_CODE); | ||
694 | } | ||
695 | } | ||
696 | |||
697 | |||
698 | int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) { | ||
699 | return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp); | ||
700 | } | ||
701 | |||
702 | |||
703 | |||
704 | int fft(int code, KCVEC(X), CVEC(R)) { | ||
705 | REQUIRES(Xn == Rn,BAD_SIZE); | ||
706 | DEBUGMSG("fft"); | ||
707 | int s = Xn; | ||
708 | gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s); | ||
709 | gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s); | ||
710 | gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn); | ||
711 | gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn); | ||
712 | gsl_blas_dcopy(&X.vector,&R.vector); | ||
713 | if(code==0) { | ||
714 | gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace); | ||
715 | } else { | ||
716 | gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace); | ||
717 | } | ||
718 | gsl_fft_complex_wavetable_free (wavetable); | ||
719 | gsl_fft_complex_workspace_free (workspace); | ||
720 | OK | ||
721 | } | ||
722 | |||
723 | |||
724 | int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr) | ||
725 | { | ||
726 | gsl_function F; | ||
727 | F.function = f; | ||
728 | F.params = 0; | ||
729 | |||
730 | if(code==0) return gsl_deriv_central (&F, x, h, result, abserr); | ||
731 | |||
732 | if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr); | ||
733 | |||
734 | if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr); | ||
735 | |||
736 | return 0; | ||
737 | } | ||
738 | |||
739 | |||
740 | int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec, | ||
741 | double *result, double*error) { | ||
742 | DEBUGMSG("integrate_qng"); | ||
743 | gsl_function F; | ||
744 | F.function = f; | ||
745 | F.params = NULL; | ||
746 | size_t neval; | ||
747 | int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval); | ||
748 | CHECK(res,res); | ||
749 | OK | ||
750 | } | ||
751 | |||
752 | int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w, | ||
753 | double *result, double* error) { | ||
754 | DEBUGMSG("integrate_qags"); | ||
755 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
756 | gsl_function F; | ||
757 | F.function = f; | ||
758 | F.params = NULL; | ||
759 | int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error); | ||
760 | CHECK(res,res); | ||
761 | gsl_integration_workspace_free (wk); | ||
762 | OK | ||
763 | } | ||
764 | |||
765 | int integrate_qagi(double f(double,void*), double aprec, double prec, int w, | ||
766 | double *result, double* error) { | ||
767 | DEBUGMSG("integrate_qagi"); | ||
768 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
769 | gsl_function F; | ||
770 | F.function = f; | ||
771 | F.params = NULL; | ||
772 | int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error); | ||
773 | CHECK(res,res); | ||
774 | gsl_integration_workspace_free (wk); | ||
775 | OK | ||
776 | } | ||
777 | |||
778 | |||
779 | int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w, | ||
780 | double *result, double* error) { | ||
781 | DEBUGMSG("integrate_qagiu"); | ||
782 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
783 | gsl_function F; | ||
784 | F.function = f; | ||
785 | F.params = NULL; | ||
786 | int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error); | ||
787 | CHECK(res,res); | ||
788 | gsl_integration_workspace_free (wk); | ||
789 | OK | ||
790 | } | ||
791 | |||
792 | |||
793 | int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w, | ||
794 | double *result, double* error) { | ||
795 | DEBUGMSG("integrate_qagil"); | ||
796 | gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w); | ||
797 | gsl_function F; | ||
798 | F.function = f; | ||
799 | F.params = NULL; | ||
800 | int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error); | ||
801 | CHECK(res,res); | ||
802 | gsl_integration_workspace_free (wk); | ||
803 | OK | ||
804 | } | ||
805 | |||
806 | int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec, | ||
807 | int w, double *result, double* error, int *neval) { | ||
808 | DEBUGMSG("integrate_cquad"); | ||
809 | gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w); | ||
810 | gsl_function F; | ||
811 | F.function = f; | ||
812 | F.params = NULL; | ||
813 | size_t * sneval = NULL; | ||
814 | int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval); | ||
815 | *neval = *sneval; | ||
816 | CHECK(res,res); | ||
817 | gsl_integration_cquad_workspace_free (wk); | ||
818 | OK | ||
819 | } | ||
820 | |||
821 | |||
822 | int polySolve(KRVEC(a), CVEC(z)) { | ||
823 | DEBUGMSG("polySolve"); | ||
824 | REQUIRES(an>1,BAD_SIZE); | ||
825 | gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an); | ||
826 | int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp); | ||
827 | CHECK(res,res); | ||
828 | gsl_poly_complex_workspace_free (w); | ||
829 | OK; | ||
830 | } | ||
831 | |||
832 | int vector_fscanf(char*filename, RVEC(a)) { | ||
833 | DEBUGMSG("gsl_vector_fscanf"); | ||
834 | DVVIEW(a); | ||
835 | FILE * f = fopen(filename,"r"); | ||
836 | CHECK(!f,BAD_FILE); | ||
837 | int res = gsl_vector_fscanf(f,V(a)); | ||
838 | CHECK(res,res); | ||
839 | fclose (f); | ||
840 | OK | ||
841 | } | ||
842 | |||
843 | int vector_fprintf(char*filename, char*fmt, RVEC(a)) { | ||
844 | DEBUGMSG("gsl_vector_fprintf"); | ||
845 | DVVIEW(a); | ||
846 | FILE * f = fopen(filename,"w"); | ||
847 | CHECK(!f,BAD_FILE); | ||
848 | int res = gsl_vector_fprintf(f,V(a),fmt); | ||
849 | CHECK(res,res); | ||
850 | fclose (f); | ||
851 | OK | ||
852 | } | ||
853 | |||
854 | int vector_fread(char*filename, RVEC(a)) { | ||
855 | DEBUGMSG("gsl_vector_fread"); | ||
856 | DVVIEW(a); | ||
857 | FILE * f = fopen(filename,"r"); | ||
858 | CHECK(!f,BAD_FILE); | ||
859 | int res = gsl_vector_fread(f,V(a)); | ||
860 | CHECK(res,res); | ||
861 | fclose (f); | ||
862 | OK | ||
863 | } | ||
864 | |||
865 | int vector_fwrite(char*filename, RVEC(a)) { | ||
866 | DEBUGMSG("gsl_vector_fwrite"); | ||
867 | DVVIEW(a); | ||
868 | FILE * f = fopen(filename,"w"); | ||
869 | CHECK(!f,BAD_FILE); | ||
870 | int res = gsl_vector_fwrite(f,V(a)); | ||
871 | CHECK(res,res); | ||
872 | fclose (f); | ||
873 | OK | ||
874 | } | ||
875 | |||
876 | int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) { | ||
877 | DEBUGMSG("matrix_fprintf"); | ||
878 | FILE * f = fopen(filename,"w"); | ||
879 | CHECK(!f,BAD_FILE); | ||
880 | int i,j,sr,sc; | ||
881 | if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} | ||
882 | #define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) | ||
883 | for (i=0; i<mr; i++) { | ||
884 | for (j=0; j<mc-1; j++) { | ||
885 | fprintf(f,fmt,AT(m,i,j)); | ||
886 | fprintf(f," "); | ||
887 | } | ||
888 | fprintf(f,fmt,AT(m,i,j)); | ||
889 | fprintf(f,"\n"); | ||
890 | } | ||
891 | fclose (f); | ||
892 | OK | ||
893 | } | ||
894 | |||
895 | //--------------------------------------------------------------- | ||
896 | |||
897 | typedef double Trawfun(int, double*); | ||
898 | |||
899 | double only_f_aux_min(const gsl_vector*x, void *pars) { | ||
900 | Trawfun * f = (Trawfun*) pars; | ||
901 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
902 | int k; | ||
903 | for(k=0;k<x->size;k++) { | ||
904 | p[k] = gsl_vector_get(x,k); | ||
905 | } | ||
906 | double res = f(x->size,p); | ||
907 | free(p); | ||
908 | return res; | ||
909 | } | ||
910 | |||
911 | double only_f_aux_root(double x, void *pars); | ||
912 | int uniMinimize(int method, double f(double), | ||
913 | double epsrel, int maxit, double min, | ||
914 | double xl, double xu, RMAT(sol)) { | ||
915 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
916 | DEBUGMSG("minimize_only_f"); | ||
917 | gsl_function my_func; | ||
918 | my_func.function = only_f_aux_root; | ||
919 | my_func.params = f; | ||
920 | size_t iter = 0; | ||
921 | int status; | ||
922 | const gsl_min_fminimizer_type *T; | ||
923 | gsl_min_fminimizer *s; | ||
924 | // Starting point | ||
925 | switch(method) { | ||
926 | case 0 : {T = gsl_min_fminimizer_goldensection; break; } | ||
927 | case 1 : {T = gsl_min_fminimizer_brent; break; } | ||
928 | case 2 : {T = gsl_min_fminimizer_quad_golden; break; } | ||
929 | default: ERROR(BAD_CODE); | ||
930 | } | ||
931 | s = gsl_min_fminimizer_alloc (T); | ||
932 | gsl_min_fminimizer_set (s, &my_func, min, xl, xu); | ||
933 | do { | ||
934 | double current_min, current_lo, current_hi; | ||
935 | status = gsl_min_fminimizer_iterate (s); | ||
936 | current_min = gsl_min_fminimizer_x_minimum (s); | ||
937 | current_lo = gsl_min_fminimizer_x_lower (s); | ||
938 | current_hi = gsl_min_fminimizer_x_upper (s); | ||
939 | solp[iter*solc] = iter + 1; | ||
940 | solp[iter*solc+1] = current_min; | ||
941 | solp[iter*solc+2] = current_lo; | ||
942 | solp[iter*solc+3] = current_hi; | ||
943 | iter++; | ||
944 | if (status) /* check if solver is stuck */ | ||
945 | break; | ||
946 | |||
947 | status = | ||
948 | gsl_min_test_interval (current_lo, current_hi, 0, epsrel); | ||
949 | } | ||
950 | while (status == GSL_CONTINUE && iter < maxit); | ||
951 | int i; | ||
952 | for (i=iter; i<solr; i++) { | ||
953 | solp[i*solc+0] = iter; | ||
954 | solp[i*solc+1]=0.; | ||
955 | solp[i*solc+2]=0.; | ||
956 | solp[i*solc+3]=0.; | ||
957 | } | ||
958 | gsl_min_fminimizer_free(s); | ||
959 | OK | ||
960 | } | ||
961 | |||
962 | |||
963 | |||
964 | // this version returns info about intermediate steps | ||
965 | int minimize(int method, double f(int, double*), double tolsize, int maxit, | ||
966 | KRVEC(xi), KRVEC(sz), RMAT(sol)) { | ||
967 | REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE); | ||
968 | DEBUGMSG("minimizeList (nmsimplex)"); | ||
969 | gsl_multimin_function my_func; | ||
970 | // extract function from pars | ||
971 | my_func.f = only_f_aux_min; | ||
972 | my_func.n = xin; | ||
973 | my_func.params = f; | ||
974 | size_t iter = 0; | ||
975 | int status; | ||
976 | double size; | ||
977 | const gsl_multimin_fminimizer_type *T; | ||
978 | gsl_multimin_fminimizer *s = NULL; | ||
979 | // Initial vertex size vector | ||
980 | KDVVIEW(sz); | ||
981 | // Starting point | ||
982 | KDVVIEW(xi); | ||
983 | // Minimizer nmsimplex, without derivatives | ||
984 | switch(method) { | ||
985 | case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; } | ||
986 | #ifdef GSL110 | ||
987 | case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; } | ||
988 | #else | ||
989 | case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; } | ||
990 | #endif | ||
991 | default: ERROR(BAD_CODE); | ||
992 | } | ||
993 | s = gsl_multimin_fminimizer_alloc (T, my_func.n); | ||
994 | gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz)); | ||
995 | do { | ||
996 | status = gsl_multimin_fminimizer_iterate (s); | ||
997 | size = gsl_multimin_fminimizer_size (s); | ||
998 | |||
999 | solp[iter*solc+0] = iter+1; | ||
1000 | solp[iter*solc+1] = s->fval; | ||
1001 | solp[iter*solc+2] = size; | ||
1002 | |||
1003 | int k; | ||
1004 | for(k=0;k<xin;k++) { | ||
1005 | solp[iter*solc+k+3] = gsl_vector_get(s->x,k); | ||
1006 | } | ||
1007 | iter++; | ||
1008 | if (status) break; | ||
1009 | status = gsl_multimin_test_size (size, tolsize); | ||
1010 | } while (status == GSL_CONTINUE && iter < maxit); | ||
1011 | int i,j; | ||
1012 | for (i=iter; i<solr; i++) { | ||
1013 | solp[i*solc+0] = iter; | ||
1014 | for(j=1;j<solc;j++) { | ||
1015 | solp[i*solc+j]=0.; | ||
1016 | } | ||
1017 | } | ||
1018 | gsl_multimin_fminimizer_free(s); | ||
1019 | OK | ||
1020 | } | ||
1021 | |||
1022 | // working with the gradient | ||
1023 | |||
1024 | typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf; | ||
1025 | |||
1026 | double f_aux_min(const gsl_vector*x, void *pars) { | ||
1027 | Tfdf * fdf = ((Tfdf*) pars); | ||
1028 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
1029 | int k; | ||
1030 | for(k=0;k<x->size;k++) { | ||
1031 | p[k] = gsl_vector_get(x,k); | ||
1032 | } | ||
1033 | double res = fdf->f(x->size,p); | ||
1034 | free(p); | ||
1035 | return res; | ||
1036 | } | ||
1037 | |||
1038 | |||
1039 | void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { | ||
1040 | Tfdf * fdf = ((Tfdf*) pars); | ||
1041 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
1042 | double* q = (double*)calloc(g->size,sizeof(double)); | ||
1043 | int k; | ||
1044 | for(k=0;k<x->size;k++) { | ||
1045 | p[k] = gsl_vector_get(x,k); | ||
1046 | } | ||
1047 | |||
1048 | fdf->df(x->size,p,g->size,q); | ||
1049 | |||
1050 | for(k=0;k<x->size;k++) { | ||
1051 | gsl_vector_set(g,k,q[k]); | ||
1052 | } | ||
1053 | free(p); | ||
1054 | free(q); | ||
1055 | } | ||
1056 | |||
1057 | void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) { | ||
1058 | *f = f_aux_min(x,pars); | ||
1059 | df_aux_min(x,pars,g); | ||
1060 | } | ||
1061 | |||
1062 | |||
1063 | int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*), | ||
1064 | double initstep, double minimpar, double tolgrad, int maxit, | ||
1065 | KRVEC(xi), RMAT(sol)) { | ||
1066 | REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); | ||
1067 | DEBUGMSG("minimizeWithDeriv (conjugate_fr)"); | ||
1068 | gsl_multimin_function_fdf my_func; | ||
1069 | // extract function from pars | ||
1070 | my_func.f = f_aux_min; | ||
1071 | my_func.df = df_aux_min; | ||
1072 | my_func.fdf = fdf_aux_min; | ||
1073 | my_func.n = xin; | ||
1074 | Tfdf stfdf; | ||
1075 | stfdf.f = f; | ||
1076 | stfdf.df = df; | ||
1077 | my_func.params = &stfdf; | ||
1078 | size_t iter = 0; | ||
1079 | int status; | ||
1080 | const gsl_multimin_fdfminimizer_type *T; | ||
1081 | gsl_multimin_fdfminimizer *s = NULL; | ||
1082 | // Starting point | ||
1083 | KDVVIEW(xi); | ||
1084 | // conjugate gradient fr | ||
1085 | switch(method) { | ||
1086 | case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; } | ||
1087 | case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; } | ||
1088 | case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; } | ||
1089 | case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; } | ||
1090 | case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; } | ||
1091 | default: ERROR(BAD_CODE); | ||
1092 | } | ||
1093 | s = gsl_multimin_fdfminimizer_alloc (T, my_func.n); | ||
1094 | gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar); | ||
1095 | do { | ||
1096 | status = gsl_multimin_fdfminimizer_iterate (s); | ||
1097 | solp[iter*solc+0] = iter+1; | ||
1098 | solp[iter*solc+1] = s->f; | ||
1099 | int k; | ||
1100 | for(k=0;k<xin;k++) { | ||
1101 | solp[iter*solc+k+2] = gsl_vector_get(s->x,k); | ||
1102 | } | ||
1103 | iter++; | ||
1104 | if (status) break; | ||
1105 | status = gsl_multimin_test_gradient (s->gradient, tolgrad); | ||
1106 | } while (status == GSL_CONTINUE && iter < maxit); | ||
1107 | int i,j; | ||
1108 | for (i=iter; i<solr; i++) { | ||
1109 | solp[i*solc+0] = iter; | ||
1110 | for(j=1;j<solc;j++) { | ||
1111 | solp[i*solc+j]=0.; | ||
1112 | } | ||
1113 | } | ||
1114 | gsl_multimin_fdfminimizer_free(s); | ||
1115 | OK | ||
1116 | } | ||
1117 | |||
1118 | //--------------------------------------------------------------- | ||
1119 | |||
1120 | double only_f_aux_root(double x, void *pars) { | ||
1121 | double (*f)(double) = (double (*)(double)) pars; | ||
1122 | return f(x); | ||
1123 | } | ||
1124 | |||
1125 | int root(int method, double f(double), | ||
1126 | double epsrel, int maxit, | ||
1127 | double xl, double xu, RMAT(sol)) { | ||
1128 | REQUIRES(solr == maxit && solc == 4,BAD_SIZE); | ||
1129 | DEBUGMSG("root_only_f"); | ||
1130 | gsl_function my_func; | ||
1131 | // extract function from pars | ||
1132 | my_func.function = only_f_aux_root; | ||
1133 | my_func.params = f; | ||
1134 | size_t iter = 0; | ||
1135 | int status; | ||
1136 | const gsl_root_fsolver_type *T; | ||
1137 | gsl_root_fsolver *s; | ||
1138 | // Starting point | ||
1139 | switch(method) { | ||
1140 | case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; } | ||
1141 | case 1 : {T = gsl_root_fsolver_falsepos; break; } | ||
1142 | case 2 : {T = gsl_root_fsolver_brent; break; } | ||
1143 | default: ERROR(BAD_CODE); | ||
1144 | } | ||
1145 | s = gsl_root_fsolver_alloc (T); | ||
1146 | gsl_root_fsolver_set (s, &my_func, xl, xu); | ||
1147 | do { | ||
1148 | double best, current_lo, current_hi; | ||
1149 | status = gsl_root_fsolver_iterate (s); | ||
1150 | best = gsl_root_fsolver_root (s); | ||
1151 | current_lo = gsl_root_fsolver_x_lower (s); | ||
1152 | current_hi = gsl_root_fsolver_x_upper (s); | ||
1153 | solp[iter*solc] = iter + 1; | ||
1154 | solp[iter*solc+1] = best; | ||
1155 | solp[iter*solc+2] = current_lo; | ||
1156 | solp[iter*solc+3] = current_hi; | ||
1157 | iter++; | ||
1158 | if (status) /* check if solver is stuck */ | ||
1159 | break; | ||
1160 | |||
1161 | status = | ||
1162 | gsl_root_test_interval (current_lo, current_hi, 0, epsrel); | ||
1163 | } | ||
1164 | while (status == GSL_CONTINUE && iter < maxit); | ||
1165 | int i; | ||
1166 | for (i=iter; i<solr; i++) { | ||
1167 | solp[i*solc+0] = iter; | ||
1168 | solp[i*solc+1]=0.; | ||
1169 | solp[i*solc+2]=0.; | ||
1170 | solp[i*solc+3]=0.; | ||
1171 | } | ||
1172 | gsl_root_fsolver_free(s); | ||
1173 | OK | ||
1174 | } | ||
1175 | |||
1176 | typedef struct { | ||
1177 | double (*f)(double); | ||
1178 | double (*jf)(double); | ||
1179 | } uniTfjf; | ||
1180 | |||
1181 | double f_aux_uni(double x, void *pars) { | ||
1182 | uniTfjf * fjf = ((uniTfjf*) pars); | ||
1183 | return (fjf->f)(x); | ||
1184 | } | ||
1185 | |||
1186 | double jf_aux_uni(double x, void * pars) { | ||
1187 | uniTfjf * fjf = ((uniTfjf*) pars); | ||
1188 | return (fjf->jf)(x); | ||
1189 | } | ||
1190 | |||
1191 | void fjf_aux_uni(double x, void * pars, double * f, double * g) { | ||
1192 | *f = f_aux_uni(x,pars); | ||
1193 | *g = jf_aux_uni(x,pars); | ||
1194 | } | ||
1195 | |||
1196 | int rootj(int method, double f(double), | ||
1197 | double df(double), | ||
1198 | double epsrel, int maxit, | ||
1199 | double x, RMAT(sol)) { | ||
1200 | REQUIRES(solr == maxit && solc == 2,BAD_SIZE); | ||
1201 | DEBUGMSG("root_fjf"); | ||
1202 | gsl_function_fdf my_func; | ||
1203 | // extract function from pars | ||
1204 | my_func.f = f_aux_uni; | ||
1205 | my_func.df = jf_aux_uni; | ||
1206 | my_func.fdf = fjf_aux_uni; | ||
1207 | uniTfjf stfjf; | ||
1208 | stfjf.f = f; | ||
1209 | stfjf.jf = df; | ||
1210 | my_func.params = &stfjf; | ||
1211 | size_t iter = 0; | ||
1212 | int status; | ||
1213 | const gsl_root_fdfsolver_type *T; | ||
1214 | gsl_root_fdfsolver *s; | ||
1215 | // Starting point | ||
1216 | switch(method) { | ||
1217 | case 0 : {T = gsl_root_fdfsolver_newton;; break; } | ||
1218 | case 1 : {T = gsl_root_fdfsolver_secant; break; } | ||
1219 | case 2 : {T = gsl_root_fdfsolver_steffenson; break; } | ||
1220 | default: ERROR(BAD_CODE); | ||
1221 | } | ||
1222 | s = gsl_root_fdfsolver_alloc (T); | ||
1223 | |||
1224 | gsl_root_fdfsolver_set (s, &my_func, x); | ||
1225 | |||
1226 | do { | ||
1227 | double x0; | ||
1228 | status = gsl_root_fdfsolver_iterate (s); | ||
1229 | x0 = x; | ||
1230 | x = gsl_root_fdfsolver_root(s); | ||
1231 | solp[iter*solc+0] = iter+1; | ||
1232 | solp[iter*solc+1] = x; | ||
1233 | |||
1234 | iter++; | ||
1235 | if (status) /* check if solver is stuck */ | ||
1236 | break; | ||
1237 | |||
1238 | status = | ||
1239 | gsl_root_test_delta (x, x0, 0, epsrel); | ||
1240 | } | ||
1241 | while (status == GSL_CONTINUE && iter < maxit); | ||
1242 | |||
1243 | int i; | ||
1244 | for (i=iter; i<solr; i++) { | ||
1245 | solp[i*solc+0] = iter; | ||
1246 | solp[i*solc+1]=0.; | ||
1247 | } | ||
1248 | gsl_root_fdfsolver_free(s); | ||
1249 | OK | ||
1250 | } | ||
1251 | |||
1252 | |||
1253 | //--------------------------------------------------------------- | ||
1254 | |||
1255 | typedef void TrawfunV(int, double*, int, double*); | ||
1256 | |||
1257 | int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) { | ||
1258 | TrawfunV * f = (TrawfunV*) pars; | ||
1259 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
1260 | double* q = (double*)calloc(y->size,sizeof(double)); | ||
1261 | int k; | ||
1262 | for(k=0;k<x->size;k++) { | ||
1263 | p[k] = gsl_vector_get(x,k); | ||
1264 | } | ||
1265 | f(x->size,p,y->size,q); | ||
1266 | for(k=0;k<y->size;k++) { | ||
1267 | gsl_vector_set(y,k,q[k]); | ||
1268 | } | ||
1269 | free(p); | ||
1270 | free(q); | ||
1271 | return 0; //hmmm | ||
1272 | } | ||
1273 | |||
1274 | int multiroot(int method, void f(int, double*, int, double*), | ||
1275 | double epsabs, int maxit, | ||
1276 | KRVEC(xi), RMAT(sol)) { | ||
1277 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); | ||
1278 | DEBUGMSG("root_only_f"); | ||
1279 | gsl_multiroot_function my_func; | ||
1280 | // extract function from pars | ||
1281 | my_func.f = only_f_aux_multiroot; | ||
1282 | my_func.n = xin; | ||
1283 | my_func.params = f; | ||
1284 | size_t iter = 0; | ||
1285 | int status; | ||
1286 | const gsl_multiroot_fsolver_type *T; | ||
1287 | gsl_multiroot_fsolver *s; | ||
1288 | // Starting point | ||
1289 | KDVVIEW(xi); | ||
1290 | switch(method) { | ||
1291 | case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; } | ||
1292 | case 1 : {T = gsl_multiroot_fsolver_hybrid; break; } | ||
1293 | case 2 : {T = gsl_multiroot_fsolver_dnewton; break; } | ||
1294 | case 3 : {T = gsl_multiroot_fsolver_broyden; break; } | ||
1295 | default: ERROR(BAD_CODE); | ||
1296 | } | ||
1297 | s = gsl_multiroot_fsolver_alloc (T, my_func.n); | ||
1298 | gsl_multiroot_fsolver_set (s, &my_func, V(xi)); | ||
1299 | |||
1300 | do { | ||
1301 | status = gsl_multiroot_fsolver_iterate (s); | ||
1302 | |||
1303 | solp[iter*solc+0] = iter+1; | ||
1304 | |||
1305 | int k; | ||
1306 | for(k=0;k<xin;k++) { | ||
1307 | solp[iter*solc+k+1] = gsl_vector_get(s->x,k); | ||
1308 | } | ||
1309 | for(k=xin;k<2*xin;k++) { | ||
1310 | solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); | ||
1311 | } | ||
1312 | |||
1313 | iter++; | ||
1314 | if (status) /* check if solver is stuck */ | ||
1315 | break; | ||
1316 | |||
1317 | status = | ||
1318 | gsl_multiroot_test_residual (s->f, epsabs); | ||
1319 | } | ||
1320 | while (status == GSL_CONTINUE && iter < maxit); | ||
1321 | |||
1322 | int i,j; | ||
1323 | for (i=iter; i<solr; i++) { | ||
1324 | solp[i*solc+0] = iter; | ||
1325 | for(j=1;j<solc;j++) { | ||
1326 | solp[i*solc+j]=0.; | ||
1327 | } | ||
1328 | } | ||
1329 | gsl_multiroot_fsolver_free(s); | ||
1330 | OK | ||
1331 | } | ||
1332 | |||
1333 | // working with the jacobian | ||
1334 | |||
1335 | typedef struct {int (*f)(int, double*, int, double *); | ||
1336 | int (*jf)(int, double*, int, int, double*);} Tfjf; | ||
1337 | |||
1338 | int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) { | ||
1339 | Tfjf * fjf = ((Tfjf*) pars); | ||
1340 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
1341 | double* q = (double*)calloc(y->size,sizeof(double)); | ||
1342 | int k; | ||
1343 | for(k=0;k<x->size;k++) { | ||
1344 | p[k] = gsl_vector_get(x,k); | ||
1345 | } | ||
1346 | (fjf->f)(x->size,p,y->size,q); | ||
1347 | for(k=0;k<y->size;k++) { | ||
1348 | gsl_vector_set(y,k,q[k]); | ||
1349 | } | ||
1350 | free(p); | ||
1351 | free(q); | ||
1352 | return 0; | ||
1353 | } | ||
1354 | |||
1355 | int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) { | ||
1356 | Tfjf * fjf = ((Tfjf*) pars); | ||
1357 | double* p = (double*)calloc(x->size,sizeof(double)); | ||
1358 | double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double)); | ||
1359 | int i,j,k; | ||
1360 | for(k=0;k<x->size;k++) { | ||
1361 | p[k] = gsl_vector_get(x,k); | ||
1362 | } | ||
1363 | |||
1364 | (fjf->jf)(x->size,p,jac->size1,jac->size2,q); | ||
1365 | |||
1366 | k=0; | ||
1367 | for(i=0;i<jac->size1;i++) { | ||
1368 | for(j=0;j<jac->size2;j++){ | ||
1369 | gsl_matrix_set(jac,i,j,q[k++]); | ||
1370 | } | ||
1371 | } | ||
1372 | free(p); | ||
1373 | free(q); | ||
1374 | return 0; | ||
1375 | } | ||
1376 | |||
1377 | int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { | ||
1378 | f_aux(x,pars,f); | ||
1379 | jf_aux(x,pars,g); | ||
1380 | return 0; | ||
1381 | } | ||
1382 | |||
1383 | int multirootj(int method, int f(int, double*, int, double*), | ||
1384 | int jac(int, double*, int, int, double*), | ||
1385 | double epsabs, int maxit, | ||
1386 | KRVEC(xi), RMAT(sol)) { | ||
1387 | REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); | ||
1388 | DEBUGMSG("root_fjf"); | ||
1389 | gsl_multiroot_function_fdf my_func; | ||
1390 | // extract function from pars | ||
1391 | my_func.f = f_aux; | ||
1392 | my_func.df = jf_aux; | ||
1393 | my_func.fdf = fjf_aux; | ||
1394 | my_func.n = xin; | ||
1395 | Tfjf stfjf; | ||
1396 | stfjf.f = f; | ||
1397 | stfjf.jf = jac; | ||
1398 | my_func.params = &stfjf; | ||
1399 | size_t iter = 0; | ||
1400 | int status; | ||
1401 | const gsl_multiroot_fdfsolver_type *T; | ||
1402 | gsl_multiroot_fdfsolver *s; | ||
1403 | // Starting point | ||
1404 | KDVVIEW(xi); | ||
1405 | switch(method) { | ||
1406 | case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; } | ||
1407 | case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; } | ||
1408 | case 2 : {T = gsl_multiroot_fdfsolver_newton; break; } | ||
1409 | case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; } | ||
1410 | default: ERROR(BAD_CODE); | ||
1411 | } | ||
1412 | s = gsl_multiroot_fdfsolver_alloc (T, my_func.n); | ||
1413 | |||
1414 | gsl_multiroot_fdfsolver_set (s, &my_func, V(xi)); | ||
1415 | |||
1416 | do { | ||
1417 | status = gsl_multiroot_fdfsolver_iterate (s); | ||
1418 | |||
1419 | solp[iter*solc+0] = iter+1; | ||
1420 | |||
1421 | int k; | ||
1422 | for(k=0;k<xin;k++) { | ||
1423 | solp[iter*solc+k+1] = gsl_vector_get(s->x,k); | ||
1424 | } | ||
1425 | for(k=xin;k<2*xin;k++) { | ||
1426 | solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin); | ||
1427 | } | ||
1428 | |||
1429 | iter++; | ||
1430 | if (status) /* check if solver is stuck */ | ||
1431 | break; | ||
1432 | |||
1433 | status = | ||
1434 | gsl_multiroot_test_residual (s->f, epsabs); | ||
1435 | } | ||
1436 | while (status == GSL_CONTINUE && iter < maxit); | ||
1437 | |||
1438 | int i,j; | ||
1439 | for (i=iter; i<solr; i++) { | ||
1440 | solp[i*solc+0] = iter; | ||
1441 | for(j=1;j<solc;j++) { | ||
1442 | solp[i*solc+j]=0.; | ||
1443 | } | ||
1444 | } | ||
1445 | gsl_multiroot_fdfsolver_free(s); | ||
1446 | OK | ||
1447 | } | ||
1448 | |||
1449 | //-------------- non linear least squares fitting ------------------- | ||
1450 | |||
1451 | int nlfit(int method, int f(int, double*, int, double*), | ||
1452 | int jac(int, double*, int, int, double*), | ||
1453 | double epsabs, double epsrel, int maxit, int p, | ||
1454 | KRVEC(xi), RMAT(sol)) { | ||
1455 | REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE); | ||
1456 | DEBUGMSG("nlfit"); | ||
1457 | const gsl_multifit_fdfsolver_type *T; | ||
1458 | gsl_multifit_fdfsolver *s; | ||
1459 | gsl_multifit_function_fdf my_f; | ||
1460 | // extract function from pars | ||
1461 | my_f.f = f_aux; | ||
1462 | my_f.df = jf_aux; | ||
1463 | my_f.fdf = fjf_aux; | ||
1464 | my_f.n = p; | ||
1465 | my_f.p = xin; // !!!! | ||
1466 | Tfjf stfjf; | ||
1467 | stfjf.f = f; | ||
1468 | stfjf.jf = jac; | ||
1469 | my_f.params = &stfjf; | ||
1470 | size_t iter = 0; | ||
1471 | int status; | ||
1472 | |||
1473 | KDVVIEW(xi); | ||
1474 | //DMVIEW(cov); | ||
1475 | |||
1476 | switch(method) { | ||
1477 | case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; } | ||
1478 | case 1 : { T = gsl_multifit_fdfsolver_lmder; break; } | ||
1479 | default: ERROR(BAD_CODE); | ||
1480 | } | ||
1481 | |||
1482 | s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p); | ||
1483 | gsl_multifit_fdfsolver_set (s, &my_f, V(xi)); | ||
1484 | |||
1485 | do { status = gsl_multifit_fdfsolver_iterate (s); | ||
1486 | |||
1487 | solp[iter*solc+0] = iter+1; | ||
1488 | solp[iter*solc+1] = gsl_blas_dnrm2 (s->f); | ||
1489 | |||
1490 | int k; | ||
1491 | for(k=0;k<xin;k++) { | ||
1492 | solp[iter*solc+k+2] = gsl_vector_get(s->x,k); | ||
1493 | } | ||
1494 | |||
1495 | iter++; | ||
1496 | if (status) /* check if solver is stuck */ | ||
1497 | break; | ||
1498 | |||
1499 | status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel); | ||
1500 | } | ||
1501 | while (status == GSL_CONTINUE && iter < maxit); | ||
1502 | |||
1503 | int i,j; | ||
1504 | for (i=iter; i<solr; i++) { | ||
1505 | solp[i*solc+0] = iter; | ||
1506 | for(j=1;j<solc;j++) { | ||
1507 | solp[i*solc+j]=0.; | ||
1508 | } | ||
1509 | } | ||
1510 | |||
1511 | //gsl_multifit_covar (s->J, 0.0, M(cov)); | ||
1512 | |||
1513 | gsl_multifit_fdfsolver_free (s); | ||
1514 | OK | ||
1515 | } | ||
1516 | |||
1517 | |||
1518 | ////////////////////////////////////////////////////// | ||
1519 | |||
1520 | |||
1521 | #define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK } | ||
1522 | |||
1523 | int random_vector(int seed, int code, RVEC(r)) { | ||
1524 | DEBUGMSG("random_vector") | ||
1525 | static gsl_rng * gen = NULL; | ||
1526 | if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);} | ||
1527 | gsl_rng_set (gen, seed); | ||
1528 | int k; | ||
1529 | switch (code) { | ||
1530 | RAN(0,gsl_rng_uniform) | ||
1531 | RAN(1,gsl_ran_ugaussian) | ||
1532 | default: ERROR(BAD_CODE); | ||
1533 | } | ||
1534 | } | ||
1535 | #undef RAN | ||
1536 | |||
1537 | ////////////////////////////////////////////////////// | ||
1538 | |||
1539 | #include "gsl-ode.c" | ||
1540 | |||
1541 | ////////////////////////////////////////////////////// | ||