summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-14 19:46:44 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-14 19:46:44 +0200
commit0f9575462eb37a7c9985583ca33ae315f6e6431d (patch)
treec8f4c6519c1a4b6969e4ae80fed9e410b852f7d7 /packages/hmatrix/src
parent494fe8f3e66f52244c4e0aedfd00eb15b9caceef (diff)
vectorized operations in base
Diffstat (limited to 'packages/hmatrix/src')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Vector.hs34
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-aux.c596
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-vector.c642
3 files changed, 642 insertions, 630 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
index 29c8bb7..3591289 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs
@@ -13,7 +13,6 @@
13module Numeric.GSL.Vector ( 13module Numeric.GSL.Vector (
14 sumF, sumR, sumQ, sumC, 14 sumF, sumR, sumQ, sumC,
15 prodF, prodR, prodQ, prodC, 15 prodF, prodR, prodQ, prodC,
16 dotF, dotR, dotQ, dotC,
17 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ, 16 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
18 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, 17 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
19 FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, 18 FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
@@ -148,39 +147,6 @@ foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV
148foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV 147foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV
149foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV 148foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV
150 149
151-- | dot product
152dotF :: Vector Float -> Vector Float -> Float
153dotF x y = unsafePerformIO $ do
154 r <- createVector 1
155 app3 c_dotF vec x vec y vec r "dotF"
156 return $ r @> 0
157
158-- | dot product
159dotR :: Vector Double -> Vector Double -> Double
160dotR x y = unsafePerformIO $ do
161 r <- createVector 1
162 app3 c_dotR vec x vec y vec r "dotR"
163 return $ r @> 0
164
165-- | dot product
166dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float
167dotQ x y = unsafePerformIO $ do
168 r <- createVector 1
169 app3 c_dotQ vec x vec y vec r "dotQ"
170 return $ r @> 0
171
172-- | dot product
173dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double
174dotC x y = unsafePerformIO $ do
175 r <- createVector 1
176 app3 c_dotC vec x vec y vec r "dotC"
177 return $ r @> 0
178
179foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF
180foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV
181foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV
182foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV
183
184------------------------------------------------------------------ 150------------------------------------------------------------------
185 151
186toScalarAux fun code v = unsafePerformIO $ do 152toScalarAux fun code v = unsafePerformIO $ do
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
index 410d157..5da94ca 100644
--- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
@@ -104,602 +104,6 @@ void no_abort_on_error() {
104} 104}
105 105
106 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 107
704int fft(int code, KCVEC(X), CVEC(R)) { 108int fft(int code, KCVEC(X), CVEC(R)) {
705 REQUIRES(Xn == Rn,BAD_SIZE); 109 REQUIRES(Xn == Rn,BAD_SIZE);
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
91int 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
101int 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
111int 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
126int 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
141int 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
151int 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
161int 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
178int 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
196int 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
214int 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
233int 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
247int 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
262inline 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
272inline 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
282inline 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
289inline 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 }
305int 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
331int 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
358int 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
387int mapC(int code, KCVEC(x), CVEC(r)) {
388 return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
389}
390
391
392gsl_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
410gsl_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 }
435int 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
464int mapQ(int code, KQVEC(x), QVEC(r)) {
465 return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
466}
467
468
469int 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
485int 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
501int 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
517int 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
522int 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
538int 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 }
545int 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
567int 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
589int 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
611int 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 }
617int 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
639int 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-09 15:34:17 +0000