diff options
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/gsl-vector.c')
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/gsl-vector.c | 642 |
1 files changed, 642 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-vector.c b/packages/hmatrix/src/Numeric/GSL/gsl-vector.c new file mode 100644 index 0000000..40a086a --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/gsl-vector.c | |||
@@ -0,0 +1,642 @@ | |||
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_complex_math.h> | ||
27 | #include <string.h> | ||
28 | #include <stdio.h> | ||
29 | |||
30 | #define MACRO(B) do {B} while (0) | ||
31 | #define ERROR(CODE) MACRO(return CODE;) | ||
32 | #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) | ||
33 | #define OK return 0; | ||
34 | |||
35 | #define MIN(A,B) ((A)<(B)?(A):(B)) | ||
36 | #define MAX(A,B) ((A)>(B)?(A):(B)) | ||
37 | |||
38 | #ifdef DBG | ||
39 | #define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M); | ||
40 | #else | ||
41 | #define DEBUGMSG(M) | ||
42 | #endif | ||
43 | |||
44 | #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) | ||
45 | |||
46 | #ifdef DBG | ||
47 | #define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n"); | ||
48 | #else | ||
49 | #define DEBUGMAT(MSG,X) | ||
50 | #endif | ||
51 | |||
52 | #ifdef DBG | ||
53 | #define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n"); | ||
54 | #else | ||
55 | #define DEBUGVEC(MSG,X) | ||
56 | #endif | ||
57 | |||
58 | #define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n) | ||
59 | #define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c) | ||
60 | #define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n) | ||
61 | #define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c) | ||
62 | #define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n) | ||
63 | #define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c) | ||
64 | #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) | ||
65 | #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) | ||
66 | |||
67 | #define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) | ||
68 | #define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) | ||
69 | #define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) | ||
70 | #define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) | ||
71 | #define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) | ||
72 | #define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) | ||
73 | #define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) | ||
74 | #define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) | ||
75 | |||
76 | #define V(a) (&a.vector) | ||
77 | #define M(a) (&a.matrix) | ||
78 | |||
79 | #define GCVEC(A) int A##n, gsl_complex*A##p | ||
80 | #define KGCVEC(A) int A##n, const gsl_complex*A##p | ||
81 | |||
82 | #define GQVEC(A) int A##n, gsl_complex_float*A##p | ||
83 | #define KGQVEC(A) int A##n, const gsl_complex_float*A##p | ||
84 | |||
85 | #define BAD_SIZE 2000 | ||
86 | #define BAD_CODE 2001 | ||
87 | #define MEM 2002 | ||
88 | #define BAD_FILE 2003 | ||
89 | |||
90 | |||
91 | int sumF(KFVEC(x),FVEC(r)) { | ||
92 | DEBUGMSG("sumF"); | ||
93 | REQUIRES(rn==1,BAD_SIZE); | ||
94 | int i; | ||
95 | float res = 0; | ||
96 | for (i = 0; i < xn; i++) res += xp[i]; | ||
97 | rp[0] = res; | ||
98 | OK | ||
99 | } | ||
100 | |||
101 | int sumR(KRVEC(x),RVEC(r)) { | ||
102 | DEBUGMSG("sumR"); | ||
103 | REQUIRES(rn==1,BAD_SIZE); | ||
104 | int i; | ||
105 | double res = 0; | ||
106 | for (i = 0; i < xn; i++) res += xp[i]; | ||
107 | rp[0] = res; | ||
108 | OK | ||
109 | } | ||
110 | |||
111 | int sumQ(KQVEC(x),QVEC(r)) { | ||
112 | DEBUGMSG("sumQ"); | ||
113 | REQUIRES(rn==1,BAD_SIZE); | ||
114 | int i; | ||
115 | gsl_complex_float res; | ||
116 | res.dat[0] = 0; | ||
117 | res.dat[1] = 0; | ||
118 | for (i = 0; i < xn; i++) { | ||
119 | res.dat[0] += xp[i].dat[0]; | ||
120 | res.dat[1] += xp[i].dat[1]; | ||
121 | } | ||
122 | rp[0] = res; | ||
123 | OK | ||
124 | } | ||
125 | |||
126 | int sumC(KCVEC(x),CVEC(r)) { | ||
127 | DEBUGMSG("sumC"); | ||
128 | REQUIRES(rn==1,BAD_SIZE); | ||
129 | int i; | ||
130 | gsl_complex res; | ||
131 | res.dat[0] = 0; | ||
132 | res.dat[1] = 0; | ||
133 | for (i = 0; i < xn; i++) { | ||
134 | res.dat[0] += xp[i].dat[0]; | ||
135 | res.dat[1] += xp[i].dat[1]; | ||
136 | } | ||
137 | rp[0] = res; | ||
138 | OK | ||
139 | } | ||
140 | |||
141 | int prodF(KFVEC(x),FVEC(r)) { | ||
142 | DEBUGMSG("prodF"); | ||
143 | REQUIRES(rn==1,BAD_SIZE); | ||
144 | int i; | ||
145 | float res = 1; | ||
146 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
147 | rp[0] = res; | ||
148 | OK | ||
149 | } | ||
150 | |||
151 | int prodR(KRVEC(x),RVEC(r)) { | ||
152 | DEBUGMSG("prodR"); | ||
153 | REQUIRES(rn==1,BAD_SIZE); | ||
154 | int i; | ||
155 | double res = 1; | ||
156 | for (i = 0; i < xn; i++) res *= xp[i]; | ||
157 | rp[0] = res; | ||
158 | OK | ||
159 | } | ||
160 | |||
161 | int prodQ(KQVEC(x),QVEC(r)) { | ||
162 | DEBUGMSG("prodQ"); | ||
163 | REQUIRES(rn==1,BAD_SIZE); | ||
164 | int i; | ||
165 | gsl_complex_float res; | ||
166 | float temp; | ||
167 | res.dat[0] = 1; | ||
168 | res.dat[1] = 0; | ||
169 | for (i = 0; i < xn; i++) { | ||
170 | temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; | ||
171 | res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; | ||
172 | res.dat[0] = temp; | ||
173 | } | ||
174 | rp[0] = res; | ||
175 | OK | ||
176 | } | ||
177 | |||
178 | int prodC(KCVEC(x),CVEC(r)) { | ||
179 | DEBUGMSG("prodC"); | ||
180 | REQUIRES(rn==1,BAD_SIZE); | ||
181 | int i; | ||
182 | gsl_complex res; | ||
183 | double temp; | ||
184 | res.dat[0] = 1; | ||
185 | res.dat[1] = 0; | ||
186 | for (i = 0; i < xn; i++) { | ||
187 | temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; | ||
188 | res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; | ||
189 | res.dat[0] = temp; | ||
190 | } | ||
191 | rp[0] = res; | ||
192 | OK | ||
193 | } | ||
194 | |||
195 | |||
196 | int toScalarR(int code, KRVEC(x), RVEC(r)) { | ||
197 | REQUIRES(rn==1,BAD_SIZE); | ||
198 | DEBUGMSG("toScalarR"); | ||
199 | KDVVIEW(x); | ||
200 | double res; | ||
201 | switch(code) { | ||
202 | case 0: { res = gsl_blas_dnrm2(V(x)); break; } | ||
203 | case 1: { res = gsl_blas_dasum(V(x)); break; } | ||
204 | case 2: { res = gsl_vector_max_index(V(x)); break; } | ||
205 | case 3: { res = gsl_vector_max(V(x)); break; } | ||
206 | case 4: { res = gsl_vector_min_index(V(x)); break; } | ||
207 | case 5: { res = gsl_vector_min(V(x)); break; } | ||
208 | default: ERROR(BAD_CODE); | ||
209 | } | ||
210 | rp[0] = res; | ||
211 | OK | ||
212 | } | ||
213 | |||
214 | int toScalarF(int code, KFVEC(x), FVEC(r)) { | ||
215 | REQUIRES(rn==1,BAD_SIZE); | ||
216 | DEBUGMSG("toScalarF"); | ||
217 | KFVVIEW(x); | ||
218 | float res; | ||
219 | switch(code) { | ||
220 | case 0: { res = gsl_blas_snrm2(V(x)); break; } | ||
221 | case 1: { res = gsl_blas_sasum(V(x)); break; } | ||
222 | case 2: { res = gsl_vector_float_max_index(V(x)); break; } | ||
223 | case 3: { res = gsl_vector_float_max(V(x)); break; } | ||
224 | case 4: { res = gsl_vector_float_min_index(V(x)); break; } | ||
225 | case 5: { res = gsl_vector_float_min(V(x)); break; } | ||
226 | default: ERROR(BAD_CODE); | ||
227 | } | ||
228 | rp[0] = res; | ||
229 | OK | ||
230 | } | ||
231 | |||
232 | |||
233 | int toScalarC(int code, KCVEC(x), RVEC(r)) { | ||
234 | REQUIRES(rn==1,BAD_SIZE); | ||
235 | DEBUGMSG("toScalarC"); | ||
236 | KCVVIEW(x); | ||
237 | double res; | ||
238 | switch(code) { | ||
239 | case 0: { res = gsl_blas_dznrm2(V(x)); break; } | ||
240 | case 1: { res = gsl_blas_dzasum(V(x)); break; } | ||
241 | default: ERROR(BAD_CODE); | ||
242 | } | ||
243 | rp[0] = res; | ||
244 | OK | ||
245 | } | ||
246 | |||
247 | int toScalarQ(int code, KQVEC(x), FVEC(r)) { | ||
248 | REQUIRES(rn==1,BAD_SIZE); | ||
249 | DEBUGMSG("toScalarQ"); | ||
250 | KQVVIEW(x); | ||
251 | float res; | ||
252 | switch(code) { | ||
253 | case 0: { res = gsl_blas_scnrm2(V(x)); break; } | ||
254 | case 1: { res = gsl_blas_scasum(V(x)); break; } | ||
255 | default: ERROR(BAD_CODE); | ||
256 | } | ||
257 | rp[0] = res; | ||
258 | OK | ||
259 | } | ||
260 | |||
261 | |||
262 | inline double sign(double x) { | ||
263 | if(x>0) { | ||
264 | return +1.0; | ||
265 | } else if (x<0) { | ||
266 | return -1.0; | ||
267 | } else { | ||
268 | return 0.0; | ||
269 | } | ||
270 | } | ||
271 | |||
272 | inline float float_sign(float x) { | ||
273 | if(x>0) { | ||
274 | return +1.0; | ||
275 | } else if (x<0) { | ||
276 | return -1.0; | ||
277 | } else { | ||
278 | return 0.0; | ||
279 | } | ||
280 | } | ||
281 | |||
282 | inline gsl_complex complex_abs(gsl_complex z) { | ||
283 | gsl_complex r; | ||
284 | r.dat[0] = gsl_complex_abs(z); | ||
285 | r.dat[1] = 0; | ||
286 | return r; | ||
287 | } | ||
288 | |||
289 | inline gsl_complex complex_signum(gsl_complex z) { | ||
290 | gsl_complex r; | ||
291 | double mag; | ||
292 | if (z.dat[0] == 0 && z.dat[1] == 0) { | ||
293 | r.dat[0] = 0; | ||
294 | r.dat[1] = 0; | ||
295 | } else { | ||
296 | mag = gsl_complex_abs(z); | ||
297 | r.dat[0] = z.dat[0]/mag; | ||
298 | r.dat[1] = z.dat[1]/mag; | ||
299 | } | ||
300 | return r; | ||
301 | } | ||
302 | |||
303 | #define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK } | ||
304 | #define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK } | ||
305 | int mapR(int code, KRVEC(x), RVEC(r)) { | ||
306 | int k; | ||
307 | REQUIRES(xn == rn,BAD_SIZE); | ||
308 | DEBUGMSG("mapR"); | ||
309 | switch (code) { | ||
310 | OP(0,sin) | ||
311 | OP(1,cos) | ||
312 | OP(2,tan) | ||
313 | OP(3,fabs) | ||
314 | OP(4,asin) | ||
315 | OP(5,acos) | ||
316 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
317 | OP(7,sinh) | ||
318 | OP(8,cosh) | ||
319 | OP(9,tanh) | ||
320 | OP(10,gsl_asinh) | ||
321 | OP(11,gsl_acosh) | ||
322 | OP(12,gsl_atanh) | ||
323 | OP(13,exp) | ||
324 | OP(14,log) | ||
325 | OP(15,sign) | ||
326 | OP(16,sqrt) | ||
327 | default: ERROR(BAD_CODE); | ||
328 | } | ||
329 | } | ||
330 | |||
331 | int mapF(int code, KFVEC(x), FVEC(r)) { | ||
332 | int k; | ||
333 | REQUIRES(xn == rn,BAD_SIZE); | ||
334 | DEBUGMSG("mapF"); | ||
335 | switch (code) { | ||
336 | OP(0,sin) | ||
337 | OP(1,cos) | ||
338 | OP(2,tan) | ||
339 | OP(3,fabs) | ||
340 | OP(4,asin) | ||
341 | OP(5,acos) | ||
342 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
343 | OP(7,sinh) | ||
344 | OP(8,cosh) | ||
345 | OP(9,tanh) | ||
346 | OP(10,gsl_asinh) | ||
347 | OP(11,gsl_acosh) | ||
348 | OP(12,gsl_atanh) | ||
349 | OP(13,exp) | ||
350 | OP(14,log) | ||
351 | OP(15,sign) | ||
352 | OP(16,sqrt) | ||
353 | default: ERROR(BAD_CODE); | ||
354 | } | ||
355 | } | ||
356 | |||
357 | |||
358 | int mapCAux(int code, KGCVEC(x), GCVEC(r)) { | ||
359 | int k; | ||
360 | REQUIRES(xn == rn,BAD_SIZE); | ||
361 | DEBUGMSG("mapC"); | ||
362 | switch (code) { | ||
363 | OP(0,gsl_complex_sin) | ||
364 | OP(1,gsl_complex_cos) | ||
365 | OP(2,gsl_complex_tan) | ||
366 | OP(3,complex_abs) | ||
367 | OP(4,gsl_complex_arcsin) | ||
368 | OP(5,gsl_complex_arccos) | ||
369 | OP(6,gsl_complex_arctan) | ||
370 | OP(7,gsl_complex_sinh) | ||
371 | OP(8,gsl_complex_cosh) | ||
372 | OP(9,gsl_complex_tanh) | ||
373 | OP(10,gsl_complex_arcsinh) | ||
374 | OP(11,gsl_complex_arccosh) | ||
375 | OP(12,gsl_complex_arctanh) | ||
376 | OP(13,gsl_complex_exp) | ||
377 | OP(14,gsl_complex_log) | ||
378 | OP(15,complex_signum) | ||
379 | OP(16,gsl_complex_sqrt) | ||
380 | |||
381 | // gsl_complex_arg | ||
382 | // gsl_complex_abs | ||
383 | default: ERROR(BAD_CODE); | ||
384 | } | ||
385 | } | ||
386 | |||
387 | int mapC(int code, KCVEC(x), CVEC(r)) { | ||
388 | return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp); | ||
389 | } | ||
390 | |||
391 | |||
392 | gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a) | ||
393 | { | ||
394 | gsl_complex c; | ||
395 | gsl_complex r; | ||
396 | |||
397 | gsl_complex_float float_r; | ||
398 | |||
399 | c.dat[0] = a.dat[0]; | ||
400 | c.dat[1] = a.dat[1]; | ||
401 | |||
402 | r = (*cf)(c); | ||
403 | |||
404 | float_r.dat[0] = r.dat[0]; | ||
405 | float_r.dat[1] = r.dat[1]; | ||
406 | |||
407 | return float_r; | ||
408 | } | ||
409 | |||
410 | gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), | ||
411 | gsl_complex_float a,gsl_complex_float b) | ||
412 | { | ||
413 | gsl_complex c1; | ||
414 | gsl_complex c2; | ||
415 | gsl_complex r; | ||
416 | |||
417 | gsl_complex_float float_r; | ||
418 | |||
419 | c1.dat[0] = a.dat[0]; | ||
420 | c1.dat[1] = a.dat[1]; | ||
421 | |||
422 | c2.dat[0] = b.dat[0]; | ||
423 | c2.dat[1] = b.dat[1]; | ||
424 | |||
425 | r = (*cf)(c1,c2); | ||
426 | |||
427 | float_r.dat[0] = r.dat[0]; | ||
428 | float_r.dat[1] = r.dat[1]; | ||
429 | |||
430 | return float_r; | ||
431 | } | ||
432 | |||
433 | #define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK } | ||
434 | #define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK } | ||
435 | int mapQAux(int code, KGQVEC(x), GQVEC(r)) { | ||
436 | int k; | ||
437 | REQUIRES(xn == rn,BAD_SIZE); | ||
438 | DEBUGMSG("mapQ"); | ||
439 | switch (code) { | ||
440 | OPC(0,gsl_complex_sin) | ||
441 | OPC(1,gsl_complex_cos) | ||
442 | OPC(2,gsl_complex_tan) | ||
443 | OPC(3,complex_abs) | ||
444 | OPC(4,gsl_complex_arcsin) | ||
445 | OPC(5,gsl_complex_arccos) | ||
446 | OPC(6,gsl_complex_arctan) | ||
447 | OPC(7,gsl_complex_sinh) | ||
448 | OPC(8,gsl_complex_cosh) | ||
449 | OPC(9,gsl_complex_tanh) | ||
450 | OPC(10,gsl_complex_arcsinh) | ||
451 | OPC(11,gsl_complex_arccosh) | ||
452 | OPC(12,gsl_complex_arctanh) | ||
453 | OPC(13,gsl_complex_exp) | ||
454 | OPC(14,gsl_complex_log) | ||
455 | OPC(15,complex_signum) | ||
456 | OPC(16,gsl_complex_sqrt) | ||
457 | |||
458 | // gsl_complex_arg | ||
459 | // gsl_complex_abs | ||
460 | default: ERROR(BAD_CODE); | ||
461 | } | ||
462 | } | ||
463 | |||
464 | int mapQ(int code, KQVEC(x), QVEC(r)) { | ||
465 | return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp); | ||
466 | } | ||
467 | |||
468 | |||
469 | int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { | ||
470 | int k; | ||
471 | double val = *pval; | ||
472 | REQUIRES(xn == rn,BAD_SIZE); | ||
473 | DEBUGMSG("mapValR"); | ||
474 | switch (code) { | ||
475 | OPV(0,val*xp[k]) | ||
476 | OPV(1,val/xp[k]) | ||
477 | OPV(2,val+xp[k]) | ||
478 | OPV(3,val-xp[k]) | ||
479 | OPV(4,pow(val,xp[k])) | ||
480 | OPV(5,pow(xp[k],val)) | ||
481 | default: ERROR(BAD_CODE); | ||
482 | } | ||
483 | } | ||
484 | |||
485 | int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) { | ||
486 | int k; | ||
487 | float val = *pval; | ||
488 | REQUIRES(xn == rn,BAD_SIZE); | ||
489 | DEBUGMSG("mapValF"); | ||
490 | switch (code) { | ||
491 | OPV(0,val*xp[k]) | ||
492 | OPV(1,val/xp[k]) | ||
493 | OPV(2,val+xp[k]) | ||
494 | OPV(3,val-xp[k]) | ||
495 | OPV(4,pow(val,xp[k])) | ||
496 | OPV(5,pow(xp[k],val)) | ||
497 | default: ERROR(BAD_CODE); | ||
498 | } | ||
499 | } | ||
500 | |||
501 | int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) { | ||
502 | int k; | ||
503 | gsl_complex val = *pval; | ||
504 | REQUIRES(xn == rn,BAD_SIZE); | ||
505 | DEBUGMSG("mapValC"); | ||
506 | switch (code) { | ||
507 | OPV(0,gsl_complex_mul(val,xp[k])) | ||
508 | OPV(1,gsl_complex_div(val,xp[k])) | ||
509 | OPV(2,gsl_complex_add(val,xp[k])) | ||
510 | OPV(3,gsl_complex_sub(val,xp[k])) | ||
511 | OPV(4,gsl_complex_pow(val,xp[k])) | ||
512 | OPV(5,gsl_complex_pow(xp[k],val)) | ||
513 | default: ERROR(BAD_CODE); | ||
514 | } | ||
515 | } | ||
516 | |||
517 | int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) { | ||
518 | return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp); | ||
519 | } | ||
520 | |||
521 | |||
522 | int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) { | ||
523 | int k; | ||
524 | gsl_complex_float val = *pval; | ||
525 | REQUIRES(xn == rn,BAD_SIZE); | ||
526 | DEBUGMSG("mapValQ"); | ||
527 | switch (code) { | ||
528 | OPCA(0,gsl_complex_mul,val,xp[k]) | ||
529 | OPCA(1,gsl_complex_div,val,xp[k]) | ||
530 | OPCA(2,gsl_complex_add,val,xp[k]) | ||
531 | OPCA(3,gsl_complex_sub,val,xp[k]) | ||
532 | OPCA(4,gsl_complex_pow,val,xp[k]) | ||
533 | OPCA(5,gsl_complex_pow,xp[k],val) | ||
534 | default: ERROR(BAD_CODE); | ||
535 | } | ||
536 | } | ||
537 | |||
538 | int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) { | ||
539 | return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp); | ||
540 | } | ||
541 | |||
542 | |||
543 | #define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK } | ||
544 | #define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK } | ||
545 | int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) { | ||
546 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
547 | int k; | ||
548 | switch(code) { | ||
549 | OPZE(4,"zipR Pow",pow) | ||
550 | OPZE(5,"zipR ATan2",atan2) | ||
551 | } | ||
552 | KDVVIEW(a); | ||
553 | KDVVIEW(b); | ||
554 | DVVIEW(r); | ||
555 | gsl_vector_memcpy(V(r),V(a)); | ||
556 | int res; | ||
557 | switch(code) { | ||
558 | OPZV(0,"zipR Add",gsl_vector_add) | ||
559 | OPZV(1,"zipR Sub",gsl_vector_sub) | ||
560 | OPZV(2,"zipR Mul",gsl_vector_mul) | ||
561 | OPZV(3,"zipR Div",gsl_vector_div) | ||
562 | default: ERROR(BAD_CODE); | ||
563 | } | ||
564 | } | ||
565 | |||
566 | |||
567 | int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) { | ||
568 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
569 | int k; | ||
570 | switch(code) { | ||
571 | OPZE(4,"zipF Pow",pow) | ||
572 | OPZE(5,"zipF ATan2",atan2) | ||
573 | } | ||
574 | KFVVIEW(a); | ||
575 | KFVVIEW(b); | ||
576 | FVVIEW(r); | ||
577 | gsl_vector_float_memcpy(V(r),V(a)); | ||
578 | int res; | ||
579 | switch(code) { | ||
580 | OPZV(0,"zipF Add",gsl_vector_float_add) | ||
581 | OPZV(1,"zipF Sub",gsl_vector_float_sub) | ||
582 | OPZV(2,"zipF Mul",gsl_vector_float_mul) | ||
583 | OPZV(3,"zipF Div",gsl_vector_float_div) | ||
584 | default: ERROR(BAD_CODE); | ||
585 | } | ||
586 | } | ||
587 | |||
588 | |||
589 | int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) { | ||
590 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
591 | int k; | ||
592 | switch(code) { | ||
593 | OPZE(0,"zipC Add",gsl_complex_add) | ||
594 | OPZE(1,"zipC Sub",gsl_complex_sub) | ||
595 | OPZE(2,"zipC Mul",gsl_complex_mul) | ||
596 | OPZE(3,"zipC Div",gsl_complex_div) | ||
597 | OPZE(4,"zipC Pow",gsl_complex_pow) | ||
598 | //OPZE(5,"zipR ATan2",atan2) | ||
599 | } | ||
600 | //KCVVIEW(a); | ||
601 | //KCVVIEW(b); | ||
602 | //CVVIEW(r); | ||
603 | //gsl_vector_memcpy(V(r),V(a)); | ||
604 | //int res; | ||
605 | switch(code) { | ||
606 | default: ERROR(BAD_CODE); | ||
607 | } | ||
608 | } | ||
609 | |||
610 | |||
611 | int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { | ||
612 | return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp); | ||
613 | } | ||
614 | |||
615 | |||
616 | #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 } | ||
617 | int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) { | ||
618 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
619 | int k; | ||
620 | switch(code) { | ||
621 | OPCZE(0,"zipQ Add",gsl_complex_add) | ||
622 | OPCZE(1,"zipQ Sub",gsl_complex_sub) | ||
623 | OPCZE(2,"zipQ Mul",gsl_complex_mul) | ||
624 | OPCZE(3,"zipQ Div",gsl_complex_div) | ||
625 | OPCZE(4,"zipQ Pow",gsl_complex_pow) | ||
626 | //OPZE(5,"zipR ATan2",atan2) | ||
627 | } | ||
628 | //KCVVIEW(a); | ||
629 | //KCVVIEW(b); | ||
630 | //CVVIEW(r); | ||
631 | //gsl_vector_memcpy(V(r),V(a)); | ||
632 | //int res; | ||
633 | switch(code) { | ||
634 | default: ERROR(BAD_CODE); | ||
635 | } | ||
636 | } | ||
637 | |||
638 | |||
639 | int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) { | ||
640 | return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp); | ||
641 | } | ||
642 | |||