diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-14 19:46:44 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-14 19:46:44 +0200 |
commit | 0f9575462eb37a7c9985583ca33ae315f6e6431d (patch) | |
tree | c8f4c6519c1a4b6969e4ae80fed9e410b852f7d7 /packages/base/src/C | |
parent | 494fe8f3e66f52244c4e0aedfd00eb15b9caceef (diff) |
vectorized operations in base
Diffstat (limited to 'packages/base/src/C')
-rw-r--r-- | packages/base/src/C/vector-aux.c | 676 |
1 files changed, 676 insertions, 0 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c new file mode 100644 index 0000000..07ccf76 --- /dev/null +++ b/packages/base/src/C/vector-aux.c | |||
@@ -0,0 +1,676 @@ | |||
1 | #include "lapack-aux.h" | ||
2 | |||
3 | #define V(x) x##n,x##p | ||
4 | |||
5 | #include <string.h> | ||
6 | #include <math.h> | ||
7 | #include <stdio.h> | ||
8 | |||
9 | #define MACRO(B) do {B} while (0) | ||
10 | #define ERROR(CODE) MACRO(return CODE;) | ||
11 | #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) | ||
12 | #define OK return 0; | ||
13 | |||
14 | #define MIN(A,B) ((A)<(B)?(A):(B)) | ||
15 | #define MAX(A,B) ((A)>(B)?(A):(B)) | ||
16 | |||
17 | #ifdef DBG | ||
18 | #define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); | ||
19 | #else | ||
20 | #define DEBUGMSG(M) | ||
21 | #endif | ||
22 | |||
23 | #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) | ||
24 | |||
25 | #ifdef DBG | ||
26 | #define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); | ||
27 | #else | ||
28 | #define DEBUGMAT(MSG,X) | ||
29 | #endif | ||
30 | |||
31 | #ifdef DBG | ||
32 | #define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); | ||
33 | #else | ||
34 | #define DEBUGVEC(MSG,X) | ||
35 | #endif | ||
36 | |||
37 | |||
38 | #define BAD_SIZE 2000 | ||
39 | #define BAD_CODE 2001 | ||
40 | #define MEM 2002 | ||
41 | #define BAD_FILE 2003 | ||
42 | |||
43 | |||
44 | int sumF(KFVEC(x),FVEC(r)) { | ||
45 | DEBUGMSG("sumF"); | ||
46 | printf("hello!\n"); | ||
47 | REQUIRES(rn==1,BAD_SIZE); | ||
48 | int i; | ||
49 | float res = 0; | ||
50 | for (i = 0; i < xn; i++) res += xp[i]; | ||
51 | rp[0] = res; | ||
52 | OK | ||
53 | } | ||
54 | |||
55 | int sumR(KDVEC(x),DVEC(r)) { | ||
56 | DEBUGMSG("sumR"); | ||
57 | REQUIRES(rn==1,BAD_SIZE); | ||
58 | int i; | ||
59 | double res = 0; | ||
60 | for (i = 0; i < xn; i++) res += xp[i]; | ||
61 | rp[0] = res; | ||
62 | OK | ||
63 | } | ||
64 | |||
65 | |||
66 | int sumQ(KQVEC(x),QVEC(r)) { | ||
67 | DEBUGMSG("sumQ"); | ||
68 | REQUIRES(rn==1,BAD_SIZE); | ||
69 | int i; | ||
70 | complex res; | ||
71 | res.r = 0; | ||
72 | res.i = 0; | ||
73 | for (i = 0; i < xn; i++) { | ||
74 | res.r += xp[i].r; | ||
75 | res.i += xp[i].i; | ||
76 | } | ||
77 | rp[0] = res; | ||
78 | OK | ||
79 | } | ||
80 | |||
81 | int sumC(KCVEC(x),CVEC(r)) { | ||
82 | DEBUGMSG("sumC"); | ||
83 | REQUIRES(rn==1,BAD_SIZE); | ||
84 | int i; | ||
85 | doublecomplex res; | ||
86 | res.r = 0; | ||
87 | res.i = 0; | ||
88 | for (i = 0; i < xn; i++) { | ||
89 | res.r += xp[i].r; | ||
90 | res.i += xp[i].i; | ||
91 | } | ||
92 | rp[0] = res; | ||
93 | OK | ||
94 | } | ||
95 | |||
96 | |||
97 | int prodF(KFVEC(x),FVEC(r)) { | ||
98 | DEBUGMSG("prodF"); | ||
99 | REQUIRES(rn==1,BAD_SIZE); | ||
100 | int i; | ||
101 | float res = 1; | ||
102 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
103 | rp[0] = res; | ||
104 | OK | ||
105 | } | ||
106 | |||
107 | int prodR(KDVEC(x),DVEC(r)) { | ||
108 | DEBUGMSG("prodR"); | ||
109 | REQUIRES(rn==1,BAD_SIZE); | ||
110 | int i; | ||
111 | double res = 1; | ||
112 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
113 | rp[0] = res; | ||
114 | OK | ||
115 | } | ||
116 | |||
117 | |||
118 | int prodQ(KQVEC(x),QVEC(r)) { | ||
119 | DEBUGMSG("prodQ"); | ||
120 | REQUIRES(rn==1,BAD_SIZE); | ||
121 | int i; | ||
122 | complex res; | ||
123 | float temp; | ||
124 | res.r = 1; | ||
125 | res.i = 0; | ||
126 | for (i = 0; i < xn; i++) { | ||
127 | temp = res.r * xp[i].r - res.i * xp[i].i; | ||
128 | res.i = res.r * xp[i].i + res.i * xp[i].r; | ||
129 | res.r = temp; | ||
130 | } | ||
131 | rp[0] = res; | ||
132 | OK | ||
133 | } | ||
134 | |||
135 | int prodC(KCVEC(x),CVEC(r)) { | ||
136 | DEBUGMSG("prodC"); | ||
137 | REQUIRES(rn==1,BAD_SIZE); | ||
138 | int i; | ||
139 | doublecomplex res; | ||
140 | double temp; | ||
141 | res.r = 1; | ||
142 | res.i = 0; | ||
143 | for (i = 0; i < xn; i++) { | ||
144 | temp = res.r * xp[i].r - res.i * xp[i].i; | ||
145 | res.i = res.r * xp[i].i + res.i * xp[i].r; | ||
146 | res.r = temp; | ||
147 | } | ||
148 | rp[0] = res; | ||
149 | OK | ||
150 | } | ||
151 | |||
152 | |||
153 | double dnrm2_(integer*, const double*, integer*); | ||
154 | double dasum_(integer*, const double*, integer*); | ||
155 | |||
156 | double vector_max(KDVEC(x)) { | ||
157 | double r = xp[0]; | ||
158 | int k; | ||
159 | for (k = 1; k<xn; k++) { | ||
160 | if(xp[k]>r) { | ||
161 | r = xp[k]; | ||
162 | } | ||
163 | } | ||
164 | return r; | ||
165 | } | ||
166 | |||
167 | double vector_min(KDVEC(x)) { | ||
168 | double r = xp[0]; | ||
169 | int k; | ||
170 | for (k = 1; k<xn; k++) { | ||
171 | if(xp[k]<r) { | ||
172 | r = xp[k]; | ||
173 | } | ||
174 | } | ||
175 | return r; | ||
176 | } | ||
177 | |||
178 | double vector_max_index(KDVEC(x)) { | ||
179 | int k, r = 0; | ||
180 | for (k = 1; k<xn; k++) { | ||
181 | if(xp[k]>xp[0]) { | ||
182 | r = k; | ||
183 | } | ||
184 | } | ||
185 | return r; | ||
186 | } | ||
187 | |||
188 | double vector_min_index(KDVEC(x)) { | ||
189 | int k, r = 0; | ||
190 | for (k = 1; k<xn; k++) { | ||
191 | if(xp[k]<xp[0]) { | ||
192 | r = k; | ||
193 | } | ||
194 | } | ||
195 | return r; | ||
196 | } | ||
197 | |||
198 | int toScalarR(int code, KDVEC(x), DVEC(r)) { | ||
199 | REQUIRES(rn==1,BAD_SIZE); | ||
200 | DEBUGMSG("toScalarR"); | ||
201 | double res; | ||
202 | integer one = 1; | ||
203 | integer n = rn; | ||
204 | switch(code) { | ||
205 | case 0: { res = dnrm2_(&n,xp,&one); break; } | ||
206 | case 1: { res = dasum_(&n,xp,&one); break; } | ||
207 | case 2: { res = vector_max_index(V(x)); break; } | ||
208 | case 3: { res = vector_max(V(x)); break; } | ||
209 | case 4: { res = vector_min_index(V(x)); break; } | ||
210 | case 5: { res = vector_min(V(x)); break; } | ||
211 | default: ERROR(BAD_CODE); | ||
212 | } | ||
213 | rp[0] = res; | ||
214 | OK | ||
215 | } | ||
216 | |||
217 | |||
218 | float snrm2_(integer*, const float*, integer*); | ||
219 | float sasum_(integer*, const float*, integer*); | ||
220 | |||
221 | float vector_max_f(KFVEC(x)) { | ||
222 | float r = xp[0]; | ||
223 | int k; | ||
224 | for (k = 1; k<xn; k++) { | ||
225 | if(xp[k]>r) { | ||
226 | r = xp[k]; | ||
227 | } | ||
228 | } | ||
229 | return r; | ||
230 | } | ||
231 | |||
232 | float vector_min_f(KFVEC(x)) { | ||
233 | float r = xp[0]; | ||
234 | int k; | ||
235 | for (k = 1; k<xn; k++) { | ||
236 | if(xp[k]<r) { | ||
237 | r = xp[k]; | ||
238 | } | ||
239 | } | ||
240 | return r; | ||
241 | } | ||
242 | |||
243 | float vector_max_index_f(KFVEC(x)) { | ||
244 | int k, r = 0; | ||
245 | for (k = 1; k<xn; k++) { | ||
246 | if(xp[k]>xp[0]) { | ||
247 | r = k; | ||
248 | } | ||
249 | } | ||
250 | return r; | ||
251 | } | ||
252 | |||
253 | float vector_min_index_f(KFVEC(x)) { | ||
254 | int k, r = 0; | ||
255 | for (k = 1; k<xn; k++) { | ||
256 | if(xp[k]<xp[0]) { | ||
257 | r = k; | ||
258 | } | ||
259 | } | ||
260 | return r; | ||
261 | } | ||
262 | |||
263 | |||
264 | int toScalarF(int code, KFVEC(x), FVEC(r)) { | ||
265 | REQUIRES(rn==1,BAD_SIZE); | ||
266 | DEBUGMSG("toScalarF"); | ||
267 | float res; | ||
268 | integer one = 1; | ||
269 | integer n = rn; | ||
270 | switch(code) { | ||
271 | case 0: { res = snrm2_(&n,xp,&one); break; } | ||
272 | case 1: { res = sasum_(&n,xp,&one); break; } | ||
273 | case 2: { res = vector_max_index_f(V(x)); break; } | ||
274 | case 3: { res = vector_max_f(V(x)); break; } | ||
275 | case 4: { res = vector_min_index_f(V(x)); break; } | ||
276 | case 5: { res = vector_min_f(V(x)); break; } | ||
277 | default: ERROR(BAD_CODE); | ||
278 | } | ||
279 | rp[0] = res; | ||
280 | OK | ||
281 | } | ||
282 | |||
283 | double dznrm2_(integer*, const doublecomplex*, integer*); | ||
284 | double dzasum_(integer*, const doublecomplex*, integer*); | ||
285 | |||
286 | int toScalarC(int code, KCVEC(x), DVEC(r)) { | ||
287 | REQUIRES(rn==1,BAD_SIZE); | ||
288 | DEBUGMSG("toScalarC"); | ||
289 | double res; | ||
290 | integer one = 1; | ||
291 | integer n = rn; | ||
292 | switch(code) { | ||
293 | case 0: { res = dznrm2_(&n,xp,&one); break; } | ||
294 | case 1: { res = dzasum_(&n,xp,&one); break; } | ||
295 | default: ERROR(BAD_CODE); | ||
296 | } | ||
297 | rp[0] = res; | ||
298 | OK | ||
299 | } | ||
300 | |||
301 | |||
302 | double scnrm2_(integer*, const complex*, integer*); | ||
303 | double scasum_(integer*, const complex*, integer*); | ||
304 | |||
305 | int toScalarQ(int code, KQVEC(x), FVEC(r)) { | ||
306 | REQUIRES(rn==1,BAD_SIZE); | ||
307 | DEBUGMSG("toScalarQ"); | ||
308 | float res; | ||
309 | integer one = 1; | ||
310 | integer n = rn; | ||
311 | switch(code) { | ||
312 | case 0: { res = scnrm2_(&n,xp,&one); break; } | ||
313 | case 1: { res = scasum_(&n,xp,&one); 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 | |||
342 | #define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK } | ||
343 | #define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK } | ||
344 | int mapR(int code, KDVEC(x), DVEC(r)) { | ||
345 | int k; | ||
346 | REQUIRES(xn == rn,BAD_SIZE); | ||
347 | DEBUGMSG("mapR"); | ||
348 | switch (code) { | ||
349 | OP(0,sin) | ||
350 | OP(1,cos) | ||
351 | OP(2,tan) | ||
352 | OP(3,fabs) | ||
353 | OP(4,asin) | ||
354 | OP(5,acos) | ||
355 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
356 | OP(7,sinh) | ||
357 | OP(8,cosh) | ||
358 | OP(9,tanh) | ||
359 | // OP(10,gsl_asinh) | ||
360 | // OP(11,gsl_acosh) | ||
361 | // OP(12,gsl_atanh) | ||
362 | OP(13,exp) | ||
363 | OP(14,log) | ||
364 | OP(15,sign) | ||
365 | OP(16,sqrt) | ||
366 | default: ERROR(BAD_CODE); | ||
367 | } | ||
368 | } | ||
369 | |||
370 | int mapF(int code, KFVEC(x), FVEC(r)) { | ||
371 | int k; | ||
372 | REQUIRES(xn == rn,BAD_SIZE); | ||
373 | DEBUGMSG("mapF"); | ||
374 | switch (code) { | ||
375 | OP(0,sin) | ||
376 | OP(1,cos) | ||
377 | OP(2,tan) | ||
378 | OP(3,fabs) | ||
379 | OP(4,asin) | ||
380 | OP(5,acos) | ||
381 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
382 | OP(7,sinh) | ||
383 | OP(8,cosh) | ||
384 | OP(9,tanh) | ||
385 | // OP(10,gsl_asinh) | ||
386 | // OP(11,gsl_acosh) | ||
387 | // OP(12,gsl_atanh) | ||
388 | OP(13,exp) | ||
389 | OP(14,log) | ||
390 | OP(15,sign) | ||
391 | OP(16,sqrt) | ||
392 | default: ERROR(BAD_CODE); | ||
393 | } | ||
394 | } | ||
395 | |||
396 | |||
397 | inline double abs_complex(doublecomplex z) { | ||
398 | return sqrt(z.r*z.r + z.i*z.i); | ||
399 | } | ||
400 | |||
401 | inline doublecomplex complex_abs_complex(doublecomplex z) { | ||
402 | doublecomplex r; | ||
403 | r.r = abs_complex(z); | ||
404 | r.i = 0; | ||
405 | return r; | ||
406 | } | ||
407 | |||
408 | inline doublecomplex complex_signum_complex(doublecomplex z) { | ||
409 | doublecomplex r; | ||
410 | double mag; | ||
411 | if (z.r == 0 && z.i == 0) { | ||
412 | r.r = 0; | ||
413 | r.i = 0; | ||
414 | } else { | ||
415 | mag = abs_complex(z); | ||
416 | r.r = z.r/mag; | ||
417 | r.i = z.i/mag; | ||
418 | } | ||
419 | return r; | ||
420 | } | ||
421 | |||
422 | |||
423 | |||
424 | int mapC(int code, KCVEC(x), CVEC(r)) { | ||
425 | int k; | ||
426 | REQUIRES(xn == rn,BAD_SIZE); | ||
427 | DEBUGMSG("mapC"); | ||
428 | switch (code) { | ||
429 | /* | ||
430 | OP(0,gsl_complex_sin) | ||
431 | OP(1,gsl_complex_cos) | ||
432 | OP(2,gsl_complex_tan) | ||
433 | */ | ||
434 | OP(3,complex_abs_complex) | ||
435 | /* | ||
436 | OP(4,gsl_complex_arcsin) | ||
437 | OP(5,gsl_complex_arccos) | ||
438 | OP(6,gsl_complex_arctan) | ||
439 | OP(7,gsl_complex_sinh) | ||
440 | OP(8,gsl_complex_cosh) | ||
441 | OP(9,gsl_complex_tanh) | ||
442 | OP(10,gsl_complex_arcsinh) | ||
443 | OP(11,gsl_complex_arccosh) | ||
444 | OP(12,gsl_complex_arctanh) | ||
445 | OP(13,gsl_complex_exp) | ||
446 | OP(14,gsl_complex_log) | ||
447 | */ | ||
448 | OP(15,complex_signum_complex) | ||
449 | |||
450 | // OP(16,gsl_complex_sqrt) | ||
451 | default: ERROR(BAD_CODE); | ||
452 | } | ||
453 | } | ||
454 | |||
455 | |||
456 | |||
457 | inline complex complex_f_math_fun(doublecomplex (*cf)(doublecomplex), complex a) | ||
458 | { | ||
459 | doublecomplex c; | ||
460 | doublecomplex r; | ||
461 | |||
462 | complex float_r; | ||
463 | |||
464 | c.r = a.r; | ||
465 | c.i = a.i; | ||
466 | |||
467 | r = (*cf)(c); | ||
468 | |||
469 | float_r.r = r.r; | ||
470 | float_r.i = r.i; | ||
471 | |||
472 | return float_r; | ||
473 | } | ||
474 | |||
475 | inline complex complex_f_math_op(doublecomplex (*cf)(doublecomplex,doublecomplex), | ||
476 | complex a,complex b) | ||
477 | { | ||
478 | doublecomplex c1,c2,r; | ||
479 | |||
480 | complex float_r; | ||
481 | |||
482 | c1.r = a.r; | ||
483 | c1.i = a.i; | ||
484 | |||
485 | c2.r = b.r; | ||
486 | c2.i = b.i; | ||
487 | |||
488 | r = (*cf)(c1,c2); | ||
489 | |||
490 | float_r.r = r.r; | ||
491 | float_r.i = r.i; | ||
492 | |||
493 | return float_r; | ||
494 | } | ||
495 | |||
496 | #define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK } | ||
497 | #define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_op(&F,A,B); OK } | ||
498 | int mapQ(int code, KQVEC(x), QVEC(r)) { | ||
499 | int k; | ||
500 | REQUIRES(xn == rn,BAD_SIZE); | ||
501 | DEBUGMSG("mapQ"); | ||
502 | switch (code) { | ||
503 | /* | ||
504 | OPC(0,gsl_complex_sin) | ||
505 | OPC(1,gsl_complex_cos) | ||
506 | OPC(2,gsl_complex_tan) | ||
507 | */ | ||
508 | OPC(3,complex_abs_complex) | ||
509 | /* | ||
510 | OPC(4,gsl_complex_arcsin) | ||
511 | OPC(5,gsl_complex_arccos) | ||
512 | OPC(6,gsl_complex_arctan) | ||
513 | OPC(7,gsl_complex_sinh) | ||
514 | OPC(8,gsl_complex_cosh) | ||
515 | OPC(9,gsl_complex_tanh) | ||
516 | OPC(10,gsl_complex_arcsinh) | ||
517 | OPC(11,gsl_complex_arccosh) | ||
518 | OPC(12,gsl_complex_arctanh) | ||
519 | OPC(13,gsl_complex_exp) | ||
520 | OPC(14,gsl_complex_log) | ||
521 | */ | ||
522 | OPC(15,complex_signum_complex) | ||
523 | |||
524 | // OPC(16,gsl_complex_sqrt) | ||
525 | default: ERROR(BAD_CODE); | ||
526 | } | ||
527 | } | ||
528 | |||
529 | |||
530 | int mapValR(int code, double* pval, KDVEC(x), DVEC(r)) { | ||
531 | int k; | ||
532 | double val = *pval; | ||
533 | REQUIRES(xn == rn,BAD_SIZE); | ||
534 | DEBUGMSG("mapValR"); | ||
535 | switch (code) { | ||
536 | OPV(0,val*xp[k]) | ||
537 | OPV(1,val/xp[k]) | ||
538 | OPV(2,val+xp[k]) | ||
539 | OPV(3,val-xp[k]) | ||
540 | OPV(4,pow(val,xp[k])) | ||
541 | OPV(5,pow(xp[k],val)) | ||
542 | default: ERROR(BAD_CODE); | ||
543 | } | ||
544 | } | ||
545 | |||
546 | int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) { | ||
547 | int k; | ||
548 | float val = *pval; | ||
549 | REQUIRES(xn == rn,BAD_SIZE); | ||
550 | DEBUGMSG("mapValF"); | ||
551 | switch (code) { | ||
552 | OPV(0,val*xp[k]) | ||
553 | OPV(1,val/xp[k]) | ||
554 | OPV(2,val+xp[k]) | ||
555 | OPV(3,val-xp[k]) | ||
556 | OPV(4,pow(val,xp[k])) | ||
557 | OPV(5,pow(xp[k],val)) | ||
558 | default: ERROR(BAD_CODE); | ||
559 | } | ||
560 | } | ||
561 | |||
562 | |||
563 | |||
564 | inline doublecomplex complex_add(doublecomplex a, doublecomplex b) { | ||
565 | doublecomplex r; | ||
566 | r.r = a.r+b.r; | ||
567 | r.i = a.i+b.i; | ||
568 | return r; | ||
569 | } | ||
570 | |||
571 | |||
572 | int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) { | ||
573 | int k; | ||
574 | doublecomplex val = *pval; | ||
575 | REQUIRES(xn == rn,BAD_SIZE); | ||
576 | DEBUGMSG("mapValC"); | ||
577 | switch (code) { | ||
578 | // OPV(0,gsl_complex_mul(val,xp[k])) | ||
579 | // OPV(1,gsl_complex_div(val,xp[k])) | ||
580 | OPV(2,complex_add(val,xp[k])) | ||
581 | // OPV(3,gsl_complex_sub(val,xp[k])) | ||
582 | // OPV(4,gsl_complex_pow(val,xp[k])) | ||
583 | // OPV(5,gsl_complex_pow(xp[k],val)) | ||
584 | default: ERROR(BAD_CODE); | ||
585 | } | ||
586 | } | ||
587 | |||
588 | |||
589 | int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) { | ||
590 | int k; | ||
591 | complex val = *pval; | ||
592 | REQUIRES(xn == rn,BAD_SIZE); | ||
593 | DEBUGMSG("mapValQ"); | ||
594 | switch (code) { | ||
595 | // OPCA(0,gsl_complex_mul,val,xp[k]) | ||
596 | // OPCA(1,gsl_complex_div,val,xp[k]) | ||
597 | OPCA(2,complex_add,val,xp[k]) | ||
598 | // OPCA(3,gsl_complex_sub,val,xp[k]) | ||
599 | // OPCA(4,gsl_complex_pow,val,xp[k]) | ||
600 | // OPCA(5,gsl_complex_pow,xp[k],val) | ||
601 | default: ERROR(BAD_CODE); | ||
602 | } | ||
603 | } | ||
604 | |||
605 | |||
606 | |||
607 | #define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK } | ||
608 | #define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK } | ||
609 | #define OPZO(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = ap[k] O bp[k]; OK } | ||
610 | |||
611 | int zipR(int code, KDVEC(a), KDVEC(b), DVEC(r)) { | ||
612 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
613 | int k; | ||
614 | switch(code) { | ||
615 | OPZO(0,"zipR Add",+) | ||
616 | OPZO(1,"zipR Sub",-) | ||
617 | OPZO(2,"zipR Mul",*) | ||
618 | OPZO(3,"zipR Div",/) | ||
619 | OPZE(4,"zipR Pow", pow) | ||
620 | OPZE(5,"zipR ATan2",atan2) | ||
621 | default: ERROR(BAD_CODE); | ||
622 | } | ||
623 | } | ||
624 | |||
625 | int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) { | ||
626 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
627 | int k; | ||
628 | switch(code) { | ||
629 | OPZO(0,"zipR Add",+) | ||
630 | OPZO(1,"zipR Sub",-) | ||
631 | OPZO(2,"zipR Mul",*) | ||
632 | OPZO(3,"zipR Div",/) | ||
633 | OPZE(4,"zipR Pow", pow) | ||
634 | OPZE(5,"zipR ATan2",atan2) | ||
635 | default: ERROR(BAD_CODE); | ||
636 | } | ||
637 | } | ||
638 | |||
639 | |||
640 | |||
641 | |||
642 | int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { | ||
643 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
644 | int k; | ||
645 | switch(code) { | ||
646 | OPZE(0,"zipC Add",complex_add) | ||
647 | // OPZE(1,"zipC Sub",gsl_complex_sub) | ||
648 | // OPZE(2,"zipC Mul",gsl_complex_mul) | ||
649 | // OPZE(3,"zipC Div",gsl_complex_div) | ||
650 | // OPZE(4,"zipC Pow",gsl_complex_pow) | ||
651 | // //OPZE(5,"zipR ATan2",atan2) | ||
652 | default: ERROR(BAD_CODE); | ||
653 | } | ||
654 | } | ||
655 | |||
656 | |||
657 | |||
658 | |||
659 | |||
660 | #define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_f_math_op(&E,ap[k],bp[k]); OK } | ||
661 | |||
662 | int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) { | ||
663 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
664 | int k; | ||
665 | switch(code) { | ||
666 | OPCZE(0,"zipQ Add",complex_add) | ||
667 | // OPCZE(1,"zipQ Sub",gsl_complex_sub) | ||
668 | // OPCZE(2,"zipQ Mul",gsl_complex_mul) | ||
669 | // OPCZE(3,"zipQ Div",gsl_complex_div) | ||
670 | // OPCZE(4,"zipQ Pow",gsl_complex_pow) | ||
671 | // //OPZE(5,"zipR ATan2",atan2) | ||
672 | default: ERROR(BAD_CODE); | ||
673 | } | ||
674 | } | ||
675 | |||
676 | |||