summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-15 18:51:03 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-15 18:51:03 +0200
commit69d1fc1588532b48a946c1501f92ed56600baf4d (patch)
tree783de9d3dc07a5d1915b044a38982952125faf2e /packages
parentf34e120d08154dc97130eca2bc569d3cdd199270 (diff)
moved complex mapVal and zip
Diffstat (limited to 'packages')
-rw-r--r--packages/base/src/C/vector-aux.c104
-rw-r--r--packages/hmatrix/hmatrix.cabal3
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Vector.hs58
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-vector.c349
4 files changed, 42 insertions, 472 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index be2dc3a..7cdc750 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -29,19 +29,6 @@ typedef float complex TCF;
29 29
30#define CHECK(RES,CODE) MACRO(if(RES) return CODE;) 30#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
31 31
32#ifdef DBG
33#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
34#else
35#define DEBUGMAT(MSG,X)
36#endif
37
38#ifdef DBG
39#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
40#else
41#define DEBUGVEC(MSG,X)
42#endif
43
44
45#define BAD_SIZE 2000 32#define BAD_SIZE 2000
46#define BAD_CODE 2001 33#define BAD_CODE 2001
47#define MEM 2002 34#define MEM 2002
@@ -358,7 +345,7 @@ int mapR(int code, KDVEC(x), DVEC(r)) {
358 OP(3,fabs) 345 OP(3,fabs)
359 OP(4,asin) 346 OP(4,asin)
360 OP(5,acos) 347 OP(5,acos)
361 OP(6,atan) /* atan2 using vectorZip */ 348 OP(6,atan)
362 OP(7,sinh) 349 OP(7,sinh)
363 OP(8,cosh) 350 OP(8,cosh)
364 OP(9,tanh) 351 OP(9,tanh)
@@ -384,7 +371,7 @@ int mapF(int code, KFVEC(x), FVEC(r)) {
384 OP(3,fabs) 371 OP(3,fabs)
385 OP(4,asin) 372 OP(4,asin)
386 OP(5,acos) 373 OP(5,acos)
387 OP(6,atan) /* atan2 using vectorZip */ 374 OP(6,atan)
388 OP(7,sinh) 375 OP(7,sinh)
389 OP(8,cosh) 376 OP(8,cosh)
390 OP(9,tanh) 377 OP(9,tanh)
@@ -474,29 +461,8 @@ inline complex complex_f_math_fun(doublecomplex (*cf)(doublecomplex), complex a)
474 return float_r; 461 return float_r;
475} 462}
476 463
477inline complex complex_f_math_op(doublecomplex (*cf)(doublecomplex,doublecomplex),
478 complex a,complex b)
479{
480 doublecomplex c1,c2,r;
481
482 complex float_r;
483
484 c1.r = a.r;
485 c1.i = a.i;
486
487 c2.r = b.r;
488 c2.i = b.i;
489
490 r = (*cf)(c1,c2);
491
492 float_r.r = r.r;
493 float_r.i = r.i;
494
495 return float_r;
496}
497 464
498#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK } 465#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK }
499#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_op(&F,A,B); OK }
500int mapQ(int code, KQVEC(x), QVEC(r)) { 466int mapQ(int code, KQVEC(x), QVEC(r)) {
501 TCF* x2p = (TCF*)xp; 467 TCF* x2p = (TCF*)xp;
502 TCF* r2p = (TCF*)rp; 468 TCF* r2p = (TCF*)rp;
@@ -567,36 +533,40 @@ inline doublecomplex complex_add(doublecomplex a, doublecomplex b) {
567 return r; 533 return r;
568} 534}
569 535
570 536#define OPVb(C,E) case C: { for(k=0;k<xn;k++) r2p[k] = E; OK }
571int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) { 537int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) {
538 TCD* x2p = (TCD*)xp;
539 TCD* r2p = (TCD*)rp;
572 int k; 540 int k;
573 doublecomplex val = *pval; 541 TCD val = * (TCD*)pval;
574 REQUIRES(xn == rn,BAD_SIZE); 542 REQUIRES(xn == rn,BAD_SIZE);
575 DEBUGMSG("mapValC"); 543 DEBUGMSG("mapValC");
576 switch (code) { 544 switch (code) {
577// OPV(0,gsl_complex_mul(val,xp[k])) 545 OPVb(0,val*x2p[k])
578// OPV(1,gsl_complex_div(val,xp[k])) 546 OPVb(1,val/x2p[k])
579 OPV(2,complex_add(val,xp[k])) 547 OPVb(2,val+x2p[k])
580// OPV(3,gsl_complex_sub(val,xp[k])) 548 OPVb(3,val-x2p[k])
581// OPV(4,gsl_complex_pow(val,xp[k])) 549 OPVb(4,cpow(val,x2p[k]))
582// OPV(5,gsl_complex_pow(xp[k],val)) 550 OPVb(5,cpow(x2p[k],val))
583 default: ERROR(BAD_CODE); 551 default: ERROR(BAD_CODE);
584 } 552 }
585} 553}
586 554
587 555
588int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) { 556int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) {
557 TCF* x2p = (TCF*)xp;
558 TCF* r2p = (TCF*)rp;
589 int k; 559 int k;
590 complex val = *pval; 560 TCF val = *(TCF*)pval;
591 REQUIRES(xn == rn,BAD_SIZE); 561 REQUIRES(xn == rn,BAD_SIZE);
592 DEBUGMSG("mapValQ"); 562 DEBUGMSG("mapValQ");
593 switch (code) { 563 switch (code) {
594// OPCA(0,gsl_complex_mul,val,xp[k]) 564 OPVb(0,val*x2p[k])
595// OPCA(1,gsl_complex_div,val,xp[k]) 565 OPVb(1,val/x2p[k])
596 OPCA(2,complex_add,val,xp[k]) 566 OPVb(2,val+x2p[k])
597// OPCA(3,gsl_complex_sub,val,xp[k]) 567 OPVb(3,val-x2p[k])
598// OPCA(4,gsl_complex_pow,val,xp[k]) 568 OPVb(4,cpow(val,x2p[k]))
599// OPCA(5,gsl_complex_pow,xp[k],val) 569 OPVb(5,cpow(x2p[k],val))
600 default: ERROR(BAD_CODE); 570 default: ERROR(BAD_CODE);
601 } 571 }
602} 572}
@@ -637,17 +607,20 @@ REQUIRES(an == bn && an == rn, BAD_SIZE);
637 607
638 608
639 609
640 610#define OPZOb(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = a2p[k] O b2p[k]; OK }
611#define OPZEb(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) r2p[k] = E(a2p[k],b2p[k]); OK }
641int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { 612int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
642 REQUIRES(an == bn && an == rn, BAD_SIZE); 613 REQUIRES(an == bn && an == rn, BAD_SIZE);
614 TCD* a2p = (TCD*)ap;
615 TCD* b2p = (TCD*)bp;
616 TCD* r2p = (TCD*)rp;
643 int k; 617 int k;
644 switch(code) { 618 switch(code) {
645 OPZE(0,"zipC Add",complex_add) 619 OPZOb(0,"zipC Add",+)
646// OPZE(1,"zipC Sub",gsl_complex_sub) 620 OPZOb(1,"zipC Sub",-)
647// OPZE(2,"zipC Mul",gsl_complex_mul) 621 OPZOb(2,"zipC Mul",*)
648// OPZE(3,"zipC Div",gsl_complex_div) 622 OPZOb(3,"zipC Div",/)
649// OPZE(4,"zipC Pow",gsl_complex_pow) 623 OPZEb(4,"zipC Pow",cpow)
650// //OPZE(5,"zipR ATan2",atan2)
651 default: ERROR(BAD_CODE); 624 default: ERROR(BAD_CODE);
652 } 625 }
653} 626}
@@ -660,14 +633,17 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
660 633
661int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) { 634int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
662 REQUIRES(an == bn && an == rn, BAD_SIZE); 635 REQUIRES(an == bn && an == rn, BAD_SIZE);
636 TCF* a2p = (TCF*)ap;
637 TCF* b2p = (TCF*)bp;
638 TCF* r2p = (TCF*)rp;
639
663 int k; 640 int k;
664 switch(code) { 641 switch(code) {
665 OPCZE(0,"zipQ Add",complex_add) 642 OPZOb(0,"zipC Add",+)
666// OPCZE(1,"zipQ Sub",gsl_complex_sub) 643 OPZOb(1,"zipC Sub",-)
667// OPCZE(2,"zipQ Mul",gsl_complex_mul) 644 OPZOb(2,"zipC Mul",*)
668// OPCZE(3,"zipQ Div",gsl_complex_div) 645 OPZOb(3,"zipC Div",/)
669// OPCZE(4,"zipQ Pow",gsl_complex_pow) 646 OPZEb(4,"zipC Pow",cpowf)
670// //OPZE(5,"zipR ATan2",atan2)
671 default: ERROR(BAD_CODE); 647 default: ERROR(BAD_CODE);
672 } 648 }
673} 649}
diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal
index eaf4262..8496663 100644
--- a/packages/hmatrix/hmatrix.cabal
+++ b/packages/hmatrix/hmatrix.cabal
@@ -109,8 +109,7 @@ library
109 Numeric.LinearAlgebra.Util.Convolution 109 Numeric.LinearAlgebra.Util.Convolution
110 110
111 111
112 C-sources: src/Numeric/GSL/gsl-vector.c 112 C-sources: src/Numeric/GSL/gsl-aux.c
113 src/Numeric/GSL/gsl-aux.c
114 113
115 cc-options: -O4 -msse2 -Wall 114 cc-options: -O4 -msse2 -Wall
116 115
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
index 5c34f70..e133c2c 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs
@@ -24,73 +24,17 @@ module Numeric.GSL.Vector (
24 24
25import Data.Packed 25import Data.Packed
26import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 26import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
27import Numeric.Vectorized( 27import Numeric.Vectorized
28 sumF, sumR, sumQ, sumC,
29 prodF, prodR, prodQ, prodC,
30 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
31 FunCodeV(..), vectorMapR, vectorMapF, vectorMapC, vectorMapQ,
32 FunCodeSV(..), vectorMapValR, vectorMapValF,
33 FunCodeVV(..), vectorZipR, vectorZipF
34 )
35 28
36import Data.Complex 29import Data.Complex
37import Foreign.Marshal.Alloc(free) 30import Foreign.Marshal.Alloc(free)
38import Foreign.Marshal.Array(newArray)
39import Foreign.Ptr(Ptr) 31import Foreign.Ptr(Ptr)
40import Foreign.C.Types 32import Foreign.C.Types
41import Foreign.C.String(newCString) 33import Foreign.C.String(newCString)
42import System.IO.Unsafe(unsafePerformIO) 34import System.IO.Unsafe(unsafePerformIO)
43import Control.Monad(when)
44 35
45fromei x = fromIntegral (fromEnum x) :: CInt 36fromei x = fromIntegral (fromEnum x) :: CInt
46 37
47------------------------------------------------------------------
48
49vectorMapAux fun code v = unsafePerformIO $ do
50 r <- createVector (dim v)
51 app2 (fun (fromei code)) vec v vec r "vectorMapAux"
52 return r
53
54vectorMapValAux fun code val v = unsafePerformIO $ do
55 r <- createVector (dim v)
56 pval <- newArray [val]
57 app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux"
58 free pval
59 return r
60
61vectorZipAux fun code u v = unsafePerformIO $ do
62 r <- createVector (dim u)
63 when (dim u > 0) $ app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux"
64 return r
65
66---------------------------------------------------------------------
67
68-- | map of complex vectors with given function
69vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
70vectorMapValC = vectorMapValAux c_vectorMapValC
71
72foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV
73
74-- | map of complex vectors with given function
75vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)
76vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)
77
78foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV
79
80-------------------------------------------------------------------
81
82-- | elementwise operation on complex vectors
83vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
84vectorZipC = vectorZipAux c_vectorZipC
85
86foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV
87
88-- | elementwise operation on complex vectors
89vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)
90vectorZipQ = vectorZipAux c_vectorZipQ
91
92foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV
93
94----------------------------------------------------------------------- 38-----------------------------------------------------------------------
95 39
96data RandDist = Uniform -- ^ uniform distribution in [0,1) 40data RandDist = Uniform -- ^ uniform distribution in [0,1)
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-vector.c b/packages/hmatrix/src/Numeric/GSL/gsl-vector.c
deleted file mode 100644
index f00424a..0000000
--- a/packages/hmatrix/src/Numeric/GSL/gsl-vector.c
+++ /dev/null
@@ -1,349 +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_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
92inline double sign(double x) {
93 if(x>0) {
94 return +1.0;
95 } else if (x<0) {
96 return -1.0;
97 } else {
98 return 0.0;
99 }
100}
101
102inline float float_sign(float x) {
103 if(x>0) {
104 return +1.0;
105 } else if (x<0) {
106 return -1.0;
107 } else {
108 return 0.0;
109 }
110}
111
112inline gsl_complex complex_abs(gsl_complex z) {
113 gsl_complex r;
114 r.dat[0] = gsl_complex_abs(z);
115 r.dat[1] = 0;
116 return r;
117}
118
119inline gsl_complex complex_signum(gsl_complex z) {
120 gsl_complex r;
121 double mag;
122 if (z.dat[0] == 0 && z.dat[1] == 0) {
123 r.dat[0] = 0;
124 r.dat[1] = 0;
125 } else {
126 mag = gsl_complex_abs(z);
127 r.dat[0] = z.dat[0]/mag;
128 r.dat[1] = z.dat[1]/mag;
129 }
130 return r;
131}
132
133#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
134#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
135
136
137int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
138 int k;
139 REQUIRES(xn == rn,BAD_SIZE);
140 DEBUGMSG("mapC");
141 switch (code) {
142 OP(0,gsl_complex_sin)
143 OP(1,gsl_complex_cos)
144 OP(2,gsl_complex_tan)
145 OP(3,complex_abs)
146 OP(4,gsl_complex_arcsin)
147 OP(5,gsl_complex_arccos)
148 OP(6,gsl_complex_arctan)
149 OP(7,gsl_complex_sinh)
150 OP(8,gsl_complex_cosh)
151 OP(9,gsl_complex_tanh)
152 OP(10,gsl_complex_arcsinh)
153 OP(11,gsl_complex_arccosh)
154 OP(12,gsl_complex_arctanh)
155 OP(13,gsl_complex_exp)
156 OP(14,gsl_complex_log)
157 OP(15,complex_signum)
158 OP(16,gsl_complex_sqrt)
159
160 // gsl_complex_arg
161 // gsl_complex_abs
162 default: ERROR(BAD_CODE);
163 }
164}
165
166int mapC(int code, KCVEC(x), CVEC(r)) {
167 return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
168}
169
170
171gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
172{
173 gsl_complex c;
174 gsl_complex r;
175
176 gsl_complex_float float_r;
177
178 c.dat[0] = a.dat[0];
179 c.dat[1] = a.dat[1];
180
181 r = (*cf)(c);
182
183 float_r.dat[0] = r.dat[0];
184 float_r.dat[1] = r.dat[1];
185
186 return float_r;
187}
188
189gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex),
190 gsl_complex_float a,gsl_complex_float b)
191{
192 gsl_complex c1;
193 gsl_complex c2;
194 gsl_complex r;
195
196 gsl_complex_float float_r;
197
198 c1.dat[0] = a.dat[0];
199 c1.dat[1] = a.dat[1];
200
201 c2.dat[0] = b.dat[0];
202 c2.dat[1] = b.dat[1];
203
204 r = (*cf)(c1,c2);
205
206 float_r.dat[0] = r.dat[0];
207 float_r.dat[1] = r.dat[1];
208
209 return float_r;
210}
211
212#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
213#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
214int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
215 int k;
216 REQUIRES(xn == rn,BAD_SIZE);
217 DEBUGMSG("mapQ");
218 switch (code) {
219 OPC(0,gsl_complex_sin)
220 OPC(1,gsl_complex_cos)
221 OPC(2,gsl_complex_tan)
222 OPC(3,complex_abs)
223 OPC(4,gsl_complex_arcsin)
224 OPC(5,gsl_complex_arccos)
225 OPC(6,gsl_complex_arctan)
226 OPC(7,gsl_complex_sinh)
227 OPC(8,gsl_complex_cosh)
228 OPC(9,gsl_complex_tanh)
229 OPC(10,gsl_complex_arcsinh)
230 OPC(11,gsl_complex_arccosh)
231 OPC(12,gsl_complex_arctanh)
232 OPC(13,gsl_complex_exp)
233 OPC(14,gsl_complex_log)
234 OPC(15,complex_signum)
235 OPC(16,gsl_complex_sqrt)
236
237 // gsl_complex_arg
238 // gsl_complex_abs
239 default: ERROR(BAD_CODE);
240 }
241}
242
243int mapQ(int code, KQVEC(x), QVEC(r)) {
244 return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
245}
246
247
248
249int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
250 int k;
251 gsl_complex val = *pval;
252 REQUIRES(xn == rn,BAD_SIZE);
253 DEBUGMSG("mapValC");
254 switch (code) {
255 OPV(0,gsl_complex_mul(val,xp[k]))
256 OPV(1,gsl_complex_div(val,xp[k]))
257 OPV(2,gsl_complex_add(val,xp[k]))
258 OPV(3,gsl_complex_sub(val,xp[k]))
259 OPV(4,gsl_complex_pow(val,xp[k]))
260 OPV(5,gsl_complex_pow(xp[k],val))
261 default: ERROR(BAD_CODE);
262 }
263}
264
265int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
266 return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
267}
268
269
270int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
271 int k;
272 gsl_complex_float val = *pval;
273 REQUIRES(xn == rn,BAD_SIZE);
274 DEBUGMSG("mapValQ");
275 switch (code) {
276 OPCA(0,gsl_complex_mul,val,xp[k])
277 OPCA(1,gsl_complex_div,val,xp[k])
278 OPCA(2,gsl_complex_add,val,xp[k])
279 OPCA(3,gsl_complex_sub,val,xp[k])
280 OPCA(4,gsl_complex_pow,val,xp[k])
281 OPCA(5,gsl_complex_pow,xp[k],val)
282 default: ERROR(BAD_CODE);
283 }
284}
285
286int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
287 return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
288}
289
290
291#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
292#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
293
294
295
296int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
297 REQUIRES(an == bn && an == rn, BAD_SIZE);
298 int k;
299 switch(code) {
300 OPZE(0,"zipC Add",gsl_complex_add)
301 OPZE(1,"zipC Sub",gsl_complex_sub)
302 OPZE(2,"zipC Mul",gsl_complex_mul)
303 OPZE(3,"zipC Div",gsl_complex_div)
304 OPZE(4,"zipC Pow",gsl_complex_pow)
305 //OPZE(5,"zipR ATan2",atan2)
306 }
307 //KCVVIEW(a);
308 //KCVVIEW(b);
309 //CVVIEW(r);
310 //gsl_vector_memcpy(V(r),V(a));
311 //int res;
312 switch(code) {
313 default: ERROR(BAD_CODE);
314 }
315}
316
317
318int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
319 return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
320}
321
322
323#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 }
324int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
325 REQUIRES(an == bn && an == rn, BAD_SIZE);
326 int k;
327 switch(code) {
328 OPCZE(0,"zipQ Add",gsl_complex_add)
329 OPCZE(1,"zipQ Sub",gsl_complex_sub)
330 OPCZE(2,"zipQ Mul",gsl_complex_mul)
331 OPCZE(3,"zipQ Div",gsl_complex_div)
332 OPCZE(4,"zipQ Pow",gsl_complex_pow)
333 //OPZE(5,"zipR ATan2",atan2)
334 }
335 //KCVVIEW(a);
336 //KCVVIEW(b);
337 //CVVIEW(r);
338 //gsl_vector_memcpy(V(r),V(a));
339 //int res;
340 switch(code) {
341 default: ERROR(BAD_CODE);
342 }
343}
344
345
346int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
347 return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
348}
349