summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL/gsl-vector.c
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL/gsl-vector.c')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-vector.c349
1 files changed, 0 insertions, 349 deletions
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