summaryrefslogtreecommitdiff
path: root/packages/hmatrix
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix')
-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
3 files changed, 2 insertions, 408 deletions
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