summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/Numeric/GSL/Vector.hs34
-rw-r--r--lib/Numeric/GSL/gsl-aux.c50
-rw-r--r--lib/Numeric/LinearAlgebra/Linear.hs105
3 files changed, 139 insertions, 50 deletions
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 14ba0ff..e8df27c 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -15,6 +15,7 @@
15 15
16module Numeric.GSL.Vector ( 16module Numeric.GSL.Vector (
17 sumF, sumR, sumQ, sumC, 17 sumF, sumR, sumQ, sumC,
18 prodF, prodR, prodQ, prodC,
18 dotF, dotR, dotQ, dotC, 19 dotF, dotR, dotQ, dotC,
19 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ, 20 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
20 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, 21 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
@@ -111,6 +112,39 @@ foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV
111foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV 112foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV
112foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV 113foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV
113 114
115-- | product of elements
116prodF :: Vector Float -> Float
117prodF x = unsafePerformIO $ do
118 r <- createVector 1
119 app2 c_prodF vec x vec r "prodF"
120 return $ r @> 0
121
122-- | product of elements
123prodR :: Vector Double -> Double
124prodR x = unsafePerformIO $ do
125 r <- createVector 1
126 app2 c_prodR vec x vec r "prodR"
127 return $ r @> 0
128
129-- | product of elements
130prodQ :: Vector (Complex Float) -> Complex Float
131prodQ x = unsafePerformIO $ do
132 r <- createVector 1
133 app2 c_prodQ vec x vec r "prodQ"
134 return $ r @> 0
135
136-- | product of elements
137prodC :: Vector (Complex Double) -> Complex Double
138prodC x = unsafePerformIO $ do
139 r <- createVector 1
140 app2 c_prodC vec x vec r "prodC"
141 return $ r @> 0
142
143foreign import ccall safe "gsl-aux.h prodF" c_prodF :: TFF
144foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV
145foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV
146foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV
147
114-- | dot product 148-- | dot product
115dotF :: Vector Float -> Vector Float -> Float 149dotF :: Vector Float -> Vector Float -> Float
116dotF x y = unsafePerformIO $ do 150dotF x y = unsafePerformIO $ do
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 689989d..3e9e2cb 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -153,6 +153,56 @@ int sumC(KCVEC(x),CVEC(r)) {
153 OK 153 OK
154} 154}
155 155
156int prodF(KFVEC(x),FVEC(r)) {
157 DEBUGMSG("prodF");
158 REQUIRES(rn==1,BAD_SIZE);
159 int i;
160 float res = 0;
161 for (i = 0; i < xn; i++) res *= xp[i];
162 rp[0] = res;
163 OK
164}
165
166int prodR(KRVEC(x),RVEC(r)) {
167 DEBUGMSG("prodR");
168 REQUIRES(rn==1,BAD_SIZE);
169 int i;
170 double res = 0;
171 for (i = 0; i < xn; i++) res *= xp[i];
172 rp[0] = res;
173 OK
174}
175
176int prodQ(KQVEC(x),QVEC(r)) {
177 DEBUGMSG("prodQ");
178 REQUIRES(rn==1,BAD_SIZE);
179 int i;
180 gsl_complex_float res;
181 res.dat[0] = 0;
182 res.dat[1] = 0;
183 for (i = 0; i < xn; i++) {
184 res.dat[0] = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
185 res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
186 }
187 rp[0] = res;
188 OK
189}
190
191int prodC(KCVEC(x),CVEC(r)) {
192 DEBUGMSG("prodC");
193 REQUIRES(rn==1,BAD_SIZE);
194 int i;
195 gsl_complex res;
196 res.dat[0] = 0;
197 res.dat[1] = 0;
198 for (i = 0; i < xn; i++) {
199 res.dat[0] = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
200 res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
201 }
202 rp[0] = res;
203 OK
204}
205
156int dotF(KFVEC(x), KFVEC(y), FVEC(r)) { 206int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
157 DEBUGMSG("dotF"); 207 DEBUGMSG("dotF");
158 REQUIRES(xn==yn,BAD_SIZE); 208 REQUIRES(xn==yn,BAD_SIZE);
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 4d7f608..778b976 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -39,65 +39,70 @@ import Control.Monad(ap)
39-- | basic Vector functions 39-- | basic Vector functions
40class Num e => Vectors a e where 40class Num e => Vectors a e where
41 -- the C functions sumX are twice as fast as using foldVector 41 -- the C functions sumX are twice as fast as using foldVector
42 vectorSum :: a e -> e 42 vectorSum :: a e -> e
43 absSum :: a e -> e 43 vectorProd :: a e -> e
44 vectorMin :: a e -> e 44 absSum :: a e -> e
45 vectorMax :: a e -> e 45 vectorMin :: a e -> e
46 minIdx :: a e -> Int 46 vectorMax :: a e -> e
47 maxIdx :: a e -> Int 47 minIdx :: a e -> Int
48 dot :: a e -> a e -> e 48 maxIdx :: a e -> Int
49 norm1 :: a e -> e 49 dot :: a e -> a e -> e
50 norm2 :: a e -> e 50 norm1 :: a e -> e
51 normInf :: a e -> e 51 norm2 :: a e -> e
52 normInf :: a e -> e
52 53
53 54
54instance Vectors Vector Float where 55instance Vectors Vector Float where
55 vectorSum = sumF 56 vectorSum = sumF
56 norm2 = toScalarF Norm2 57 vectorProd = prodF
57 absSum = toScalarF AbsSum 58 norm2 = toScalarF Norm2
58 vectorMin = toScalarF Min 59 absSum = toScalarF AbsSum
59 vectorMax = toScalarF Max 60 vectorMin = toScalarF Min
60 minIdx = round . toScalarF MinIdx 61 vectorMax = toScalarF Max
61 maxIdx = round . toScalarF MaxIdx 62 minIdx = round . toScalarF MinIdx
62 dot = dotF 63 maxIdx = round . toScalarF MaxIdx
63 norm1 = toScalarF AbsSum 64 dot = dotF
64 normInf = vectorMax . vectorMapF Abs 65 norm1 = toScalarF AbsSum
66 normInf = vectorMax . vectorMapF Abs
65 67
66instance Vectors Vector Double where 68instance Vectors Vector Double where
67 vectorSum = sumR 69 vectorSum = sumR
68 norm2 = toScalarR Norm2 70 vectorProd = prodR
69 absSum = toScalarR AbsSum 71 norm2 = toScalarR Norm2
70 vectorMin = toScalarR Min 72 absSum = toScalarR AbsSum
71 vectorMax = toScalarR Max 73 vectorMin = toScalarR Min
72 minIdx = round . toScalarR MinIdx 74 vectorMax = toScalarR Max
73 maxIdx = round . toScalarR MaxIdx 75 minIdx = round . toScalarR MinIdx
74 dot = dotR 76 maxIdx = round . toScalarR MaxIdx
75 norm1 = toScalarR AbsSum 77 dot = dotR
76 normInf = vectorMax . vectorMapR Abs 78 norm1 = toScalarR AbsSum
79 normInf = vectorMax . vectorMapR Abs
77 80
78instance Vectors Vector (Complex Float) where 81instance Vectors Vector (Complex Float) where
79 vectorSum = sumQ 82 vectorSum = sumQ
80 norm2 = (:+ 0) . toScalarQ Norm2 83 vectorProd = prodQ
81 absSum = (:+ 0) . toScalarQ AbsSum 84 norm2 = (:+ 0) . toScalarQ Norm2
82 vectorMin = ap (@>) minIdx 85 absSum = (:+ 0) . toScalarQ AbsSum
83 vectorMax = ap (@>) maxIdx 86 vectorMin = ap (@>) minIdx
84 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 87 vectorMax = ap (@>) maxIdx
85 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 88 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
86 dot = dotQ 89 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
87 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs 90 dot = dotQ
88 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs 91 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
92 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
89 93
90instance Vectors Vector (Complex Double) where 94instance Vectors Vector (Complex Double) where
91 vectorSum = sumC 95 vectorSum = sumC
92 norm2 = (:+ 0) . toScalarC Norm2 96 vectorProd = prodC
93 absSum = (:+ 0) . toScalarC AbsSum 97 norm2 = (:+ 0) . toScalarC Norm2
94 vectorMin = ap (@>) minIdx 98 absSum = (:+ 0) . toScalarC AbsSum
95 vectorMax = ap (@>) maxIdx 99 vectorMin = ap (@>) minIdx
96 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 100 vectorMax = ap (@>) maxIdx
97 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 101 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
98 dot = dotC 102 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
99 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs 103 dot = dotC
100 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs 104 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
105 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
101 106
102---------------------------------------------------- 107----------------------------------------------------
103 108