summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/gsl-aux.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/gsl-aux.c')
-rw-r--r--lib/Numeric/GSL/gsl-aux.c1541
1 files changed, 0 insertions, 1541 deletions
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
deleted file mode 100644
index 410d157..0000000
--- a/lib/Numeric/GSL/gsl-aux.c
+++ /dev/null
@@ -1,1541 +0,0 @@
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
102void no_abort_on_error() {
103 gsl_set_error_handler_off();
104}
105
106
107int 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
117int 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
127int 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
142int 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
157int 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
167int 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
177int 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
194int 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
211int 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
222int 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
233int 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
244int 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
255int 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
273int 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
292int 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
306int 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
321inline 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
331inline 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
341inline 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
348inline 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 }
364int 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
390int 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
417int 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
446int mapC(int code, KCVEC(x), CVEC(r)) {
447 return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
448}
449
450
451gsl_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
469gsl_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 }
494int 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
523int mapQ(int code, KQVEC(x), QVEC(r)) {
524 return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
525}
526
527
528int 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
544int 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
560int 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
576int 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
581int 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
597int 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 }
604int 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
626int 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
648int 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
670int 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 }
676int 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
698int 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
704int 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
724int 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
740int 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
752int 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
765int 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
779int 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
793int 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
806int 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
822int 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
832int 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
843int 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
854int 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
865int 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
876int 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
897typedef double Trawfun(int, double*);
898
899double 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
911double only_f_aux_root(double x, void *pars);
912int 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
965int 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
1024typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
1025
1026double 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
1039void 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
1057void 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
1063int 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
1120double only_f_aux_root(double x, void *pars) {
1121 double (*f)(double) = (double (*)(double)) pars;
1122 return f(x);
1123}
1124
1125int 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
1176typedef struct {
1177 double (*f)(double);
1178 double (*jf)(double);
1179} uniTfjf;
1180
1181double f_aux_uni(double x, void *pars) {
1182 uniTfjf * fjf = ((uniTfjf*) pars);
1183 return (fjf->f)(x);
1184}
1185
1186double jf_aux_uni(double x, void * pars) {
1187 uniTfjf * fjf = ((uniTfjf*) pars);
1188 return (fjf->jf)(x);
1189}
1190
1191void 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
1196int 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
1255typedef void TrawfunV(int, double*, int, double*);
1256
1257int 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
1274int 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
1335typedef struct {int (*f)(int, double*, int, double *);
1336 int (*jf)(int, double*, int, int, double*);} Tfjf;
1337
1338int 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
1355int 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
1377int 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
1383int 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
1451int 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
1523int 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//////////////////////////////////////////////////////
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-07 06:25:16 +0000