summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/base/src/Internal/Devel.hs89
-rw-r--r--packages/base/src/Internal/LAPACK.hs54
-rw-r--r--packages/base/src/Internal/Matrix.hs72
-rw-r--r--packages/base/src/Internal/Sparse.hs4
-rw-r--r--packages/base/src/Internal/Util.hs6
-rw-r--r--packages/base/src/Internal/Vector.hs12
-rw-r--r--packages/base/src/Internal/Vectorized.hs38
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Devel.hs6
-rw-r--r--packages/glpk/src/Numeric/LinearProgramming.hs9
-rw-r--r--packages/gsl/src/Numeric/GSL/Fitting.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Fourier.hs5
-rw-r--r--packages/gsl/src/Numeric/GSL/Internal.hs20
-rw-r--r--packages/gsl/src/Numeric/GSL/LinearAlgebra.hs12
-rw-r--r--packages/gsl/src/Numeric/GSL/Minimization.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Polynomials.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs1
-rw-r--r--packages/gsl/src/Numeric/GSL/Vector.hs12
-rw-r--r--packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs2
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs2
20 files changed, 200 insertions, 152 deletions
diff --git a/packages/base/src/Internal/Devel.hs b/packages/base/src/Internal/Devel.hs
index b8e04ef..4be0afd 100644
--- a/packages/base/src/Internal/Devel.hs
+++ b/packages/base/src/Internal/Devel.hs
@@ -1,4 +1,5 @@
1{-# LANGUAGE TypeOperators #-} 1{-# LANGUAGE TypeOperators #-}
2{-# LANGUAGE TypeFamilies #-}
2 3
3-- | 4-- |
4-- Module : Internal.Devel 5-- Module : Internal.Devel
@@ -16,68 +17,14 @@ import Foreign.C.Types ( CInt )
16--import Foreign.Storable.Complex () 17--import Foreign.Storable.Complex ()
17import Foreign.Ptr(Ptr) 18import Foreign.Ptr(Ptr)
18import Control.Exception as E ( SomeException, catch ) 19import Control.Exception as E ( SomeException, catch )
19 20import Internal.Vector(Vector,avec,arrvec)
21import Foreign.Storable(Storable)
20 22
21-- | postfix function application (@flip ($)@) 23-- | postfix function application (@flip ($)@)
22(//) :: x -> (x -> y) -> y 24(//) :: x -> (x -> y) -> y
23infixl 0 // 25infixl 0 //
24(//) = flip ($) 26(//) = flip ($)
25 27
26-- hmm..
27ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
28ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
29ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
30ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
31ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 f = w1 o1 $ ww5 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 . f
32ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 f = w1 o1 $ ww6 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 . f
33ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 f = w1 o1 $ ww7 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 . f
34ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 f = w1 o1 $ ww8 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 . f
35ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 f = w1 o1 $ ww9 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 . f
36
37type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
38
39type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
40type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
41type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
42type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
43type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5
44type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
45type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
46type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
47type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
48type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
49
50app1 :: f -> Adapt1 f t1
51app2 :: f -> Adapt2 f t1 r1 t2
52app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
53app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
54app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
55app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
56app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
57app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
58app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
59app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
60
61app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
62app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
63app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
64 \a1 a2 a3 -> f // a1 // a2 // a3 // check s
65app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $
66 \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
67app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $
68 \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s
69app6 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 s = ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 $
70 \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s
71app7 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 s = ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 $
72 \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s
73app8 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 s = ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 $
74 \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s
75app9 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 s = ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 $
76 \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s
77app10 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 s = ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 $
78 \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s
79
80
81 28
82-- GSL error codes are <= 1024 29-- GSL error codes are <= 1024
83-- | error codes for the auxiliary functions required by the wrappers 30-- | error codes for the auxiliary functions required by the wrappers
@@ -104,6 +51,11 @@ check msg f = do
104 when (err/=0) $ error (msg++": "++errorCode err) 51 when (err/=0) $ error (msg++": "++errorCode err)
105 return () 52 return ()
106 53
54
55-- | postfix error code check
56infixl 0 #|
57(#|) = flip check
58
107-- | Error capture and conversion to Maybe 59-- | Error capture and conversion to Maybe
108mbCatch :: IO x -> IO (Maybe x) 60mbCatch :: IO x -> IO (Maybe x)
109mbCatch act = E.catch (Just `fmap` act) f 61mbCatch act = E.catch (Just `fmap` act) f
@@ -124,4 +76,27 @@ type (:>) t r = CV t r
124type (::>) t r = OM t r 76type (::>) t r = OM t r
125type (..>) t r = CM t r 77type (..>) t r = CM t r
126 78
79class TransArray c
80 where
81 type Trans c b
82 type TransRaw c b
83 type Elem c
84 apply :: (Trans c b) -> c -> b
85 applyRaw :: (TransRaw c b) -> c -> b
86 applyArray :: (Ptr CInt -> Ptr (Elem c) -> b) -> c -> b
87 infixl 1 `apply`, `applyRaw`, `applyArray`
88
89instance Storable t => TransArray (Vector t)
90 where
91 type Trans (Vector t) b = CInt -> Ptr t -> b
92 type TransRaw (Vector t) b = CInt -> Ptr t -> b
93 type Elem (Vector t) = t
94 apply = avec
95 {-# INLINE apply #-}
96 applyRaw = avec
97 {-# INLINE applyRaw #-}
98 applyArray = arrvec
99 {-# INLINE applyArray #-}
100
101
127 102
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 8df568d..3a9abbb 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -17,7 +17,7 @@ module Internal.LAPACK where
17 17
18import Internal.Devel 18import Internal.Devel
19import Internal.Vector 19import Internal.Vector
20import Internal.Matrix 20import Internal.Matrix hiding ((#))
21import Internal.Conversion 21import Internal.Conversion
22import Internal.Element 22import Internal.Element
23import Foreign.Ptr(nullPtr) 23import Foreign.Ptr(nullPtr)
@@ -27,6 +27,16 @@ import System.IO.Unsafe(unsafePerformIO)
27 27
28----------------------------------------------------------------------------------- 28-----------------------------------------------------------------------------------
29 29
30infixl 1 #
31a # b = applyRaw a b
32{-# INLINE (#) #-}
33
34infixl 1 #!
35a #! b = apply a b
36{-# INLINE (#!) #-}
37
38-----------------------------------------------------------------------------------
39
30type TMMM t = t ..> t ..> t ..> Ok 40type TMMM t = t ..> t ..> t ..> Ok
31 41
32type F = Float 42type F = Float
@@ -49,7 +59,7 @@ multiplyAux f st a b = unsafePerformIO $ do
49 when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ 59 when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
50 show (rows a,cols a) ++ " x " ++ show (rows b, cols b) 60 show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
51 s <- createMatrix ColumnMajor (rows a) (cols b) 61 s <- createMatrix ColumnMajor (rows a) (cols b)
52 app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st 62 f (isT a) (isT b) # (tt a) # (tt b) # s #| st
53 return s 63 return s
54 64
55-- | Matrix product based on BLAS's /dgemm/. 65-- | Matrix product based on BLAS's /dgemm/.
@@ -73,7 +83,7 @@ multiplyI m a b = unsafePerformIO $ do
73 when (cols a /= rows b) $ error $ 83 when (cols a /= rows b) $ error $
74 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b 84 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
75 s <- createMatrix ColumnMajor (rows a) (cols b) 85 s <- createMatrix ColumnMajor (rows a) (cols b)
76 app3 (c_multiplyI m) omat a omat b omat s "c_multiplyI" 86 c_multiplyI m #! a #! b #! s #|"c_multiplyI"
77 return s 87 return s
78 88
79multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z 89multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
@@ -81,7 +91,7 @@ multiplyL m a b = unsafePerformIO $ do
81 when (cols a /= rows b) $ error $ 91 when (cols a /= rows b) $ error $
82 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b 92 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
83 s <- createMatrix ColumnMajor (rows a) (cols b) 93 s <- createMatrix ColumnMajor (rows a) (cols b)
84 app3 (c_multiplyL m) omat a omat b omat s "c_multiplyL" 94 c_multiplyL m #! a #! b #! s #|"c_multiplyL"
85 return s 95 return s
86 96
87----------------------------------------------------------------------------- 97-----------------------------------------------------------------------------
@@ -113,7 +123,7 @@ svdAux f st x = unsafePerformIO $ do
113 u <- createMatrix ColumnMajor r r 123 u <- createMatrix ColumnMajor r r
114 s <- createVector (min r c) 124 s <- createVector (min r c)
115 v <- createMatrix ColumnMajor c c 125 v <- createMatrix ColumnMajor c c
116 app4 f mat x mat u vec s mat v st 126 f # x # u # s # v #| st
117 return (u,s,v) 127 return (u,s,v)
118 where r = rows x 128 where r = rows x
119 c = cols x 129 c = cols x
@@ -139,7 +149,7 @@ thinSVDAux f st x = unsafePerformIO $ do
139 u <- createMatrix ColumnMajor r q 149 u <- createMatrix ColumnMajor r q
140 s <- createVector q 150 s <- createVector q
141 v <- createMatrix ColumnMajor q c 151 v <- createMatrix ColumnMajor q c
142 app4 f mat x mat u vec s mat v st 152 f # x # u # s # v #| st
143 return (u,s,v) 153 return (u,s,v)
144 where r = rows x 154 where r = rows x
145 c = cols x 155 c = cols x
@@ -164,7 +174,7 @@ svCd = svAux zgesdd "svCd" . fmat
164 174
165svAux f st x = unsafePerformIO $ do 175svAux f st x = unsafePerformIO $ do
166 s <- createVector q 176 s <- createVector q
167 app2 g mat x vec s st 177 g # x # s #| st
168 return s 178 return s
169 where r = rows x 179 where r = rows x
170 c = cols x 180 c = cols x
@@ -183,7 +193,7 @@ rightSVC = rightSVAux zgesvd "rightSVC" . fmat
183rightSVAux f st x = unsafePerformIO $ do 193rightSVAux f st x = unsafePerformIO $ do
184 s <- createVector q 194 s <- createVector q
185 v <- createMatrix ColumnMajor c c 195 v <- createMatrix ColumnMajor c c
186 app3 g mat x vec s mat v st 196 g # x # s # v #| st
187 return (s,v) 197 return (s,v)
188 where r = rows x 198 where r = rows x
189 c = cols x 199 c = cols x
@@ -202,7 +212,7 @@ leftSVC = leftSVAux zgesvd "leftSVC" . fmat
202leftSVAux f st x = unsafePerformIO $ do 212leftSVAux f st x = unsafePerformIO $ do
203 u <- createMatrix ColumnMajor r r 213 u <- createMatrix ColumnMajor r r
204 s <- createVector q 214 s <- createVector q
205 app3 g mat x mat u vec s st 215 g # x # u # s #| st
206 return (u,s) 216 return (u,s)
207 where r = rows x 217 where r = rows x
208 c = cols x 218 c = cols x
@@ -219,7 +229,7 @@ foreign import ccall unsafe "eig_l_H" zheev :: CInt -> C ..> R :> C ..> Ok
219eigAux f st m = unsafePerformIO $ do 229eigAux f st m = unsafePerformIO $ do
220 l <- createVector r 230 l <- createVector r
221 v <- createMatrix ColumnMajor r r 231 v <- createMatrix ColumnMajor r r
222 app3 g mat m vec l mat v st 232 g # m # l # v #| st
223 return (l,v) 233 return (l,v)
224 where r = rows m 234 where r = rows m
225 g ra ca pa = f ra ca pa 0 0 nullPtr 235 g ra ca pa = f ra ca pa 0 0 nullPtr
@@ -232,7 +242,7 @@ eigC = eigAux zgeev "eigC" . fmat
232 242
233eigOnlyAux f st m = unsafePerformIO $ do 243eigOnlyAux f st m = unsafePerformIO $ do
234 l <- createVector r 244 l <- createVector r
235 app2 g mat m vec l st 245 g # m # l #| st
236 return l 246 return l
237 where r = rows m 247 where r = rows m
238 g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr 248 g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr
@@ -255,7 +265,7 @@ eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
255eigRaux m = unsafePerformIO $ do 265eigRaux m = unsafePerformIO $ do
256 l <- createVector r 266 l <- createVector r
257 v <- createMatrix ColumnMajor r r 267 v <- createMatrix ColumnMajor r r
258 app3 g mat m vec l mat v "eigR" 268 g # m # l # v #| "eigR"
259 return (l,v) 269 return (l,v)
260 where r = rows m 270 where r = rows m
261 g ra ca pa = dgeev ra ca pa 0 0 nullPtr 271 g ra ca pa = dgeev ra ca pa 0 0 nullPtr
@@ -282,7 +292,7 @@ eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
282eigSHAux f st m = unsafePerformIO $ do 292eigSHAux f st m = unsafePerformIO $ do
283 l <- createVector r 293 l <- createVector r
284 v <- createMatrix ColumnMajor r r 294 v <- createMatrix ColumnMajor r r
285 app3 f mat m vec l mat v st 295 f # m # l # v #| st
286 return (l,v) 296 return (l,v)
287 where r = rows m 297 where r = rows m
288 298
@@ -332,7 +342,7 @@ foreign import ccall unsafe "cholSolveC_l" zpotrs :: TMMM C
332linearSolveSQAux g f st a b 342linearSolveSQAux g f st a b
333 | n1==n2 && n1==r = unsafePerformIO . g $ do 343 | n1==n2 && n1==r = unsafePerformIO . g $ do
334 s <- createMatrix ColumnMajor r c 344 s <- createMatrix ColumnMajor r c
335 app3 f mat a mat b mat s st 345 f # a # b # s #| st
336 return s 346 return s
337 | otherwise = error $ st ++ " of nonsquare matrix" 347 | otherwise = error $ st ++ " of nonsquare matrix"
338 where n1 = rows a 348 where n1 = rows a
@@ -371,7 +381,7 @@ foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TMMM C
371 381
372linearSolveAux f st a b = unsafePerformIO $ do 382linearSolveAux f st a b = unsafePerformIO $ do
373 r <- createMatrix ColumnMajor (max m n) nrhs 383 r <- createMatrix ColumnMajor (max m n) nrhs
374 app3 f mat a mat b mat r st 384 f # a # b # r #| st
375 return r 385 return r
376 where m = rows a 386 where m = rows a
377 n = cols a 387 n = cols a
@@ -412,7 +422,7 @@ foreign import ccall unsafe "chol_l_S" dpotrf :: TMM R
412 422
413cholAux f st a = do 423cholAux f st a = do
414 r <- createMatrix ColumnMajor n n 424 r <- createMatrix ColumnMajor n n
415 app2 f mat a mat r st 425 f # a # r #| st
416 return r 426 return r
417 where n = rows a 427 where n = rows a
418 428
@@ -450,7 +460,7 @@ qrC = qrAux zgeqr2 "qrC" . fmat
450qrAux f st a = unsafePerformIO $ do 460qrAux f st a = unsafePerformIO $ do
451 r <- createMatrix ColumnMajor m n 461 r <- createMatrix ColumnMajor m n
452 tau <- createVector mn 462 tau <- createVector mn
453 app3 f mat a vec tau mat r st 463 f # a # tau # r #| st
454 return (r,tau) 464 return (r,tau)
455 where 465 where
456 m = rows a 466 m = rows a
@@ -469,7 +479,7 @@ qrgrC = qrgrAux zungqr "qrgrC"
469 479
470qrgrAux f st n (a, tau) = unsafePerformIO $ do 480qrgrAux f st n (a, tau) = unsafePerformIO $ do
471 res <- createMatrix ColumnMajor (rows a) n 481 res <- createMatrix ColumnMajor (rows a) n
472 app3 f mat (fmat a) vec (subVector 0 n tau') mat res st 482 f # (fmat a) # (subVector 0 n tau') # res #| st
473 return res 483 return res
474 where 484 where
475 tau' = vjoin [tau, constantD 0 n] 485 tau' = vjoin [tau, constantD 0 n]
@@ -489,7 +499,7 @@ hessC = hessAux zgehrd "hessC" . fmat
489hessAux f st a = unsafePerformIO $ do 499hessAux f st a = unsafePerformIO $ do
490 r <- createMatrix ColumnMajor m n 500 r <- createMatrix ColumnMajor m n
491 tau <- createVector (mn-1) 501 tau <- createVector (mn-1)
492 app3 f mat a vec tau mat r st 502 f # a # tau # r #| st
493 return (r,tau) 503 return (r,tau)
494 where m = rows a 504 where m = rows a
495 n = cols a 505 n = cols a
@@ -510,7 +520,7 @@ schurC = schurAux zgees "schurC" . fmat
510schurAux f st a = unsafePerformIO $ do 520schurAux f st a = unsafePerformIO $ do
511 u <- createMatrix ColumnMajor n n 521 u <- createMatrix ColumnMajor n n
512 s <- createMatrix ColumnMajor n n 522 s <- createMatrix ColumnMajor n n
513 app3 f mat a mat u mat s st 523 f # a # u # s #| st
514 return (u,s) 524 return (u,s)
515 where n = rows a 525 where n = rows a
516 526
@@ -529,7 +539,7 @@ luC = luAux zgetrf "luC" . fmat
529luAux f st a = unsafePerformIO $ do 539luAux f st a = unsafePerformIO $ do
530 lu <- createMatrix ColumnMajor n m 540 lu <- createMatrix ColumnMajor n m
531 piv <- createVector (min n m) 541 piv <- createVector (min n m)
532 app3 f mat a vec piv mat lu st 542 f # a # piv # lu #| st
533 return (lu, map (pred.round) (toList piv)) 543 return (lu, map (pred.round) (toList piv))
534 where n = rows a 544 where n = rows a
535 m = cols a 545 m = cols a
@@ -552,7 +562,7 @@ lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
552lusAux f st a piv b 562lusAux f st a piv b
553 | n1==n2 && n2==n =unsafePerformIO $ do 563 | n1==n2 && n2==n =unsafePerformIO $ do
554 x <- createMatrix ColumnMajor n m 564 x <- createMatrix ColumnMajor n m
555 app4 f mat a vec piv' mat b mat x st 565 f # a # piv' # b # x #| st
556 return x 566 return x
557 | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" 567 | otherwise = error $ st ++ " on LU factorization of nonsquare matrix"
558 where n1 = rows a 568 where n1 = rows a
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index 8f8c219..db0a609 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -3,6 +3,8 @@
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE BangPatterns #-} 4{-# LANGUAGE BangPatterns #-}
5{-# LANGUAGE TypeOperators #-} 5{-# LANGUAGE TypeOperators #-}
6{-# LANGUAGE TypeFamilies #-}
7
6 8
7-- | 9-- |
8-- Module : Internal.Matrix 10-- Module : Internal.Matrix
@@ -18,7 +20,7 @@ module Internal.Matrix where
18 20
19import Internal.Vector 21import Internal.Vector
20import Internal.Devel 22import Internal.Devel
21import Internal.Vectorized 23import Internal.Vectorized hiding ((#))
22import Foreign.Marshal.Alloc ( free ) 24import Foreign.Marshal.Alloc ( free )
23import Foreign.Marshal.Array(newArray) 25import Foreign.Marshal.Array(newArray)
24import Foreign.Ptr ( Ptr ) 26import Foreign.Ptr ( Ptr )
@@ -79,8 +81,6 @@ data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int
79-- RowMajor: preferred by C, fdat may require a transposition 81-- RowMajor: preferred by C, fdat may require a transposition
80-- ColumnMajor: preferred by LAPACK, cdat may require a transposition 82-- ColumnMajor: preferred by LAPACK, cdat may require a transposition
81 83
82--cdat = xdat
83--fdat = xdat
84 84
85rows :: Matrix t -> Int 85rows :: Matrix t -> Int
86rows = irows 86rows = irows
@@ -129,6 +129,48 @@ omat a f =
129 g (fi (rows a)) (fi (cols a)) (stepRow a) (stepCol a) p 129 g (fi (rows a)) (fi (cols a)) (stepRow a) (stepCol a) p
130 f m 130 f m
131 131
132--------------------------------------------------------------------------------
133
134{-# INLINE amatr #-}
135amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
136amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c))
137 where
138 r = fromIntegral (rows x)
139 c = fromIntegral (cols x)
140
141{-# INLINE amat #-}
142amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
143amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc))
144 where
145 r = fromIntegral (rows x)
146 c = fromIntegral (cols x)
147 sr = stepRow x
148 sc = stepCol x
149
150{-# INLINE arrmat #-}
151arrmat :: Storable a => (Ptr CInt -> Ptr a -> b) -> Matrix a -> b
152arrmat f x = inlinePerformIO (unsafeWith s (\p -> unsafeWith (xdat x) (return . f p)))
153 where
154 s = fromList [fi (rows x), fi (cols x), stepRow x, stepCol x]
155
156
157instance Storable t => TransArray (Matrix t)
158 where
159 type Elem (Matrix t) = t
160 type TransRaw (Matrix t) b = CInt -> CInt -> Ptr t -> b
161 type Trans (Matrix t) b = CInt -> CInt -> CInt -> CInt -> Ptr t -> b
162 apply = amat
163 {-# INLINE apply #-}
164 applyRaw = amatr
165 {-# INLINE applyRaw #-}
166 applyArray = arrmat
167 {-# INLINE applyArray #-}
168
169infixl 1 #
170a # b = apply a b
171{-# INLINE (#) #-}
172
173--------------------------------------------------------------------------------
132 174
133{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. 175{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
134 176
@@ -139,12 +181,6 @@ fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
139flatten :: Element t => Matrix t -> Vector t 181flatten :: Element t => Matrix t -> Vector t
140flatten = xdat . cmat 182flatten = xdat . cmat
141 183
142{-
143type Mt t s = Int -> Int -> Ptr t -> s
144
145infixr 6 ::>
146type t ::> s = Mt t s
147-}
148 184
149-- | the inverse of 'Data.Packed.Matrix.fromLists' 185-- | the inverse of 'Data.Packed.Matrix.fromLists'
150toLists :: (Element t) => Matrix t -> [[t]] 186toLists :: (Element t) => Matrix t -> [[t]]
@@ -445,7 +481,7 @@ extractAux f m moder vr modec vc = do
445 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr 481 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
446 nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc 482 nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
447 r <- createMatrix RowMajor nr nc 483 r <- createMatrix RowMajor nr nc
448 app4 (f moder modec) vec vr vec vc omat m omat r "extractAux" 484 f moder modec # vr # vc # m # r #|"extract"
449 return r 485 return r
450 486
451type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt)))) 487type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))
@@ -459,7 +495,7 @@ foreign import ccall unsafe "extractL" c_extractL :: Extr Z
459 495
460--------------------------------------------------------------- 496---------------------------------------------------------------
461 497
462setRectAux f i j m r = app2 (f (fi i) (fi j)) omat m omat r "setRect" 498setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect"
463 499
464type SetRect x = I -> I -> x ::> x::> Ok 500type SetRect x = I -> I -> x ::> x::> Ok
465 501
@@ -474,7 +510,7 @@ foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z
474 510
475sortG f v = unsafePerformIO $ do 511sortG f v = unsafePerformIO $ do
476 r <- createVector (dim v) 512 r <- createVector (dim v)
477 app2 f vec v vec r "sortG" 513 f # v # r #|"sortG"
478 return r 514 return r
479 515
480sortIdxD = sortG c_sort_indexD 516sortIdxD = sortG c_sort_indexD
@@ -501,7 +537,7 @@ foreign import ccall unsafe "sort_valuesL" c_sort_valL :: Z :> Z :> Ok
501 537
502compareG f u v = unsafePerformIO $ do 538compareG f u v = unsafePerformIO $ do
503 r <- createVector (dim v) 539 r <- createVector (dim v)
504 app3 f vec u vec v vec r "compareG" 540 f # u # v # r #|"compareG"
505 return r 541 return r
506 542
507compareD = compareG c_compareD 543compareD = compareG c_compareD
@@ -518,7 +554,7 @@ foreign import ccall unsafe "compareL" c_compareL :: Z :> Z :> I :> Ok
518 554
519selectG f c u v w = unsafePerformIO $ do 555selectG f c u v w = unsafePerformIO $ do
520 r <- createVector (dim v) 556 r <- createVector (dim v)
521 app5 f vec c vec u vec v vec w vec r "selectG" 557 f # c # u # v # w # r #|"selectG"
522 return r 558 return r
523 559
524selectD = selectG c_selectD 560selectD = selectG c_selectD
@@ -541,7 +577,7 @@ foreign import ccall unsafe "chooseL" c_selectL :: Sel Z
541 577
542remapG f i j m = unsafePerformIO $ do 578remapG f i j m = unsafePerformIO $ do
543 r <- createMatrix RowMajor (rows i) (cols i) 579 r <- createMatrix RowMajor (rows i) (cols i)
544 app4 f omat i omat j omat m omat r "remapG" 580 f # i # j # m # r #|"remapG"
545 return r 581 return r
546 582
547remapD = remapG c_remapD 583remapD = remapG c_remapD
@@ -564,7 +600,7 @@ foreign import ccall unsafe "remapL" c_remapL :: Rem Z
564 600
565rowOpAux f c x i1 i2 j1 j2 m = do 601rowOpAux f c x i1 i2 j1 j2 m = do
566 px <- newArray [x] 602 px <- newArray [x]
567 app1 (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) omat m "rowOp" 603 f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp"
568 free px 604 free px
569 605
570type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok 606type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok
@@ -580,7 +616,7 @@ foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z
580 616
581-------------------------------------------------------------------------------- 617--------------------------------------------------------------------------------
582 618
583gemmg f u v m1 m2 m3 = app5 f vec u vec v omat m1 omat m2 omat m3 "gemmg" 619gemmg f u v m1 m2 m3 = f # u # v # m1 # m2 # m3 #|"gemmg"
584 620
585type Tgemm x = x :> I :> x ::> x ::> x ::> Ok 621type Tgemm x = x :> I :> x ::> x ::> x ::> Ok
586 622
@@ -608,7 +644,7 @@ saveMatrix
608saveMatrix name format m = do 644saveMatrix name format m = do
609 cname <- newCString name 645 cname <- newCString name
610 cformat <- newCString format 646 cformat <- newCString format
611 app1 (c_saveMatrix cname cformat) mat m "saveMatrix" 647 c_saveMatrix cname cformat `applyRaw` m #|"saveMatrix"
612 free cname 648 free cname
613 free cformat 649 free cformat
614 return () 650 return ()
diff --git a/packages/base/src/Internal/Sparse.hs b/packages/base/src/Internal/Sparse.hs
index b365c15..eb4ee1b 100644
--- a/packages/base/src/Internal/Sparse.hs
+++ b/packages/base/src/Internal/Sparse.hs
@@ -145,13 +145,13 @@ gmXv :: GMatrix -> Vector Double -> Vector Double
145gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do 145gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do
146 dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) 146 dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
147 r <- createVector nRows 147 r <- createVector nRows
148 app5 c_smXv vec csrVals vec csrCols vec csrRows vec v vec r "CSRXv" 148 c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv"
149 return r 149 return r
150 150
151gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do 151gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do
152 dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) 152 dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
153 r <- createVector nRows 153 r <- createVector nRows
154 app5 c_smTXv vec cscVals vec cscRows vec cscCols vec v vec r "CSCXv" 154 c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv"
155 return r 155 return r
156 156
157gmXv Diag{..} v 157gmXv Diag{..} v
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 079663d..924ca4c 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -31,7 +31,7 @@ module Internal.Util(
31 diagl, 31 diagl,
32 row, 32 row,
33 col, 33 col,
34 (&), (¦), (|||), (——), (===), (#), 34 (&), (¦), (|||), (——), (===),
35 (?), (¿), 35 (?), (¿),
36 Indexable(..), size, 36 Indexable(..), size,
37 Numeric, 37 Numeric,
@@ -185,10 +185,6 @@ infixl 2 ——
185(——) = (===) 185(——) = (===)
186 186
187 187
188(#) :: Matrix Double -> Matrix Double -> Matrix Double
189infixl 2 #
190a # b = fromBlocks [[a],[b]]
191
192-- | create a single row real matrix from a list 188-- | create a single row real matrix from a list
193-- 189--
194-- >>> row [2,3,1,8] 190-- >>> row [2,3,1,8]
diff --git a/packages/base/src/Internal/Vector.hs b/packages/base/src/Internal/Vector.hs
index 0e9161d..e5ac440 100644
--- a/packages/base/src/Internal/Vector.hs
+++ b/packages/base/src/Internal/Vector.hs
@@ -14,7 +14,7 @@ module Internal.Vector(
14 I,Z,R,C, 14 I,Z,R,C,
15 fi,ti, 15 fi,ti,
16 Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith, 16 Vector, fromList, unsafeToForeignPtr, unsafeFromForeignPtr, unsafeWith,
17 createVector, vec, 17 createVector, vec, avec, arrvec, inlinePerformIO,
18 toList, dim, (@>), at', (|>), 18 toList, dim, (@>), at', (|>),
19 vjoin, subVector, takesV, idxs, 19 vjoin, subVector, takesV, idxs,
20 buildVector, 20 buildVector,
@@ -75,6 +75,16 @@ vec x f = unsafeWith x $ \p -> do
75 f v 75 f v
76{-# INLINE vec #-} 76{-# INLINE vec #-}
77 77
78{-# INLINE avec #-}
79avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
80avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v))))
81infixl 1 `avec`
82
83{-# INLINE arrvec #-}
84arrvec :: Storable a => (Ptr CInt -> Ptr a -> b) -> Vector a -> b
85arrvec f v = inlinePerformIO (unsafeWith (idxs [1,dim v]) (\p -> unsafeWith v (return . f p)))
86
87
78 88
79-- allocates memory for a new vector 89-- allocates memory for a new vector
80createVector :: Storable a => Int -> IO (Vector a) 90createVector :: Storable a => Int -> IO (Vector a)
diff --git a/packages/base/src/Internal/Vectorized.hs b/packages/base/src/Internal/Vectorized.hs
index 5c89ac9..03bcf90 100644
--- a/packages/base/src/Internal/Vectorized.hs
+++ b/packages/base/src/Internal/Vectorized.hs
@@ -1,4 +1,5 @@
1{-# LANGUAGE TypeOperators #-} 1{-# LANGUAGE TypeOperators #-}
2{-# LANGUAGE TypeFamilies #-}
2 3
3----------------------------------------------------------------------------- 4-----------------------------------------------------------------------------
4-- | 5-- |
@@ -26,7 +27,9 @@ import Foreign.C.String
26import System.IO.Unsafe(unsafePerformIO) 27import System.IO.Unsafe(unsafePerformIO)
27import Control.Monad(when) 28import Control.Monad(when)
28 29
29 30infixl 1 #
31a # b = applyRaw a b
32{-# INLINE (#) #-}
30 33
31fromei x = fromIntegral (fromEnum x) :: CInt 34fromei x = fromIntegral (fromEnum x) :: CInt
32 35
@@ -100,7 +103,7 @@ sumL m = sumg (c_sumL m)
100 103
101sumg f x = unsafePerformIO $ do 104sumg f x = unsafePerformIO $ do
102 r <- createVector 1 105 r <- createVector 1
103 app2 f vec x vec r "sum" 106 f # x # r #| "sum"
104 return $ r @> 0 107 return $ r @> 0
105 108
106type TVV t = t :> t :> Ok 109type TVV t = t :> t :> Ok
@@ -128,14 +131,15 @@ prodQ = prodg c_prodQ
128prodC :: Vector (Complex Double) -> Complex Double 131prodC :: Vector (Complex Double) -> Complex Double
129prodC = prodg c_prodC 132prodC = prodg c_prodC
130 133
131 134prodI :: I-> Vector I -> I
132prodI = prodg . c_prodI 135prodI = prodg . c_prodI
133 136
137prodL :: Z-> Vector Z -> Z
134prodL = prodg . c_prodL 138prodL = prodg . c_prodL
135 139
136prodg f x = unsafePerformIO $ do 140prodg f x = unsafePerformIO $ do
137 r <- createVector 1 141 r <- createVector 1
138 app2 f vec x vec r "prod" 142 f # x # r #| "prod"
139 return $ r @> 0 143 return $ r @> 0
140 144
141 145
@@ -150,24 +154,24 @@ foreign import ccall unsafe "prodL" c_prodL :: Z -> TVV Z
150 154
151toScalarAux fun code v = unsafePerformIO $ do 155toScalarAux fun code v = unsafePerformIO $ do
152 r <- createVector 1 156 r <- createVector 1
153 app2 (fun (fromei code)) vec v vec r "toScalarAux" 157 fun (fromei code) # v # r #|"toScalarAux"
154 return (r @> 0) 158 return (r @> 0)
155 159
156vectorMapAux fun code v = unsafePerformIO $ do 160vectorMapAux fun code v = unsafePerformIO $ do
157 r <- createVector (dim v) 161 r <- createVector (dim v)
158 app2 (fun (fromei code)) vec v vec r "vectorMapAux" 162 fun (fromei code) # v # r #|"vectorMapAux"
159 return r 163 return r
160 164
161vectorMapValAux fun code val v = unsafePerformIO $ do 165vectorMapValAux fun code val v = unsafePerformIO $ do
162 r <- createVector (dim v) 166 r <- createVector (dim v)
163 pval <- newArray [val] 167 pval <- newArray [val]
164 app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux" 168 fun (fromei code) pval # v # r #|"vectorMapValAux"
165 free pval 169 free pval
166 return r 170 return r
167 171
168vectorZipAux fun code u v = unsafePerformIO $ do 172vectorZipAux fun code u v = unsafePerformIO $ do
169 r <- createVector (dim u) 173 r <- createVector (dim u)
170 app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux" 174 fun (fromei code) # u # v # r #|"vectorZipAux"
171 return r 175 return r
172 176
173--------------------------------------------------------------------- 177---------------------------------------------------------------------
@@ -364,7 +368,7 @@ randomVector :: Seed
364 -> Vector Double 368 -> Vector Double
365randomVector seed dist n = unsafePerformIO $ do 369randomVector seed dist n = unsafePerformIO $ do
366 r <- createVector n 370 r <- createVector n
367 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" 371 c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
368 return r 372 return r
369 373
370foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok 374foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok
@@ -373,7 +377,7 @@ foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> D
373 377
374roundVector v = unsafePerformIO $ do 378roundVector v = unsafePerformIO $ do
375 r <- createVector (dim v) 379 r <- createVector (dim v)
376 app2 c_round_vector vec v vec r "roundVector" 380 c_round_vector # v # r #|"roundVector"
377 return r 381 return r
378 382
379foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double 383foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
@@ -387,7 +391,7 @@ foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
387range :: Int -> Vector I 391range :: Int -> Vector I
388range n = unsafePerformIO $ do 392range n = unsafePerformIO $ do
389 r <- createVector n 393 r <- createVector n
390 app1 c_range_vector vec r "range" 394 c_range_vector # r #|"range"
391 return r 395 return r
392 396
393foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok 397foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok
@@ -427,7 +431,7 @@ long2intV = tog c_long2int
427 431
428tog f v = unsafePerformIO $ do 432tog f v = unsafePerformIO $ do
429 r <- createVector (dim v) 433 r <- createVector (dim v)
430 app2 f vec v vec r "tog" 434 f # v # r #|"tog"
431 return r 435 return r
432 436
433foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok 437foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok
@@ -446,7 +450,7 @@ foreign import ccall unsafe "long2int" c_long2int :: Z :> I :> Ok
446 450
447stepg f v = unsafePerformIO $ do 451stepg f v = unsafePerformIO $ do
448 r <- createVector (dim v) 452 r <- createVector (dim v)
449 app2 f vec v vec r "step" 453 f # v # r #|"step"
450 return r 454 return r
451 455
452stepD :: Vector Double -> Vector Double 456stepD :: Vector Double -> Vector Double
@@ -471,7 +475,7 @@ foreign import ccall unsafe "stepL" c_stepL :: TVV Z
471 475
472conjugateAux fun x = unsafePerformIO $ do 476conjugateAux fun x = unsafePerformIO $ do
473 v <- createVector (dim x) 477 v <- createVector (dim x)
474 app2 fun vec x vec v "conjugateAux" 478 fun # x # v #|"conjugateAux"
475 return v 479 return v
476 480
477conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) 481conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
@@ -489,7 +493,7 @@ cloneVector v = do
489 let n = dim v 493 let n = dim v
490 r <- createVector n 494 r <- createVector n
491 let f _ s _ d = copyArray d s n >> return 0 495 let f _ s _ d = copyArray d s n >> return 0
492 app2 f vec v vec r "cloneVector" 496 f # v # r #|"cloneVector"
493 return r 497 return r
494 498
495-------------------------------------------------------------------------------- 499--------------------------------------------------------------------------------
@@ -497,7 +501,7 @@ cloneVector v = do
497constantAux fun x n = unsafePerformIO $ do 501constantAux fun x n = unsafePerformIO $ do
498 v <- createVector n 502 v <- createVector n
499 px <- newArray [x] 503 px <- newArray [x]
500 app1 (fun px) vec v "constantAux" 504 fun px # v #|"constantAux"
501 free px 505 free px
502 return v 506 return v
503 507
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 2a45771..d2843c2 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -104,7 +104,7 @@ import Internal.Element
104import Internal.IO 104import Internal.IO
105import Internal.Numeric 105import Internal.Numeric
106import Internal.Container 106import Internal.Container
107import Internal.Util hiding ((&),(#)) 107import Internal.Util hiding ((&))
108import Data.Complex 108import Data.Complex
109import Internal.Sparse 109import Internal.Sparse
110import Internal.Modular 110import Internal.Modular
diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
index db4236b..f18d35b 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
@@ -23,14 +23,12 @@ module Numeric.LinearAlgebra.Devel(
23 -- | See @examples/devel@ in the repository. 23 -- | See @examples/devel@ in the repository.
24 24
25 createVector, createMatrix, 25 createVector, createMatrix,
26 vec, mat, omat, 26 TransArray(..),
27 app1, app2, app3, app4,
28 app5, app6, app7, app8, app9, app10,
29 MatrixOrder(..), orderOf, cmat, fmat, 27 MatrixOrder(..), orderOf, cmat, fmat,
30 matrixFromVector, 28 matrixFromVector,
31 unsafeFromForeignPtr, 29 unsafeFromForeignPtr,
32 unsafeToForeignPtr, 30 unsafeToForeignPtr,
33 check, (//), 31 check, (//), (#|),
34 at', atM', fi, ti, 32 at', atM', fi, ti,
35 33
36 -- * ST 34 -- * ST
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index 7bf4279..0a776fa 100644
--- a/packages/glpk/src/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -275,6 +275,11 @@ mkConstrS n objfun b1 = fromLists (ob ++ co) where
275 275
276----------------------------------------------------- 276-----------------------------------------------------
277 277
278(##) :: TransArray c => TransRaw c b -> c -> b
279infixl 1 ##
280a ## b = applyRaw a b
281{-# INLINE (##) #-}
282
278foreign import ccall unsafe "c_simplex_sparse" c_simplex_sparse 283foreign import ccall unsafe "c_simplex_sparse" c_simplex_sparse
279 :: CInt -> CInt -- rows and cols 284 :: CInt -> CInt -- rows and cols
280 -> CInt -> CInt -> Ptr Double -- coeffs 285 -> CInt -> CInt -> Ptr Double -- coeffs
@@ -285,7 +290,7 @@ foreign import ccall unsafe "c_simplex_sparse" c_simplex_sparse
285simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double 290simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
286simplexSparse m n c b = unsafePerformIO $ do 291simplexSparse m n c b = unsafePerformIO $ do
287 s <- createVector (2+n) 292 s <- createVector (2+n)
288 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse" 293 c_simplex_sparse (fi m) (fi n) ## (cmat c) ## (cmat b) ## s #|"c_simplex_sparse"
289 return s 294 return s
290 295
291foreign import ccall unsafe "c_exact_sparse" c_exact_sparse 296foreign import ccall unsafe "c_exact_sparse" c_exact_sparse
@@ -298,7 +303,7 @@ foreign import ccall unsafe "c_exact_sparse" c_exact_sparse
298exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double 303exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
299exactSparse m n c b = unsafePerformIO $ do 304exactSparse m n c b = unsafePerformIO $ do
300 s <- createVector (2+n) 305 s <- createVector (2+n)
301 app3 (c_exact_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_exact_sparse" 306 c_exact_sparse (fi m) (fi n) ## (cmat c) ## (cmat b) ## s #|"c_exact_sparse"
302 return s 307 return s
303 308
304glpFR, glpLO, glpUP, glpDB, glpFX :: Double 309glpFR, glpLO, glpUP, glpDB, glpFX :: Double
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
index db9d82f..8eb93a7 100644
--- a/packages/gsl/src/Numeric/GSL/Fitting.hs
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -87,7 +87,7 @@ nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
87 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f)) 87 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
88 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac)) 88 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
89 rawpath <- createMatrix RowMajor maxit (2+p) 89 rawpath <- createMatrix RowMajor maxit (2+p)
90 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit" 90 c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n) # xiv # rawpath #|"c_nlfit"
91 let it = round (rawpath `atIndex` (maxit-1,0)) 91 let it = round (rawpath `atIndex` (maxit-1,0))
92 path = takeRows it rawpath 92 path = takeRows it rawpath
93 [sol] = toRows $ dropRows (it-1) path 93 [sol] = toRows $ dropRows (it-1) path
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
index d824b4f..1c2c053 100644
--- a/packages/gsl/src/Numeric/GSL/Fourier.hs
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE TypeFamilies #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Fourier 4Module : Numeric.GSL.Fourier
3Copyright : (c) Alberto Ruiz 2006 5Copyright : (c) Alberto Ruiz 2006
@@ -23,7 +25,7 @@ import System.IO.Unsafe (unsafePerformIO)
23 25
24genfft code v = unsafePerformIO $ do 26genfft code v = unsafePerformIO $ do
25 r <- createVector (size v) 27 r <- createVector (size v)
26 app2 (c_fft code) vec v vec r "fft" 28 c_fft code # v # r #|"fft"
27 return r 29 return r
28 30
29foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res) 31foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res)
@@ -41,3 +43,4 @@ fft = genfft 0
41-- | The inverse of 'fft', using /gsl_fft_complex_inverse/. 43-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
42ifft :: Vector (Complex Double) -> Vector (Complex Double) 44ifft :: Vector (Complex Double) -> Vector (Complex Double)
43ifft = genfft 1 45ifft = genfft 1
46
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
index a269224..dcd3bc4 100644
--- a/packages/gsl/src/Numeric/GSL/Internal.hs
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -23,7 +23,7 @@ module Numeric.GSL.Internal(
23 createV, 23 createV,
24 createMIO, 24 createMIO,
25 module Numeric.LinearAlgebra.Devel, 25 module Numeric.LinearAlgebra.Devel,
26 check, 26 check,(#),vec, ww2,
27 Res,TV,TM,TCV,TCM 27 Res,TV,TM,TCV,TCM
28) where 28) where
29 29
@@ -35,7 +35,7 @@ import Foreign.Ptr(Ptr, FunPtr)
35import Foreign.C.Types 35import Foreign.C.Types
36import Foreign.C.String(peekCString) 36import Foreign.C.String(peekCString)
37import System.IO.Unsafe(unsafePerformIO) 37import System.IO.Unsafe(unsafePerformIO)
38import Data.Vector.Storable(unsafeWith) 38import Data.Vector.Storable as V (unsafeWith,length)
39import Control.Monad(when) 39import Control.Monad(when)
40 40
41iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) 41iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
@@ -86,12 +86,12 @@ aux_vTom f n p rr cr r = g where
86 86
87createV n fun msg = unsafePerformIO $ do 87createV n fun msg = unsafePerformIO $ do
88 r <- createVector n 88 r <- createVector n
89 app1 fun vec r msg 89 fun # r #| msg
90 return r 90 return r
91 91
92createMIO r c fun msg = do 92createMIO r c fun msg = do
93 res <- createMatrix RowMajor r c 93 res <- createMatrix RowMajor r c
94 app1 fun mat res msg 94 fun # res #| msg
95 return res 95 return res
96 96
97-------------------------------------------------------------------------------- 97--------------------------------------------------------------------------------
@@ -123,3 +123,15 @@ type TCM x = CInt -> CInt -> PC -> x
123type TVV = TV (TV Res) 123type TVV = TV (TV Res)
124type TVM = TV (TM Res) 124type TVM = TV (TM Res)
125 125
126ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
127
128vec x f = unsafeWith x $ \p -> do
129 let v g = do
130 g (fi $ V.length x) p
131 f v
132{-# INLINE vec #-}
133
134infixl 1 #
135a # b = applyRaw a b
136{-# INLINE (#) #-}
137
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
index cb78bf4..6ffe306 100644
--- a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -40,7 +40,7 @@ randomVector :: Int -- ^ seed
40 -> Vector Double 40 -> Vector Double
41randomVector seed dist n = unsafePerformIO $ do 41randomVector seed dist n = unsafePerformIO $ do
42 r <- createVector n 42 r <- createVector n
43 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" 43 c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
44 return r 44 return r
45 45
46foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV 46foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
@@ -56,7 +56,7 @@ saveMatrix filename fmt m = do
56 charname <- newCString filename 56 charname <- newCString filename
57 charfmt <- newCString fmt 57 charfmt <- newCString fmt
58 let o = if orderOf m == RowMajor then 1 else 0 58 let o = if orderOf m == RowMajor then 1 else 0
59 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" 59 matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
60 free charname 60 free charname
61 free charfmt 61 free charfmt
62 62
@@ -69,7 +69,7 @@ fscanfVector :: FilePath -> Int -> IO (Vector Double)
69fscanfVector filename n = do 69fscanfVector filename n = do
70 charname <- newCString filename 70 charname <- newCString filename
71 res <- createVector n 71 res <- createVector n
72 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" 72 gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
73 free charname 73 free charname
74 return res 74 return res
75 75
@@ -80,7 +80,7 @@ fprintfVector :: FilePath -> String -> Vector Double -> IO ()
80fprintfVector filename fmt v = do 80fprintfVector filename fmt v = do
81 charname <- newCString filename 81 charname <- newCString filename
82 charfmt <- newCString fmt 82 charfmt <- newCString fmt
83 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" 83 gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
84 free charname 84 free charname
85 free charfmt 85 free charfmt
86 86
@@ -91,7 +91,7 @@ freadVector :: FilePath -> Int -> IO (Vector Double)
91freadVector filename n = do 91freadVector filename n = do
92 charname <- newCString filename 92 charname <- newCString filename
93 res <- createVector n 93 res <- createVector n
94 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" 94 gsl_vector_fread charname # res #| "gsl_vector_fread"
95 free charname 95 free charname
96 return res 96 return res
97 97
@@ -101,7 +101,7 @@ foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
101fwriteVector :: FilePath -> Vector Double -> IO () 101fwriteVector :: FilePath -> Vector Double -> IO ()
102fwriteVector filename v = do 102fwriteVector filename v = do
103 charname <- newCString filename 103 charname <- newCString filename
104 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" 104 gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
105 free charname 105 free charname
106 106
107foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV 107foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
index 00e0619..a0e5306 100644
--- a/packages/gsl/src/Numeric/GSL/Minimization.hs
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -137,7 +137,7 @@ minimizeV :: MinimizeMethod
137minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi) 137minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
138 where v2l (v,m) = (toList v, m) 138 where v2l (v,m) = (toList v, m)
139 139
140ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 140
141 141
142minimizeV method eps maxit szv f xiv = unsafePerformIO $ do 142minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
143 let n = size xiv 143 let n = size xiv
diff --git a/packages/gsl/src/Numeric/GSL/Polynomials.hs b/packages/gsl/src/Numeric/GSL/Polynomials.hs
index 246e301..8890f8f 100644
--- a/packages/gsl/src/Numeric/GSL/Polynomials.hs
+++ b/packages/gsl/src/Numeric/GSL/Polynomials.hs
@@ -48,7 +48,7 @@ polySolve = toList . polySolve' . fromList
48polySolve' :: Vector Double -> Vector (Complex Double) 48polySolve' :: Vector Double -> Vector (Complex Double)
49polySolve' v | size v > 1 = unsafePerformIO $ do 49polySolve' v | size v > 1 = unsafePerformIO $ do
50 r <- createVector (size v-1) 50 r <- createVector (size v-1)
51 app2 c_polySolve vec v vec r "polySolve" 51 c_polySolve # v # r #| "polySolve"
52 return r 52 return r
53 | otherwise = error "polySolve on a polynomial of degree zero" 53 | otherwise = error "polySolve on a polynomial of degree zero"
54 54
diff --git a/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs
index 9f9ed97..11b22d3 100644
--- a/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs
+++ b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs
@@ -55,7 +55,6 @@ import Foreign.Ptr(Ptr, FunPtr, nullFunPtr)
55import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr) 55import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr)
56import Foreign.C.Types 56import Foreign.C.Types
57import System.IO.Unsafe(unsafePerformIO) 57import System.IO.Unsafe(unsafePerformIO)
58import Control.Applicative ((<*>), (<$>))
59 58
60import System.IO (hFlush, stdout) 59import System.IO (hFlush, stdout)
61 60
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
index 0cd99eb..fb982c5 100644
--- a/packages/gsl/src/Numeric/GSL/Vector.hs
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -34,7 +34,7 @@ randomVector :: Int -- ^ seed
34 -> Vector Double 34 -> Vector Double
35randomVector seed dist n = unsafePerformIO $ do 35randomVector seed dist n = unsafePerformIO $ do
36 r <- createVector n 36 r <- createVector n
37 app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL" 37 c_random_vector_GSL (fi seed) ((fi.fromEnum) dist) # r #|"randomVectorGSL"
38 return r 38 return r
39 39
40foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV 40foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
@@ -50,7 +50,7 @@ saveMatrix filename fmt m = do
50 charname <- newCString filename 50 charname <- newCString filename
51 charfmt <- newCString fmt 51 charfmt <- newCString fmt
52 let o = if orderOf m == RowMajor then 1 else 0 52 let o = if orderOf m == RowMajor then 1 else 0
53 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" 53 matrix_fprintf charname charfmt o # m #|"matrix_fprintf"
54 free charname 54 free charname
55 free charfmt 55 free charfmt
56 56
@@ -63,7 +63,7 @@ fscanfVector :: FilePath -> Int -> IO (Vector Double)
63fscanfVector filename n = do 63fscanfVector filename n = do
64 charname <- newCString filename 64 charname <- newCString filename
65 res <- createVector n 65 res <- createVector n
66 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" 66 gsl_vector_fscanf charname # res #|"gsl_vector_fscanf"
67 free charname 67 free charname
68 return res 68 return res
69 69
@@ -74,7 +74,7 @@ fprintfVector :: FilePath -> String -> Vector Double -> IO ()
74fprintfVector filename fmt v = do 74fprintfVector filename fmt v = do
75 charname <- newCString filename 75 charname <- newCString filename
76 charfmt <- newCString fmt 76 charfmt <- newCString fmt
77 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" 77 gsl_vector_fprintf charname charfmt # v #|"gsl_vector_fprintf"
78 free charname 78 free charname
79 free charfmt 79 free charfmt
80 80
@@ -85,7 +85,7 @@ freadVector :: FilePath -> Int -> IO (Vector Double)
85freadVector filename n = do 85freadVector filename n = do
86 charname <- newCString filename 86 charname <- newCString filename
87 res <- createVector n 87 res <- createVector n
88 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" 88 gsl_vector_fread charname # res #|"gsl_vector_fread"
89 free charname 89 free charname
90 return res 90 return res
91 91
@@ -95,7 +95,7 @@ foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
95fwriteVector :: FilePath -> Vector Double -> IO () 95fwriteVector :: FilePath -> Vector Double -> IO ()
96fwriteVector filename v = do 96fwriteVector filename v = do
97 charname <- newCString filename 97 charname <- newCString filename
98 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" 98 gsl_vector_fwrite charname # v #|"gsl_vector_fwrite"
99 free charname 99 free charname
100 100
101foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV 101foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
diff --git a/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
index d75fec2..b2ca7f0 100644
--- a/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
+++ b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
@@ -27,7 +27,7 @@ dss :: CSR -> Vector Double -> Vector Double
27dss CSR{..} b = unsafePerformIO $ do 27dss CSR{..} b = unsafePerformIO $ do
28 size b /= csrNRows ??? printf "dss: incorrect sizes: (%d,%d) x %d" csrNRows csrNCols (size b) 28 size b /= csrNRows ??? printf "dss: incorrect sizes: (%d,%d) x %d" csrNRows csrNCols (size b)
29 r <- createVector csrNCols 29 r <- createVector csrNCols
30 app5 c_dss vec csrVals vec csrCols vec csrRows vec b vec r "dss" 30 c_dss `apply` csrVals `apply` csrCols `apply` csrRows `apply` b `apply` r #|"dss"
31 return r 31 return r
32 32
33foreign import ccall unsafe "dss" 33foreign import ccall unsafe "dss"
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index b1428fb..d1fc6ec 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -29,7 +29,7 @@ module Numeric.LinearAlgebra.Tests(
29) where 29) where
30 30
31import Numeric.LinearAlgebra hiding (unitary) 31import Numeric.LinearAlgebra hiding (unitary)
32import Numeric.LinearAlgebra.Devel hiding (vec) 32import Numeric.LinearAlgebra.Devel
33import Numeric.LinearAlgebra.Static(L) 33import Numeric.LinearAlgebra.Static(L)
34import Numeric.LinearAlgebra.Tests.Instances 34import Numeric.LinearAlgebra.Tests.Instances
35import Numeric.LinearAlgebra.Tests.Properties 35import Numeric.LinearAlgebra.Tests.Properties