summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/GSL/Vector.hs163
-rw-r--r--lib/GSL/gsl-aux.c721
-rw-r--r--lib/GSL/gsl-aux.h66
3 files changed, 950 insertions, 0 deletions
diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs
new file mode 100644
index 0000000..e538b7e
--- /dev/null
+++ b/lib/GSL/Vector.hs
@@ -0,0 +1,163 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : GSL.Vector
5-- Copyright : (c) Alberto Ruiz 2007
6-- License : GPL-style
7--
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional
10-- Portability : portable (uses FFI)
11--
12-- Vector operations
13--
14-----------------------------------------------------------------------------
15
16module GSL.Vector (
17 FunCodeS(..), toScalarR,
18 FunCodeV(..), vectorMapR, vectorMapC,
19 FunCodeSV(..), vectorMapValR, vectorMapValC,
20 FunCodeVV(..), vectorZipR, vectorZipC,
21 scale, addConstant
22) where
23
24import Data.Packed.Internal
25import Complex
26import Foreign
27
28data FunCodeV = Sin
29 | Cos
30 | Tan
31 | Abs
32 | ASin
33 | ACos
34 | ATan
35 | Sinh
36 | Cosh
37 | Tanh
38 | ASinh
39 | ACosh
40 | ATanh
41 | Exp
42 | Log
43 | Sign
44 | Sqrt
45 deriving Enum
46
47data FunCodeSV = Scale
48 | Recip
49 | AddConstant
50 | Negate
51 | PowSV
52 | PowVS
53 deriving Enum
54
55data FunCodeVV = Add
56 | Sub
57 | Mul
58 | Div
59 | Pow
60 | ATan2
61 deriving Enum
62
63data FunCodeS = Norm2
64 | AbsSum
65 | MaxIdx
66 | Max
67 | MinIdx
68 | Min
69 deriving Enum
70
71
72scale :: (Num a, Field a) => a -> Vector a -> Vector a
73scale x v | isReal baseOf v = scast $ vectorMapValR Scale (scast x) (scast v)
74 | isComp baseOf v = scast $ vectorMapValC Scale (scast x) (scast v)
75 | otherwise = fromList $ map (*x) $ toList v
76
77addConstant :: (Num a, Field a) => a -> Vector a -> Vector a
78addConstant x v | isReal baseOf v = scast $ vectorMapValR AddConstant (scast x) (scast v)
79 | isComp baseOf v = scast $ vectorMapValC AddConstant (scast x) (scast v)
80 | otherwise = fromList $ map (*x) $ toList v
81
82------------------------------------------------------------------
83
84toScalarAux fun code v = unsafePerformIO $ do
85 r <- createVector 1
86 fun (fromEnum code) // vec v // vec r // check "toScalarAux" [v]
87 return (r `at` 0)
88
89vectorMapAux fun code v = unsafePerformIO $ do
90 r <- createVector (dim v)
91 fun (fromEnum code) // vec v // vec r // check "vectorMapAux" [v]
92 return r
93
94vectorMapValAux fun code val v = unsafePerformIO $ do
95 r <- createVector (dim v)
96 pval <- newArray [val]
97 fun (fromEnum code) pval // vec v // vec r // check "vectorMapValAux" [v]
98 free pval
99 return r
100
101vectorZipAux fun code u v = unsafePerformIO $ do
102 r <- createVector (dim u)
103 fun (fromEnum code) // vec u // vec v // vec r // check "vectorZipAux" [u,v]
104 return r
105
106---------------------------------------------------------------------
107
108-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
109toScalarR :: FunCodeS -> Vector Double -> Double
110toScalarR oper = toScalarAux c_toScalarR (fromEnum oper)
111
112foreign import ccall safe "gsl-aux.h toScalarR"
113 c_toScalarR :: Int -> Double :> Double :> IO Int
114
115------------------------------------------------------------------
116
117-- | map of real vectors with given function
118vectorMapR :: FunCodeV -> Vector Double -> Vector Double
119vectorMapR = vectorMapAux c_vectorMapR
120
121foreign import ccall safe "gsl-aux.h mapR"
122 c_vectorMapR :: Int -> Double :> Double :> IO Int
123
124-- | map of complex vectors with given function
125vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
126vectorMapC oper = vectorMapAux c_vectorMapC (fromEnum oper)
127
128foreign import ccall safe "gsl-aux.h mapC"
129 c_vectorMapC :: Int -> Complex Double :> Complex Double :> IO Int
130
131-------------------------------------------------------------------
132
133-- | map of real vectors with given function
134vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
135vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromEnum oper)
136
137foreign import ccall safe "gsl-aux.h mapValR"
138 c_vectorMapValR :: Int -> Ptr Double -> Double :> Double :> IO Int
139
140-- | map of complex vectors with given function
141vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
142vectorMapValC = vectorMapValAux c_vectorMapValC
143
144foreign import ccall safe "gsl-aux.h mapValC"
145 c_vectorMapValC :: Int -> Ptr (Complex Double) -> Complex Double :> Complex Double :> IO Int
146
147-------------------------------------------------------------------
148
149-- | elementwise operation on real vectors
150vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
151vectorZipR = vectorZipAux c_vectorZipR
152
153foreign import ccall safe "gsl-aux.h zipR"
154 c_vectorZipR :: Int -> Double :> Double :> Double :> IO Int
155
156-- | elementwise operation on complex vectors
157vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
158vectorZipC = vectorZipAux c_vectorZipC
159
160foreign import ccall safe "gsl-aux.h zipC"
161 c_vectorZipC :: Int -> Complex Double :> Complex Double :> Complex Double :> IO Int
162
163
diff --git a/lib/GSL/gsl-aux.c b/lib/GSL/gsl-aux.c
new file mode 100644
index 0000000..fc95e6f
--- /dev/null
+++ b/lib/GSL/gsl-aux.c
@@ -0,0 +1,721 @@
1#include "gsl-aux.h"
2#include <gsl/gsl_blas.h>
3#include <gsl/gsl_linalg.h>
4#include <gsl/gsl_matrix.h>
5#include <gsl/gsl_math.h>
6#include <gsl/gsl_errno.h>
7#include <gsl/gsl_fft_complex.h>
8#include <gsl/gsl_eigen.h>
9#include <gsl/gsl_integration.h>
10#include <gsl/gsl_deriv.h>
11#include <gsl/gsl_poly.h>
12#include <gsl/gsl_multimin.h>
13#include <gsl/gsl_complex.h>
14#include <gsl/gsl_complex_math.h>
15#include <string.h>
16#include <stdio.h>
17
18#define MACRO(B) do {B} while (0)
19#define ERROR(CODE) MACRO(return CODE;)
20#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
21#define OK return 0;
22
23#define MIN(A,B) ((A)<(B)?(A):(B))
24#define MAX(A,B) ((A)>(B)?(A):(B))
25
26#ifdef DBG
27#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
28#else
29#define DEBUGMSG(M)
30#endif
31
32#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
33
34#ifdef DBG
35#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
36#else
37#define DEBUGMAT(MSG,X)
38#endif
39
40#ifdef DBG
41#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
42#else
43#define DEBUGVEC(MSG,X)
44#endif
45
46#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
47#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
48#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
49#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
50#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
51#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
52#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
53#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
54
55#define V(a) (&a.vector)
56#define M(a) (&a.matrix)
57
58#define GCVEC(A) int A##n, gsl_complex*A##p
59#define KGCVEC(A) int A##n, const gsl_complex*A##p
60
61#define BAD_SIZE 1000
62#define BAD_CODE 1001
63#define MEM 1002
64#define BAD_FILE 1003
65
66
67int toScalarR(int code, KRVEC(x), RVEC(r)) {
68 REQUIRES(rn==1,BAD_SIZE);
69 DEBUGMSG("toScalarR");
70 KDVVIEW(x);
71 double res;
72 switch(code) {
73 case 0: { res = gsl_blas_dnrm2(V(x)); break; }
74 case 1: { res = gsl_blas_dasum(V(x)); break; }
75 case 2: { res = gsl_vector_max_index(V(x)); break; }
76 case 3: { res = gsl_vector_max(V(x)); break; }
77 case 4: { res = gsl_vector_min_index(V(x)); break; }
78 case 5: { res = gsl_vector_min(V(x)); break; }
79 default: ERROR(BAD_CODE);
80 }
81 rp[0] = res;
82 OK
83}
84
85
86inline double sign(double x) {
87 if(x>0) {
88 return +1.0;
89 } else if (x<0) {
90 return -1.0;
91 } else {
92 return 0.0;
93 }
94}
95
96
97#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
98#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
99int mapR(int code, KRVEC(x), RVEC(r)) {
100 int k;
101 REQUIRES(xn == rn,BAD_SIZE);
102 DEBUGMSG("mapR");
103 switch (code) {
104 OP(0,sin)
105 OP(1,cos)
106 OP(2,tan)
107 OP(3,fabs)
108 OP(4,asin)
109 OP(5,acos)
110 OP(6,atan) /* atan2 mediante vectorZip */
111 OP(7,sinh)
112 OP(8,cosh)
113 OP(9,tanh)
114 OP(10,asinh)
115 OP(11,acosh)
116 OP(12,atanh)
117 OP(13,exp)
118 OP(14,log)
119 OP(15,sign)
120 OP(16,sqrt)
121 default: ERROR(BAD_CODE);
122 }
123}
124
125int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
126 int k;
127 REQUIRES(xn == rn,BAD_SIZE);
128 DEBUGMSG("mapC");
129 switch (code) {
130 OP(0,gsl_complex_sin)
131 OP(1,gsl_complex_cos)
132 OP(2,gsl_complex_tan)
133
134 OP(4,gsl_complex_arcsin)
135 OP(5,gsl_complex_arccos)
136 OP(6,gsl_complex_arctan)
137 OP(7,gsl_complex_sinh)
138 OP(8,gsl_complex_cosh)
139 OP(9,gsl_complex_tanh)
140 OP(10,gsl_complex_arcsinh)
141 OP(11,gsl_complex_arccosh)
142 OP(12,gsl_complex_arctanh)
143 OP(13,gsl_complex_exp)
144 OP(14,gsl_complex_log)
145
146 OP(16,gsl_complex_sqrt)
147
148 // gsl_complex_arg
149 // gsl_complex_abs
150 default: ERROR(BAD_CODE);
151 }
152}
153
154int mapC(int code, KCVEC(x), CVEC(r)) {
155 return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
156}
157
158
159/*
160int scaleR(double* alpha, KRVEC(x), RVEC(r)) {
161 REQUIRES(xn == rn,BAD_SIZE);
162 DEBUGMSG("scaleR");
163 KDVVIEW(x);
164 DVVIEW(r);
165 CHECK( gsl_vector_memcpy(V(r),V(x)) , MEM);
166 int res = gsl_vector_scale(V(r),*alpha);
167 CHECK(res,res);
168 OK
169}
170
171int scaleC(gsl_complex *alpha, KCVEC(x), CVEC(r)) {
172 REQUIRES(xn == rn,BAD_SIZE);
173 DEBUGMSG("scaleC");
174 //KCVVIEW(x);
175 CVVIEW(r);
176 gsl_vector_const_view vrx = gsl_vector_const_view_array((double*)xp,xn*2);
177 gsl_vector_view vrr = gsl_vector_view_array((double*)rp,rn*2);
178 CHECK(gsl_vector_memcpy(V(vrr),V(vrx)) , MEM);
179 gsl_blas_zscal(*alpha,V(r)); // void !
180 int res = 0;
181 CHECK(res,res);
182 OK
183}
184
185int addConstantR(double offs, KRVEC(x), RVEC(r)) {
186 REQUIRES(xn == rn,BAD_SIZE);
187 DEBUGMSG("addConstantR");
188 KDVVIEW(x);
189 DVVIEW(r);
190 CHECK(gsl_vector_memcpy(V(r),V(x)), MEM);
191 int res = gsl_vector_add_constant(V(r),offs);
192 CHECK(res,res);
193 OK
194}
195
196*/
197
198
199
200int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
201 int k;
202 double val = *pval;
203 REQUIRES(xn == rn,BAD_SIZE);
204 DEBUGMSG("mapValR");
205 switch (code) {
206 OPV(0,val*xp[k])
207 OPV(1,val/xp[k])
208 OPV(2,val+xp[k])
209 OPV(3,val-xp[k])
210 OPV(4,pow(val,xp[k]))
211 OPV(5,pow(xp[k],val))
212 default: ERROR(BAD_CODE);
213 }
214}
215
216int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
217 int k;
218 gsl_complex val = *pval;
219 REQUIRES(xn == rn,BAD_SIZE);
220 DEBUGMSG("mapValC");
221 switch (code) {
222 OPV(0,gsl_complex_mul(val,xp[k]))
223 OPV(1,gsl_complex_div(val,xp[k]))
224 OPV(2,gsl_complex_add(val,xp[k]))
225 OPV(3,gsl_complex_sub(val,xp[k]))
226 OPV(4,gsl_complex_pow(val,xp[k]))
227 OPV(5,gsl_complex_pow(val,xp[k]))
228 default: ERROR(BAD_CODE);
229 }
230}
231
232int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
233 return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
234}
235
236
237#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
238#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
239int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
240 REQUIRES(an == bn && an == rn, BAD_SIZE);
241 int k;
242 switch(code) {
243 OPZE(4,"zipR Pow",pow)
244 OPZE(5,"zipR ATan2",atan2)
245 }
246 KDVVIEW(a);
247 KDVVIEW(b);
248 DVVIEW(r);
249 gsl_vector_memcpy(V(r),V(a));
250 int res;
251 switch(code) {
252 OPZV(0,"zipR Add",gsl_vector_add)
253 OPZV(1,"zipR Sub",gsl_vector_sub)
254 OPZV(2,"zipR Mul",gsl_vector_mul)
255 OPZV(3,"zipR Div",gsl_vector_div)
256 default: ERROR(BAD_CODE);
257 }
258}
259
260
261int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
262 REQUIRES(an == bn && an == rn, BAD_SIZE);
263 int k;
264 switch(code) {
265 OPZE(0,"zipC Add",gsl_complex_add)
266 OPZE(1,"zipC Sub",gsl_complex_sub)
267 OPZE(2,"zipC Mul",gsl_complex_mul)
268 OPZE(3,"zipC Div",gsl_complex_div)
269 OPZE(4,"zipC Pow",gsl_complex_pow)
270 //OPZE(5,"zipR ATan2",atan2)
271 }
272 //KCVVIEW(a);
273 //KCVVIEW(b);
274 //CVVIEW(r);
275 //gsl_vector_memcpy(V(r),V(a));
276 //int res;
277 switch(code) {
278 default: ERROR(BAD_CODE);
279 }
280}
281
282
283int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
284 return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
285}
286
287
288
289
290int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)) {
291 REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
292 DEBUGMSG("luSolveR");
293 KDMVIEW(a);
294 KDMVIEW(b);
295 DMVIEW(r);
296 int res;
297 gsl_matrix *LU = gsl_matrix_alloc(ar,ar);
298 CHECK(!LU,MEM);
299 int s;
300 gsl_permutation * p = gsl_permutation_alloc (ar);
301 CHECK(!p,MEM);
302 CHECK(gsl_matrix_memcpy(LU,M(a)),MEM);
303 res = gsl_linalg_LU_decomp(LU, p, &s);
304 CHECK(res,res);
305 int c;
306
307 for (c=0; c<bc; c++) {
308 gsl_vector_const_view colb = gsl_matrix_const_column (M(b), c);
309 gsl_vector_view colr = gsl_matrix_column (M(r), c);
310 res = gsl_linalg_LU_solve (LU, p, V(colb), V(colr));
311 CHECK(res,res);
312 }
313 gsl_permutation_free(p);
314 gsl_matrix_free(LU);
315 OK
316}
317
318
319int luSolveC(KCMAT(a),KCMAT(b),CMAT(r)) {
320 REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
321 DEBUGMSG("luSolveC");
322 KCMVIEW(a);
323 KCMVIEW(b);
324 CMVIEW(r);
325 gsl_matrix_complex *LU = gsl_matrix_complex_alloc(ar,ar);
326 CHECK(!LU,MEM);
327 int s;
328 gsl_permutation * p = gsl_permutation_alloc (ar);
329 CHECK(!p,MEM);
330 CHECK(gsl_matrix_complex_memcpy(LU,M(a)),MEM);
331 int res;
332 res = gsl_linalg_complex_LU_decomp(LU, p, &s);
333 CHECK(res,res);
334 int c;
335 for (c=0; c<bc; c++) {
336 gsl_vector_complex_const_view colb = gsl_matrix_complex_const_column (M(b), c);
337 gsl_vector_complex_view colr = gsl_matrix_complex_column (M(r), c);
338 res = gsl_linalg_complex_LU_solve (LU, p, V(colb), V(colr));
339 CHECK(res,res);
340 }
341 gsl_permutation_free(p);
342 gsl_matrix_complex_free(LU);
343 OK
344}
345
346
347int luRaux(KRMAT(a),RVEC(b)) {
348 REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
349 DEBUGMSG("luRaux");
350 KDMVIEW(a);
351 //DVVIEW(b);
352 gsl_matrix_view LU = gsl_matrix_view_array(bp,ar,ac);
353 int s;
354 gsl_permutation * p = gsl_permutation_alloc (ar);
355 CHECK(!p,MEM);
356 CHECK(gsl_matrix_memcpy(M(LU),M(a)),MEM);
357 gsl_linalg_LU_decomp(M(LU), p, &s);
358 bp[ar*ar] = s;
359 int k;
360 for (k=0; k<ar; k++) {
361 bp[ar*ar+k+1] = gsl_permutation_get(p,k);
362 }
363 gsl_permutation_free(p);
364 OK
365}
366
367int luCaux(KCMAT(a),CVEC(b)) {
368 REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
369 DEBUGMSG("luCaux");
370 KCMVIEW(a);
371 //DVVIEW(b);
372 gsl_matrix_complex_view LU = gsl_matrix_complex_view_array((double*)bp,ar,ac);
373 int s;
374 gsl_permutation * p = gsl_permutation_alloc (ar);
375 CHECK(!p,MEM);
376 CHECK(gsl_matrix_complex_memcpy(M(LU),M(a)),MEM);
377 int res;
378 res = gsl_linalg_complex_LU_decomp(M(LU), p, &s);
379 CHECK(res,res);
380 ((double*)bp)[2*ar*ar] = s;
381 ((double*)bp)[2*ar*ar+1] = 0;
382 int k;
383 for (k=0; k<ar; k++) {
384 ((double*)bp)[2*ar*ar+2*k+2] = gsl_permutation_get(p,k);
385 ((double*)bp)[2*ar*ar+2*k+2+1] = 0;
386 }
387 gsl_permutation_free(p);
388 OK
389}
390
391int svd(KRMAT(a),RMAT(u), RVEC(s),RMAT(v)) {
392 REQUIRES(ar==ur && ac==uc && ac==sn && ac==vr && ac==vc,BAD_SIZE);
393 DEBUGMSG("svd");
394 KDMVIEW(a);
395 DMVIEW(u);
396 DVVIEW(s);
397 DMVIEW(v);
398 gsl_vector *workv = gsl_vector_alloc(ac);
399 CHECK(!workv,MEM);
400 gsl_matrix *workm = gsl_matrix_alloc(ac,ac);
401 CHECK(!workm,MEM);
402 CHECK(gsl_matrix_memcpy(M(u),M(a)),MEM);
403 // int res = gsl_linalg_SV_decomp_jacobi (&U.matrix, &V.matrix, &S.vector);
404 // doesn't work
405 //int res = gsl_linalg_SV_decomp (&U.matrix, &V.matrix, &S.vector, workv);
406 int res = gsl_linalg_SV_decomp_mod (M(u), workm, M(v), V(s), workv);
407 CHECK(res,res);
408 //gsl_matrix_transpose(M(v));
409 gsl_vector_free(workv);
410 gsl_matrix_free(workm);
411 OK
412}
413
414
415// for real symmetric matrices
416int eigensystemR(KRMAT(x),RVEC(l),RMAT(v)) {
417 REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
418 DEBUGMSG("eigensystemR (gsl_eigen_symmv)");
419 KDMVIEW(x);
420 DVVIEW(l);
421 DMVIEW(v);
422 gsl_matrix *XC = gsl_matrix_alloc(xr,xr);
423 gsl_matrix_memcpy(XC,M(x)); // needed because the argument is destroyed
424 // many thanks to Nico Mahlo for the bug report
425 gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (xc);
426 int res = gsl_eigen_symmv (XC, V(l), M(v), w);
427 CHECK(res,res);
428 gsl_eigen_symmv_free (w);
429 gsl_matrix_free(XC);
430 gsl_eigen_symmv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
431 OK
432}
433
434// for hermitian matrices
435int eigensystemC(KCMAT(x),RVEC(l),CMAT(v)) {
436 REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
437 DEBUGMSG("eigensystemC");
438 KCMVIEW(x);
439 DVVIEW(l);
440 CMVIEW(v);
441 gsl_matrix_complex *XC = gsl_matrix_complex_alloc(xr,xr);
442 gsl_matrix_complex_memcpy(XC,M(x)); // again needed because the argument is destroyed
443 gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc (xc);
444 int res = gsl_eigen_hermv (XC, V(l), M(v), w);
445 CHECK(res,res);
446 gsl_eigen_hermv_free (w);
447 gsl_matrix_complex_free(XC);
448 gsl_eigen_hermv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
449 OK
450}
451
452int QR(KRMAT(x),RMAT(q),RMAT(r)) {
453 REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
454 DEBUGMSG("QR");
455 KDMVIEW(x);
456 DMVIEW(q);
457 DMVIEW(r);
458 gsl_matrix * a = gsl_matrix_alloc(xr,xc);
459 gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
460 gsl_matrix_memcpy(a,M(x));
461 int res = gsl_linalg_QR_decomp(a,tau);
462 CHECK(res,res);
463 gsl_linalg_QR_unpack(a,tau,M(q),M(r));
464 gsl_vector_free(tau);
465 gsl_matrix_free(a);
466 OK
467}
468
469int chol(KRMAT(x),RMAT(l)) {
470 REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
471 DEBUGMSG("chol");
472 KDMVIEW(x);
473 DMVIEW(l);
474 gsl_matrix_memcpy(M(l),M(x));
475 int res = gsl_linalg_cholesky_decomp(M(l));
476 CHECK(res,res);
477 int r,c;
478 for (r=0; r<xr-1; r++) {
479 for(c=r+1; c<xc; c++) {
480 lp[r*lc+c] = 0.;
481 }
482 }
483 OK
484}
485
486
487int fft(int code, KCVEC(X), CVEC(R)) {
488 REQUIRES(Xn == Rn,BAD_SIZE);
489 DEBUGMSG("fft");
490 int s = Xn;
491 gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
492 gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
493 gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
494 gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
495 gsl_blas_dcopy(&X.vector,&R.vector);
496 if(code==0) {
497 gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
498 } else {
499 gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
500 }
501 gsl_fft_complex_wavetable_free (wavetable);
502 gsl_fft_complex_workspace_free (workspace);
503 OK
504}
505
506
507int integrate_qng(double f(double, void*), double a, double b, double prec,
508 double *result, double*error) {
509 DEBUGMSG("integrate_qng");
510 gsl_function F;
511 F.function = f;
512 F.params = NULL;
513 size_t neval;
514 int res = gsl_integration_qng (&F, a,b, 0, prec, result, error, &neval);
515 CHECK(res,res);
516 OK
517}
518
519int integrate_qags(double f(double,void*), double a, double b, double prec, int w,
520 double *result, double* error) {
521 DEBUGMSG("integrate_qags");
522 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
523 gsl_function F;
524 F.function = f;
525 F.params = NULL;
526 int res = gsl_integration_qags (&F, a,b, 0, prec, w,wk, result, error);
527 CHECK(res,res);
528 gsl_integration_workspace_free (wk);
529 OK
530}
531
532int polySolve(KRVEC(a), CVEC(z)) {
533 DEBUGMSG("polySolve");
534 REQUIRES(an>1,BAD_SIZE);
535 gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
536 int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
537 CHECK(res,res);
538 gsl_poly_complex_workspace_free (w);
539 OK;
540}
541
542int matrix_fscanf(char*filename, RMAT(a)) {
543 DEBUGMSG("gsl_matrix_fscanf");
544 //printf(filename); printf("\n");
545 DMVIEW(a);
546 FILE * f = fopen(filename,"r");
547 CHECK(!f,BAD_FILE);
548 int res = gsl_matrix_fscanf(f, M(a));
549 CHECK(res,res);
550 fclose (f);
551 OK
552}
553
554//---------------------------------------------------------------
555
556typedef double Trawfun(int, double*);
557
558double only_f_aux_min(const gsl_vector*x, void *pars) {
559 Trawfun * f = (Trawfun*) pars;
560 double* p = (double*)calloc(x->size,sizeof(double));
561 int k;
562 for(k=0;k<x->size;k++) {
563 p[k] = gsl_vector_get(x,k);
564 }
565 double res = f(x->size,p);
566 free(p);
567 return res;
568}
569
570// this version returns info about intermediate steps
571int minimize(double f(int, double*), double tolsize, int maxit,
572 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
573 REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
574 DEBUGMSG("minimizeList (nmsimplex)");
575 gsl_multimin_function my_func;
576 // extract function from pars
577 my_func.f = only_f_aux_min;
578 my_func.n = xin;
579 my_func.params = f;
580 size_t iter = 0;
581 int status;
582 double size;
583 const gsl_multimin_fminimizer_type *T;
584 gsl_multimin_fminimizer *s = NULL;
585 // Initial vertex size vector
586 KDVVIEW(sz);
587 // Starting point
588 KDVVIEW(xi);
589 // Minimizer nmsimplex, without derivatives
590 T = gsl_multimin_fminimizer_nmsimplex;
591 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
592 gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
593 do {
594 status = gsl_multimin_fminimizer_iterate (s);
595 size = gsl_multimin_fminimizer_size (s);
596
597 solp[iter*solc+0] = iter;
598 solp[iter*solc+1] = s->fval;
599 solp[iter*solc+2] = size;
600
601 int k;
602 for(k=0;k<xin;k++) {
603 solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
604 }
605 status = gsl_multimin_test_size (size, tolsize);
606 iter++;
607 } while (status == GSL_CONTINUE && iter < maxit);
608 int i,j;
609 for (i=iter; i<solr; i++) {
610 solp[i*solc+0] = iter;
611 for(j=1;j<solc;j++) {
612 solp[i*solc+j]=0.;
613 }
614 }
615 gsl_multimin_fminimizer_free(s);
616 OK
617}
618
619// working with the gradient
620
621typedef struct {double (*f)(int, double*); void (*df)(int, double*, double*);} Tfdf;
622
623double f_aux_min(const gsl_vector*x, void *pars) {
624 Tfdf * fdf = ((Tfdf*) pars);
625 double* p = (double*)calloc(x->size,sizeof(double));
626 int k;
627 for(k=0;k<x->size;k++) {
628 p[k] = gsl_vector_get(x,k);
629 }
630 double res = fdf->f(x->size,p);
631 free(p);
632 return res;
633}
634
635
636void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
637 Tfdf * fdf = ((Tfdf*) pars);
638 double* p = (double*)calloc(x->size,sizeof(double));
639 double* q = (double*)calloc(x->size,sizeof(double));
640 int k;
641 for(k=0;k<x->size;k++) {
642 p[k] = gsl_vector_get(x,k);
643 }
644
645 fdf->df(x->size,p,q);
646
647 for(k=0;k<x->size;k++) {
648 gsl_vector_set(g,k,q[k]);
649 }
650 free(p);
651 free(q);
652}
653
654void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
655 *f = f_aux_min(x,pars);
656 df_aux_min(x,pars,g);
657}
658
659// conjugate gradient
660int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),
661 double initstep, double minimpar, double tolgrad, int maxit,
662 KRVEC(xi), RMAT(sol)) {
663 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
664 DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
665 gsl_multimin_function_fdf my_func;
666 // extract function from pars
667 my_func.f = f_aux_min;
668 my_func.df = df_aux_min;
669 my_func.fdf = fdf_aux_min;
670 my_func.n = xin;
671 Tfdf stfdf;
672 stfdf.f = f;
673 stfdf.df = df;
674 my_func.params = &stfdf;
675 size_t iter = 0;
676 int status;
677 const gsl_multimin_fdfminimizer_type *T;
678 gsl_multimin_fdfminimizer *s = NULL;
679 // Starting point
680 KDVVIEW(xi);
681 // conjugate gradient fr
682 T = gsl_multimin_fdfminimizer_conjugate_fr;
683 s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
684 gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
685 do {
686 status = gsl_multimin_fdfminimizer_iterate (s);
687 solp[iter*solc+0] = iter;
688 solp[iter*solc+1] = s->f;
689 int k;
690 for(k=0;k<xin;k++) {
691 solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
692 }
693 status = gsl_multimin_test_gradient (s->gradient, tolgrad);
694 iter++;
695 } while (status == GSL_CONTINUE && iter < maxit);
696 int i,j;
697 for (i=iter; i<solr; i++) {
698 solp[i*solc+0] = iter;
699 for(j=1;j<solc;j++) {
700 solp[i*solc+j]=0.;
701 }
702 }
703 gsl_multimin_fdfminimizer_free(s);
704 OK
705}
706
707
708int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
709{
710 gsl_function F;
711 F.function = f;
712 F.params = 0;
713
714 if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
715
716 if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
717
718 if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
719
720 return 0;
721}
diff --git a/lib/GSL/gsl-aux.h b/lib/GSL/gsl-aux.h
new file mode 100644
index 0000000..89bd75e
--- /dev/null
+++ b/lib/GSL/gsl-aux.h
@@ -0,0 +1,66 @@
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
13int toScalarR(int code, KRVEC(x), RVEC(r));
14/* norm2, absdif, maximum, posmax, etc. */
15
16int mapR(int code, KRVEC(x), RVEC(r));
17int mapC(int code, KCVEC(x), CVEC(r));
18/* sin cos tan etc. */
19
20int mapValR(int code, double*, KRVEC(x), RVEC(r));
21int mapValC(int code, gsl_complex*, KCVEC(x), CVEC(r));
22
23int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r));
24int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r));
25
26
27int luSolveR(KRMAT(a),KRMAT(b),RMAT(r));
28int luSolveC(KCMAT(a),KCMAT(b),CMAT(r));
29int luRaux(KRMAT(a),RVEC(b));
30int luCaux(KCMAT(a),CVEC(b));
31
32int svd(KRMAT(x),RMAT(u), RVEC(s),RMAT(v));
33
34int eigensystemR(KRMAT(x),RVEC(l),RMAT(v));
35int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
36
37int QR(KRMAT(x),RMAT(q),RMAT(r));
38
39int chol(KRMAT(x),RMAT(l));
40
41int fft(int code, KCVEC(a), CVEC(b));
42
43int integrate_qng(double f(double, void*), double a, double b, double prec,
44 double *result, double*error);
45
46int integrate_qags(double f(double,void*), double a, double b, double prec, int w,
47 double *result, double* error);
48
49int polySolve(KRVEC(a), CVEC(z));
50
51int matrix_fscanf(char*filename, RMAT(a));
52
53int minimize(double f(int, double*), double tolsize, int maxit,
54 KRVEC(xi), KRVEC(sz), RMAT(sol));
55
56int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),
57 double initstep, double minimpar, double tolgrad, int maxit,
58 KRVEC(xi), RMAT(sol));
59
60int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr);
61
62double gsl_sf_erf(double);
63double gsl_sf_erf_Z(double);
64
65int gsl_sf_bessel_J0_e(double, double*); // hmmm...
66int gsl_sf_exp_e10_e(double, double*); // HMMMMM...