diff options
Diffstat (limited to 'lib/Numeric/GSL')
-rw-r--r-- | lib/Numeric/GSL/Vector.hs | 32 | ||||
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 114 |
2 files changed, 142 insertions, 4 deletions
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs index 28c3b06..d09323b 100644 --- a/lib/Numeric/GSL/Vector.hs +++ b/lib/Numeric/GSL/Vector.hs | |||
@@ -14,10 +14,10 @@ | |||
14 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
15 | 15 | ||
16 | module Numeric.GSL.Vector ( | 16 | module Numeric.GSL.Vector ( |
17 | FunCodeS(..), toScalarR, | 17 | FunCodeS(..), toScalarR, toScalarF, |
18 | FunCodeV(..), vectorMapR, vectorMapC, | 18 | FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, |
19 | FunCodeSV(..), vectorMapValR, vectorMapValC, | 19 | FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, |
20 | FunCodeVV(..), vectorZipR, vectorZipC, | 20 | FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, |
21 | RandDist(..), randomVector | 21 | RandDist(..), randomVector |
22 | ) where | 22 | ) where |
23 | 23 | ||
@@ -106,6 +106,12 @@ toScalarR oper = toScalarAux c_toScalarR (fromei oper) | |||
106 | 106 | ||
107 | foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV | 107 | foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV |
108 | 108 | ||
109 | -- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. | ||
110 | toScalarF :: FunCodeS -> Vector Float -> Float | ||
111 | toScalarF oper = toScalarAux c_toScalarF (fromei oper) | ||
112 | |||
113 | foreign import ccall safe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF | ||
114 | |||
109 | ------------------------------------------------------------------ | 115 | ------------------------------------------------------------------ |
110 | 116 | ||
111 | -- | map of real vectors with given function | 117 | -- | map of real vectors with given function |
@@ -120,6 +126,12 @@ vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper) | |||
120 | 126 | ||
121 | foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV | 127 | foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV |
122 | 128 | ||
129 | -- | map of real vectors with given function | ||
130 | vectorMapF :: FunCodeV -> Vector Float -> Vector Float | ||
131 | vectorMapF = vectorMapAux c_vectorMapF | ||
132 | |||
133 | foreign import ccall safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF | ||
134 | |||
123 | ------------------------------------------------------------------- | 135 | ------------------------------------------------------------------- |
124 | 136 | ||
125 | -- | map of real vectors with given function | 137 | -- | map of real vectors with given function |
@@ -134,6 +146,12 @@ vectorMapValC = vectorMapValAux c_vectorMapValC | |||
134 | 146 | ||
135 | foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV | 147 | foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV |
136 | 148 | ||
149 | -- | map of real vectors with given function | ||
150 | vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float | ||
151 | vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper) | ||
152 | |||
153 | foreign import ccall safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF | ||
154 | |||
137 | ------------------------------------------------------------------- | 155 | ------------------------------------------------------------------- |
138 | 156 | ||
139 | -- | elementwise operation on real vectors | 157 | -- | elementwise operation on real vectors |
@@ -148,6 +166,12 @@ vectorZipC = vectorZipAux c_vectorZipC | |||
148 | 166 | ||
149 | foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV | 167 | foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV |
150 | 168 | ||
169 | -- | elementwise operation on real vectors | ||
170 | vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float | ||
171 | vectorZipF = vectorZipAux c_vectorZipF | ||
172 | |||
173 | foreign import ccall safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF | ||
174 | |||
151 | ----------------------------------------------------------------------- | 175 | ----------------------------------------------------------------------- |
152 | 176 | ||
153 | data RandDist = Uniform -- ^ uniform distribution in [0,1) | 177 | data RandDist = Uniform -- ^ uniform distribution in [0,1) |
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index f5a95f4..6bb16f0 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -10,6 +10,16 @@ | |||
10 | #define KCVEC(A) int A##n, const 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 | 11 | #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p |
12 | 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 | |||
13 | #include <gsl/gsl_blas.h> | 23 | #include <gsl/gsl_blas.h> |
14 | #include <gsl/gsl_math.h> | 24 | #include <gsl/gsl_math.h> |
15 | #include <gsl/gsl_errno.h> | 25 | #include <gsl/gsl_errno.h> |
@@ -64,12 +74,24 @@ | |||
64 | #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) | 74 | #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) | 75 | #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) |
66 | 76 | ||
77 | #define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) | ||
78 | #define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) | ||
79 | #define QVVIEW(A) gsl_vector_float_complex_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) | ||
80 | #define QMVIEW(A) gsl_matrix_float_complex_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) | ||
81 | #define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) | ||
82 | #define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) | ||
83 | #define KQVVIEW(A) gsl_vector_float_complex_const_view A = gsl_vector_float_complex_const_view_array((float*)A##p,A##n) | ||
84 | #define KQMVIEW(A) gsl_matrix_float_complex_const_view A = gsl_matrix_float_complex_const_view_array((float*)A##p,A##r,A##c) | ||
85 | |||
67 | #define V(a) (&a.vector) | 86 | #define V(a) (&a.vector) |
68 | #define M(a) (&a.matrix) | 87 | #define M(a) (&a.matrix) |
69 | 88 | ||
70 | #define GCVEC(A) int A##n, gsl_complex*A##p | 89 | #define GCVEC(A) int A##n, gsl_complex*A##p |
71 | #define KGCVEC(A) int A##n, const gsl_complex*A##p | 90 | #define KGCVEC(A) int A##n, const gsl_complex*A##p |
72 | 91 | ||
92 | #define GQVEC(A) int A##n, gsl_complex_float*A##p | ||
93 | #define KGQVEC(A) int A##n, const gsl_complex_float*A##p | ||
94 | |||
73 | #define BAD_SIZE 2000 | 95 | #define BAD_SIZE 2000 |
74 | #define BAD_CODE 2001 | 96 | #define BAD_CODE 2001 |
75 | #define MEM 2002 | 97 | #define MEM 2002 |
@@ -99,6 +121,24 @@ int toScalarR(int code, KRVEC(x), RVEC(r)) { | |||
99 | OK | 121 | OK |
100 | } | 122 | } |
101 | 123 | ||
124 | int toScalarF(int code, KFVEC(x), FVEC(r)) { | ||
125 | REQUIRES(rn==1,BAD_SIZE); | ||
126 | DEBUGMSG("toScalarF"); | ||
127 | KFVVIEW(x); | ||
128 | float res; | ||
129 | switch(code) { | ||
130 | case 0: { res = gsl_blas_snrm2(V(x)); break; } | ||
131 | case 1: { res = gsl_blas_sasum(V(x)); break; } | ||
132 | case 2: { res = gsl_vector_float_max_index(V(x)); break; } | ||
133 | case 3: { res = gsl_vector_float_max(V(x)); break; } | ||
134 | case 4: { res = gsl_vector_float_min_index(V(x)); break; } | ||
135 | case 5: { res = gsl_vector_float_min(V(x)); break; } | ||
136 | default: ERROR(BAD_CODE); | ||
137 | } | ||
138 | rp[0] = res; | ||
139 | OK | ||
140 | } | ||
141 | |||
102 | 142 | ||
103 | inline double sign(double x) { | 143 | inline double sign(double x) { |
104 | if(x>0) { | 144 | if(x>0) { |
@@ -110,6 +150,16 @@ inline double sign(double x) { | |||
110 | } | 150 | } |
111 | } | 151 | } |
112 | 152 | ||
153 | inline float float_sign(float x) { | ||
154 | if(x>0) { | ||
155 | return +1.0; | ||
156 | } else if (x<0) { | ||
157 | return -1.0; | ||
158 | } else { | ||
159 | return 0.0; | ||
160 | } | ||
161 | } | ||
162 | |||
113 | inline gsl_complex complex_abs(gsl_complex z) { | 163 | inline gsl_complex complex_abs(gsl_complex z) { |
114 | gsl_complex r; | 164 | gsl_complex r; |
115 | r.dat[0] = gsl_complex_abs(z); | 165 | r.dat[0] = gsl_complex_abs(z); |
@@ -159,6 +209,32 @@ int mapR(int code, KRVEC(x), RVEC(r)) { | |||
159 | } | 209 | } |
160 | } | 210 | } |
161 | 211 | ||
212 | int mapF(int code, KFVEC(x), FVEC(r)) { | ||
213 | int k; | ||
214 | REQUIRES(xn == rn,BAD_SIZE); | ||
215 | DEBUGMSG("mapF"); | ||
216 | switch (code) { | ||
217 | OP(0,sin) | ||
218 | OP(1,cos) | ||
219 | OP(2,tan) | ||
220 | OP(3,fabs) | ||
221 | OP(4,asin) | ||
222 | OP(5,acos) | ||
223 | OP(6,atan) /* atan2 mediante vectorZip */ | ||
224 | OP(7,sinh) | ||
225 | OP(8,cosh) | ||
226 | OP(9,tanh) | ||
227 | OP(10,gsl_asinh) | ||
228 | OP(11,gsl_acosh) | ||
229 | OP(12,gsl_atanh) | ||
230 | OP(13,exp) | ||
231 | OP(14,log) | ||
232 | OP(15,sign) | ||
233 | OP(16,sqrt) | ||
234 | default: ERROR(BAD_CODE); | ||
235 | } | ||
236 | } | ||
237 | |||
162 | 238 | ||
163 | int mapCAux(int code, KGCVEC(x), GCVEC(r)) { | 239 | int mapCAux(int code, KGCVEC(x), GCVEC(r)) { |
164 | int k; | 240 | int k; |
@@ -210,6 +286,22 @@ int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { | |||
210 | } | 286 | } |
211 | } | 287 | } |
212 | 288 | ||
289 | int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) { | ||
290 | int k; | ||
291 | float val = *pval; | ||
292 | REQUIRES(xn == rn,BAD_SIZE); | ||
293 | DEBUGMSG("mapValF"); | ||
294 | switch (code) { | ||
295 | OPV(0,val*xp[k]) | ||
296 | OPV(1,val/xp[k]) | ||
297 | OPV(2,val+xp[k]) | ||
298 | OPV(3,val-xp[k]) | ||
299 | OPV(4,pow(val,xp[k])) | ||
300 | OPV(5,pow(xp[k],val)) | ||
301 | default: ERROR(BAD_CODE); | ||
302 | } | ||
303 | } | ||
304 | |||
213 | int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) { | 305 | int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) { |
214 | int k; | 306 | int k; |
215 | gsl_complex val = *pval; | 307 | gsl_complex val = *pval; |
@@ -255,6 +347,28 @@ int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) { | |||
255 | } | 347 | } |
256 | 348 | ||
257 | 349 | ||
350 | int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) { | ||
351 | REQUIRES(an == bn && an == rn, BAD_SIZE); | ||
352 | int k; | ||
353 | switch(code) { | ||
354 | OPZE(4,"zipF Pow",pow) | ||
355 | OPZE(5,"zipF ATan2",atan2) | ||
356 | } | ||
357 | KFVVIEW(a); | ||
358 | KFVVIEW(b); | ||
359 | FVVIEW(r); | ||
360 | gsl_vector_float_memcpy(V(r),V(a)); | ||
361 | int res; | ||
362 | switch(code) { | ||
363 | OPZV(0,"zipF Add",gsl_vector_float_add) | ||
364 | OPZV(1,"zipF Sub",gsl_vector_float_sub) | ||
365 | OPZV(2,"zipF Mul",gsl_vector_float_mul) | ||
366 | OPZV(3,"zipF Div",gsl_vector_float_div) | ||
367 | default: ERROR(BAD_CODE); | ||
368 | } | ||
369 | } | ||
370 | |||
371 | |||
258 | int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) { | 372 | int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) { |
259 | REQUIRES(an == bn && an == rn, BAD_SIZE); | 373 | REQUIRES(an == bn && an == rn, BAD_SIZE); |
260 | int k; | 374 | int k; |