diff options
-rw-r--r-- | packages/base/src/C/vector-aux.c | 104 | ||||
-rw-r--r-- | packages/hmatrix/hmatrix.cabal | 3 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/Vector.hs | 58 | ||||
-rw-r--r-- | packages/hmatrix/src/Numeric/GSL/gsl-vector.c | 349 |
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 | ||
477 | inline 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 } | ||
500 | int mapQ(int code, KQVEC(x), QVEC(r)) { | 466 | int 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 } | |
571 | int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) { | 537 | int 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 | ||
588 | int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) { | 556 | int 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 } | ||
641 | int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { | 612 | int 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 | ||
661 | int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) { | 634 | int 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 | ||
25 | import Data.Packed | 25 | import Data.Packed |
26 | import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) | 26 | import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) |
27 | import Numeric.Vectorized( | 27 | import 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 | ||
36 | import Data.Complex | 29 | import Data.Complex |
37 | import Foreign.Marshal.Alloc(free) | 30 | import Foreign.Marshal.Alloc(free) |
38 | import Foreign.Marshal.Array(newArray) | ||
39 | import Foreign.Ptr(Ptr) | 31 | import Foreign.Ptr(Ptr) |
40 | import Foreign.C.Types | 32 | import Foreign.C.Types |
41 | import Foreign.C.String(newCString) | 33 | import Foreign.C.String(newCString) |
42 | import System.IO.Unsafe(unsafePerformIO) | 34 | import System.IO.Unsafe(unsafePerformIO) |
43 | import Control.Monad(when) | ||
44 | 35 | ||
45 | fromei x = fromIntegral (fromEnum x) :: CInt | 36 | fromei x = fromIntegral (fromEnum x) :: CInt |
46 | 37 | ||
47 | ------------------------------------------------------------------ | ||
48 | |||
49 | vectorMapAux fun code v = unsafePerformIO $ do | ||
50 | r <- createVector (dim v) | ||
51 | app2 (fun (fromei code)) vec v vec r "vectorMapAux" | ||
52 | return r | ||
53 | |||
54 | vectorMapValAux 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 | |||
61 | vectorZipAux 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 | ||
69 | vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double) | ||
70 | vectorMapValC = vectorMapValAux c_vectorMapValC | ||
71 | |||
72 | foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV | ||
73 | |||
74 | -- | map of complex vectors with given function | ||
75 | vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float) | ||
76 | vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper) | ||
77 | |||
78 | foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV | ||
79 | |||
80 | ------------------------------------------------------------------- | ||
81 | |||
82 | -- | elementwise operation on complex vectors | ||
83 | vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) | ||
84 | vectorZipC = vectorZipAux c_vectorZipC | ||
85 | |||
86 | foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV | ||
87 | |||
88 | -- | elementwise operation on complex vectors | ||
89 | vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) | ||
90 | vectorZipQ = vectorZipAux c_vectorZipQ | ||
91 | |||
92 | foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV | ||
93 | |||
94 | ----------------------------------------------------------------------- | 38 | ----------------------------------------------------------------------- |
95 | 39 | ||
96 | data RandDist = Uniform -- ^ uniform distribution in [0,1) | 40 | data 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 | |||
92 | inline 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 | |||
102 | inline 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 | |||
112 | inline 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 | |||
119 | inline 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 | |||
137 | int 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 | |||
166 | int mapC(int code, KCVEC(x), CVEC(r)) { | ||
167 | return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp); | ||
168 | } | ||
169 | |||
170 | |||
171 | gsl_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 | |||
189 | gsl_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 } | ||
214 | int 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 | |||
243 | int 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 | |||
249 | int 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 | |||
265 | int 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 | |||
270 | int 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 | |||
286 | int 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 | |||
296 | int 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 | |||
318 | int 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 } | ||
324 | int 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 | |||
346 | int 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 | |||