summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Data/Packed/Internal
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Data/Packed/Internal')
-rw-r--r--packages/hmatrix/src/Data/Packed/Internal/Common.hs171
-rw-r--r--packages/hmatrix/src/Data/Packed/Internal/Matrix.hs491
-rw-r--r--packages/hmatrix/src/Data/Packed/Internal/Signatures.hs72
-rw-r--r--packages/hmatrix/src/Data/Packed/Internal/Vector.hs521
4 files changed, 1255 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Common.hs b/packages/hmatrix/src/Data/Packed/Internal/Common.hs
new file mode 100644
index 0000000..edef3c2
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Common.hs
@@ -0,0 +1,171 @@
1{-# LANGUAGE CPP #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Internal.Common
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-- Development utilities.
13--
14-----------------------------------------------------------------------------
15-- #hide
16
17module Data.Packed.Internal.Common(
18 Adapt,
19 app1, app2, app3, app4,
20 app5, app6, app7, app8, app9, app10,
21 (//), check, mbCatch,
22 splitEvery, common, compatdim,
23 fi,
24 table
25) where
26
27import Foreign
28import Control.Monad(when)
29import Foreign.C.String(peekCString)
30import Foreign.C.Types
31import Foreign.Storable.Complex()
32import Data.List(transpose,intersperse)
33import Control.Exception as E
34
35-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
36splitEvery :: Int -> [a] -> [[a]]
37splitEvery _ [] = []
38splitEvery k l = take k l : splitEvery k (drop k l)
39
40-- | obtains the common value of a property of a list
41common :: (Eq a) => (b->a) -> [b] -> Maybe a
42common f = commonval . map f where
43 commonval :: (Eq a) => [a] -> Maybe a
44 commonval [] = Nothing
45 commonval [a] = Just a
46 commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
47
48-- | common value with \"adaptable\" 1
49compatdim :: [Int] -> Maybe Int
50compatdim [] = Nothing
51compatdim [a] = Just a
52compatdim (a:b:xs)
53 | a==b = compatdim (b:xs)
54 | a==1 = compatdim (b:xs)
55 | b==1 = compatdim (a:xs)
56 | otherwise = Nothing
57
58-- | Formatting tool
59table :: String -> [[String]] -> String
60table sep as = unlines . map unwords' $ transpose mtp where
61 mt = transpose as
62 longs = map (maximum . map length) mt
63 mtp = zipWith (\a b -> map (pad a) b) longs mt
64 pad n str = replicate (n - length str) ' ' ++ str
65 unwords' = concat . intersperse sep
66
67-- | postfix function application (@flip ($)@)
68(//) :: x -> (x -> y) -> y
69infixl 0 //
70(//) = flip ($)
71
72-- | specialized fromIntegral
73fi :: Int -> CInt
74fi = fromIntegral
75
76-- hmm..
77ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
78ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
79ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
80ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
81ww6 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
82ww7 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
83ww8 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
84ww9 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
85ww10 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
86
87type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
88
89type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
90type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
91type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
92type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
93type 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
94type 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
95type 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
96type 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
97type 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
98type 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
99
100app1 :: f -> Adapt1 f t1
101app2 :: f -> Adapt2 f t1 r1 t2
102app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
103app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
104app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
105app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
106app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
107app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
108app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
109app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
110
111app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
112app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
113app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
114 \a1 a2 a3 -> f // a1 // a2 // a3 // check s
115app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $
116 \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
117app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $
118 \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s
119app6 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 $
120 \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s
121app7 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 $
122 \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s
123app8 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 $
124 \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s
125app9 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 $
126 \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s
127app10 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 $
128 \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s
129
130
131
132-- GSL error codes are <= 1024
133-- | error codes for the auxiliary functions required by the wrappers
134errorCode :: CInt -> String
135errorCode 2000 = "bad size"
136errorCode 2001 = "bad function code"
137errorCode 2002 = "memory problem"
138errorCode 2003 = "bad file"
139errorCode 2004 = "singular"
140errorCode 2005 = "didn't converge"
141errorCode 2006 = "the input matrix is not positive definite"
142errorCode 2007 = "not yet supported in this OS"
143errorCode n = "code "++show n
144
145
146-- | clear the fpu
147foreign import ccall unsafe "asm_finit" finit :: IO ()
148
149-- | check the error code
150check :: String -> IO CInt -> IO ()
151check msg f = do
152#if FINIT
153 finit
154#endif
155 err <- f
156 when (err/=0) $ if err > 1024
157 then (error (msg++": "++errorCode err)) -- our errors
158 else do -- GSL errors
159 ps <- gsl_strerror err
160 s <- peekCString ps
161 error (msg++": "++s)
162 return ()
163
164-- | description of GSL error codes
165foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
166
167-- | Error capture and conversion to Maybe
168mbCatch :: IO x -> IO (Maybe x)
169mbCatch act = E.catch (Just `fmap` act) f
170 where f :: SomeException -> IO (Maybe x)
171 f _ = return Nothing
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
new file mode 100644
index 0000000..9719fc0
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Matrix.hs
@@ -0,0 +1,491 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE BangPatterns #-}
5-----------------------------------------------------------------------------
6-- |
7-- Module : Data.Packed.Internal.Matrix
8-- Copyright : (c) Alberto Ruiz 2007
9-- License : GPL-style
10--
11-- Maintainer : Alberto Ruiz <aruiz@um.es>
12-- Stability : provisional
13-- Portability : portable (uses FFI)
14--
15-- Internal matrix representation
16--
17-----------------------------------------------------------------------------
18-- #hide
19
20module Data.Packed.Internal.Matrix(
21 Matrix(..), rows, cols, cdat, fdat,
22 MatrixOrder(..), orderOf,
23 createMatrix, mat,
24 cmat, fmat,
25 toLists, flatten, reshape,
26 Element(..),
27 trans,
28 fromRows, toRows, fromColumns, toColumns,
29 matrixFromVector,
30 subMatrix,
31 liftMatrix, liftMatrix2,
32 (@@>), atM',
33 saveMatrix,
34 singleton,
35 emptyM,
36 size, shSize, conformVs, conformMs, conformVTo, conformMTo
37) where
38
39import Data.Packed.Internal.Common
40import Data.Packed.Internal.Signatures
41import Data.Packed.Internal.Vector
42
43import Foreign.Marshal.Alloc(alloca, free)
44import Foreign.Marshal.Array(newArray)
45import Foreign.Ptr(Ptr, castPtr)
46import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
47import Data.Complex(Complex)
48import Foreign.C.Types
49import Foreign.C.String(newCString)
50import System.IO.Unsafe(unsafePerformIO)
51import Control.DeepSeq
52
53-----------------------------------------------------------------
54
55{- Design considerations for the Matrix Type
56 -----------------------------------------
57
58- we must easily handle both row major and column major order,
59 for bindings to LAPACK and GSL/C
60
61- we'd like to simplify redundant matrix transposes:
62 - Some of them arise from the order requirements of some functions
63 - some functions (matrix product) admit transposed arguments
64
65- maybe we don't really need this kind of simplification:
66 - more complex code
67 - some computational overhead
68 - only appreciable gain in code with a lot of redundant transpositions
69 and cheap matrix computations
70
71- we could carry both the matrix and its (lazily computed) transpose.
72 This may save some transpositions, but it is necessary to keep track of the
73 data which is actually computed to be used by functions like the matrix product
74 which admit both orders.
75
76- but if we need the transposed data and it is not in the structure, we must make
77 sure that we touch the same foreignptr that is used in the computation.
78
79- a reasonable solution is using two constructors for a matrix. Transposition just
80 "flips" the constructor. Actual data transposition is not done if followed by a
81 matrix product or another transpose.
82
83-}
84
85data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
86
87transOrder RowMajor = ColumnMajor
88transOrder ColumnMajor = RowMajor
89{- | Matrix representation suitable for GSL and LAPACK computations.
90
91The elements are stored in a continuous memory array.
92
93-}
94
95data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int
96 , icols :: {-# UNPACK #-} !Int
97 , xdat :: {-# UNPACK #-} !(Vector t)
98 , order :: !MatrixOrder }
99-- RowMajor: preferred by C, fdat may require a transposition
100-- ColumnMajor: preferred by LAPACK, cdat may require a transposition
101
102cdat = xdat
103fdat = xdat
104
105rows :: Matrix t -> Int
106rows = irows
107
108cols :: Matrix t -> Int
109cols = icols
110
111orderOf :: Matrix t -> MatrixOrder
112orderOf = order
113
114
115-- | Matrix transpose.
116trans :: Matrix t -> Matrix t
117trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o}
118
119cmat :: (Element t) => Matrix t -> Matrix t
120cmat m@Matrix{order = RowMajor} = m
121cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor}
122
123fmat :: (Element t) => Matrix t -> Matrix t
124fmat m@Matrix{order = ColumnMajor} = m
125fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor}
126
127-- C-Haskell matrix adapter
128-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
129
130mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
131mat a f =
132 unsafeWith (xdat a) $ \p -> do
133 let m g = do
134 g (fi (rows a)) (fi (cols a)) p
135 f m
136
137{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
138
139>>> flatten (ident 3)
140fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
141
142-}
143flatten :: Element t => Matrix t -> Vector t
144flatten = xdat . cmat
145
146{-
147type Mt t s = Int -> Int -> Ptr t -> s
148
149infixr 6 ::>
150type t ::> s = Mt t s
151-}
152
153-- | the inverse of 'Data.Packed.Matrix.fromLists'
154toLists :: (Element t) => Matrix t -> [[t]]
155toLists m = splitEvery (cols m) . toList . flatten $ m
156
157-- | Create a matrix from a list of vectors.
158-- All vectors must have the same dimension,
159-- or dimension 1, which is are automatically expanded.
160fromRows :: Element t => [Vector t] -> Matrix t
161fromRows [] = emptyM 0 0
162fromRows vs = case compatdim (map dim vs) of
163 Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
164 Just 0 -> emptyM r 0
165 Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
166 where
167 r = length vs
168 adapt c v
169 | c == 0 = fromList[]
170 | dim v == c = v
171 | otherwise = constantD (v@>0) c
172
173-- | extracts the rows of a matrix as a list of vectors
174toRows :: Element t => Matrix t -> [Vector t]
175toRows m
176 | c == 0 = replicate r (fromList[])
177 | otherwise = toRows' 0
178 where
179 v = flatten m
180 r = rows m
181 c = cols m
182 toRows' k | k == r*c = []
183 | otherwise = subVector k c v : toRows' (k+c)
184
185-- | Creates a matrix from a list of vectors, as columns
186fromColumns :: Element t => [Vector t] -> Matrix t
187fromColumns m = trans . fromRows $ m
188
189-- | Creates a list of vectors from the columns of a matrix
190toColumns :: Element t => Matrix t -> [Vector t]
191toColumns m = toRows . trans $ m
192
193-- | Reads a matrix position.
194(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
195infixl 9 @@>
196m@Matrix {irows = r, icols = c} @@> (i,j)
197 | safe = if i<0 || i>=r || j<0 || j>=c
198 then error "matrix indexing out of range"
199 else atM' m i j
200 | otherwise = atM' m i j
201{-# INLINE (@@>) #-}
202
203-- Unsafe matrix access without range checking
204atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j)
205atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
206{-# INLINE atM' #-}
207
208------------------------------------------------------------------
209
210matrixFromVector o r c v
211 | r * c == dim v = m
212 | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
213 where
214 m = Matrix { irows = r, icols = c, xdat = v, order = o }
215
216-- allocates memory for a new matrix
217createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
218createMatrix ord r c = do
219 p <- createVector (r*c)
220 return (matrixFromVector ord r c p)
221
222{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
223where r is the desired number of rows.)
224
225>>> reshape 4 (fromList [1..12])
226(3><4)
227 [ 1.0, 2.0, 3.0, 4.0
228 , 5.0, 6.0, 7.0, 8.0
229 , 9.0, 10.0, 11.0, 12.0 ]
230
231-}
232reshape :: Storable t => Int -> Vector t -> Matrix t
233reshape 0 v = matrixFromVector RowMajor 0 0 v
234reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
235
236singleton x = reshape 1 (fromList [x])
237
238-- | application of a vector function on the flattened matrix elements
239liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
240liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d)
241
242-- | application of a vector function on the flattened matrices elements
243liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
244liftMatrix2 f m1 m2
245 | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
246 | otherwise = case orderOf m1 of
247 RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2))
248 ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2))
249
250
251compat :: Matrix a -> Matrix b -> Bool
252compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
253
254------------------------------------------------------------------
255
256{- | Supported matrix elements.
257
258 This class provides optimized internal
259 operations for selected element types.
260 It provides unoptimised defaults for any 'Storable' type,
261 so you can create instances simply as:
262 @instance Element Foo@.
263-}
264class (Storable a) => Element a where
265 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
266 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
267 -> Matrix a -> Matrix a
268 subMatrixD = subMatrix'
269 transdata :: Int -> Vector a -> Int -> Vector a
270 transdata = transdataP -- transdata'
271 constantD :: a -> Int -> Vector a
272 constantD = constantP -- constant'
273
274
275instance Element Float where
276 transdata = transdataAux ctransF
277 constantD = constantAux cconstantF
278
279instance Element Double where
280 transdata = transdataAux ctransR
281 constantD = constantAux cconstantR
282
283instance Element (Complex Float) where
284 transdata = transdataAux ctransQ
285 constantD = constantAux cconstantQ
286
287instance Element (Complex Double) where
288 transdata = transdataAux ctransC
289 constantD = constantAux cconstantC
290
291-------------------------------------------------------------------
292
293transdata' :: Storable a => Int -> Vector a -> Int -> Vector a
294transdata' c1 v c2 =
295 if noneed
296 then v
297 else unsafePerformIO $ do
298 w <- createVector (r2*c2)
299 unsafeWith v $ \p ->
300 unsafeWith w $ \q -> do
301 let go (-1) _ = return ()
302 go !i (-1) = go (i-1) (c1-1)
303 go !i !j = do x <- peekElemOff p (i*c1+j)
304 pokeElemOff q (j*c2+i) x
305 go i (j-1)
306 go (r1-1) (c1-1)
307 return w
308 where r1 = dim v `div` c1
309 r2 = dim v `div` c2
310 noneed = dim v == 0 || r1 == 1 || c1 == 1
311
312-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int -> Vector Double #-}
313-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-}
314
315-- I don't know how to specialize...
316-- The above pragmas only seem to work on top level defs
317-- Fortunately everything seems to work using the above class
318
319-- C versions, still a little faster:
320
321transdataAux fun c1 d c2 =
322 if noneed
323 then d
324 else unsafePerformIO $ do
325 v <- createVector (dim d)
326 unsafeWith d $ \pd ->
327 unsafeWith v $ \pv ->
328 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
329 return v
330 where r1 = dim d `div` c1
331 r2 = dim d `div` c2
332 noneed = dim d == 0 || r1 == 1 || c1 == 1
333
334transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
335transdataP c1 d c2 =
336 if noneed
337 then d
338 else unsafePerformIO $ do
339 v <- createVector (dim d)
340 unsafeWith d $ \pd ->
341 unsafeWith v $ \pv ->
342 ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"
343 return v
344 where r1 = dim d `div` c1
345 r2 = dim d `div` c2
346 sz = sizeOf (d @> 0)
347 noneed = dim d == 0 || r1 == 1 || c1 == 1
348
349foreign import ccall unsafe "transF" ctransF :: TFMFM
350foreign import ccall unsafe "transR" ctransR :: TMM
351foreign import ccall unsafe "transQ" ctransQ :: TQMQM
352foreign import ccall unsafe "transC" ctransC :: TCMCM
353foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
354
355----------------------------------------------------------------------
356
357constant' v n = unsafePerformIO $ do
358 w <- createVector n
359 unsafeWith w $ \p -> do
360 let go (-1) = return ()
361 go !k = pokeElemOff p k v >> go (k-1)
362 go (n-1)
363 return w
364
365-- C versions
366
367constantAux fun x n = unsafePerformIO $ do
368 v <- createVector n
369 px <- newArray [x]
370 app1 (fun px) vec v "constantAux"
371 free px
372 return v
373
374constantF :: Float -> Int -> Vector Float
375constantF = constantAux cconstantF
376foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
377
378constantR :: Double -> Int -> Vector Double
379constantR = constantAux cconstantR
380foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
381
382constantQ :: Complex Float -> Int -> Vector (Complex Float)
383constantQ = constantAux cconstantQ
384foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
385
386constantC :: Complex Double -> Int -> Vector (Complex Double)
387constantC = constantAux cconstantC
388foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
389
390constantP :: Storable a => a -> Int -> Vector a
391constantP a n = unsafePerformIO $ do
392 let sz = sizeOf a
393 v <- createVector n
394 unsafeWith v $ \p -> do
395 alloca $ \k -> do
396 poke k a
397 cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
398 return v
399foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
400
401----------------------------------------------------------------------
402
403-- | Extracts a submatrix from a matrix.
404subMatrix :: Element a
405 => (Int,Int) -- ^ (r0,c0) starting position
406 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
407 -> Matrix a -- ^ input matrix
408 -> Matrix a -- ^ result
409subMatrix (r0,c0) (rt,ct) m
410 | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
411 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
412 | otherwise = error $ "wrong subMatrix "++
413 show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
414
415subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
416 w <- createVector (rt*ct)
417 unsafeWith v $ \p ->
418 unsafeWith w $ \q -> do
419 let go (-1) _ = return ()
420 go !i (-1) = go (i-1) (ct-1)
421 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)
422 pokeElemOff q (i*ct+j) x
423 go i (j-1)
424 go (rt-1) (ct-1)
425 return w
426
427subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor
428subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
429
430--------------------------------------------------------------------------
431
432-- | Saves a matrix as 2D ASCII table.
433saveMatrix :: FilePath
434 -> String -- ^ format (%f, %g, %e)
435 -> Matrix Double
436 -> IO ()
437saveMatrix filename fmt m = do
438 charname <- newCString filename
439 charfmt <- newCString fmt
440 let o = if orderOf m == RowMajor then 1 else 0
441 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
442 free charname
443 free charfmt
444
445foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
446
447----------------------------------------------------------------------
448
449maxZ xs = if minimum xs == 0 then 0 else maximum xs
450
451conformMs ms = map (conformMTo (r,c)) ms
452 where
453 r = maxZ (map rows ms)
454 c = maxZ (map cols ms)
455
456
457conformVs vs = map (conformVTo n) vs
458 where
459 n = maxZ (map dim vs)
460
461conformMTo (r,c) m
462 | size m == (r,c) = m
463 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
464 | size m == (r,1) = repCols c m
465 | size m == (1,c) = repRows r m
466 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
467
468conformVTo n v
469 | dim v == n = v
470 | dim v == 1 = constantD (v@>0) n
471 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
472
473repRows n x = fromRows (replicate n (flatten x))
474repCols n x = fromColumns (replicate n (flatten x))
475
476size m = (rows m, cols m)
477
478shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
479
480emptyM r c = matrixFromVector RowMajor r c (fromList[])
481
482----------------------------------------------------------------------
483
484instance (Storable t, NFData t) => NFData (Matrix t)
485 where
486 rnf m | d > 0 = rnf (v @> 0)
487 | otherwise = ()
488 where
489 d = dim v
490 v = xdat m
491
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
new file mode 100644
index 0000000..2835720
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Signatures.hs
@@ -0,0 +1,72 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Internal.Signatures
4-- Copyright : (c) Alberto Ruiz 2009
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9-- Portability : portable (uses FFI)
10--
11-- Signatures of the C functions.
12--
13-----------------------------------------------------------------------------
14
15module Data.Packed.Internal.Signatures where
16
17import Foreign.Ptr(Ptr)
18import Data.Complex(Complex)
19import Foreign.C.Types(CInt)
20
21type PF = Ptr Float --
22type PD = Ptr Double --
23type PQ = Ptr (Complex Float) --
24type PC = Ptr (Complex Double) --
25type TF = CInt -> PF -> IO CInt --
26type TFF = CInt -> PF -> TF --
27type TFV = CInt -> PF -> TV --
28type TVF = CInt -> PD -> TF --
29type TFFF = CInt -> PF -> TFF --
30type TV = CInt -> PD -> IO CInt --
31type TVV = CInt -> PD -> TV --
32type TVVV = CInt -> PD -> TVV --
33type TFM = CInt -> CInt -> PF -> IO CInt --
34type TFMFM = CInt -> CInt -> PF -> TFM --
35type TFMFMFM = CInt -> CInt -> PF -> TFMFM --
36type TM = CInt -> CInt -> PD -> IO CInt --
37type TMM = CInt -> CInt -> PD -> TM --
38type TVMM = CInt -> PD -> TMM --
39type TMVMM = CInt -> CInt -> PD -> TVMM --
40type TMMM = CInt -> CInt -> PD -> TMM --
41type TVM = CInt -> PD -> TM --
42type TVVM = CInt -> PD -> TVM --
43type TMV = CInt -> CInt -> PD -> TV --
44type TMMV = CInt -> CInt -> PD -> TMV --
45type TMVM = CInt -> CInt -> PD -> TVM --
46type TMMVM = CInt -> CInt -> PD -> TMVM --
47type TCM = CInt -> CInt -> PC -> IO CInt --
48type TCVCM = CInt -> PC -> TCM --
49type TCMCVCM = CInt -> CInt -> PC -> TCVCM --
50type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM --
51type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM --
52type TCMCM = CInt -> CInt -> PC -> TCM --
53type TVCM = CInt -> PD -> TCM --
54type TCMVCM = CInt -> CInt -> PC -> TVCM --
55type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM --
56type TCMCMCM = CInt -> CInt -> PC -> TCMCM --
57type TCV = CInt -> PC -> IO CInt --
58type TCVCV = CInt -> PC -> TCV --
59type TCVCVCV = CInt -> PC -> TCVCV --
60type TCVV = CInt -> PC -> TV --
61type TQV = CInt -> PQ -> IO CInt --
62type TQVQV = CInt -> PQ -> TQV --
63type TQVQVQV = CInt -> PQ -> TQVQV --
64type TQVF = CInt -> PQ -> TF --
65type TQM = CInt -> CInt -> PQ -> IO CInt --
66type TQMQM = CInt -> CInt -> PQ -> TQM --
67type TQMQMQM = CInt -> CInt -> PQ -> TQMQM --
68type TCMCV = CInt -> CInt -> PC -> TCV --
69type TVCV = CInt -> PD -> TCV --
70type TCVM = CInt -> PC -> TM --
71type TMCVM = CInt -> CInt -> PD -> TCVM --
72type TMMCVM = CInt -> CInt -> PD -> TMCVM --
diff --git a/packages/hmatrix/src/Data/Packed/Internal/Vector.hs b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
new file mode 100644
index 0000000..6d03438
--- /dev/null
+++ b/packages/hmatrix/src/Data/Packed/Internal/Vector.hs
@@ -0,0 +1,521 @@
1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Internal.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 implementation
13--
14-----------------------------------------------------------------------------
15
16module Data.Packed.Internal.Vector (
17 Vector, dim,
18 fromList, toList, (|>),
19 vjoin, (@>), safe, at, at', subVector, takesV,
20 mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith,
21 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
22 foldVector, foldVectorG, foldLoop, foldVectorWithIndex,
23 createVector, vec,
24 asComplex, asReal, float2DoubleV, double2FloatV,
25 stepF, stepD, condF, condD,
26 conjugateQ, conjugateC,
27 fwriteVector, freadVector, fprintfVector, fscanfVector,
28 cloneVector,
29 unsafeToForeignPtr,
30 unsafeFromForeignPtr,
31 unsafeWith
32) where
33
34import Data.Packed.Internal.Common
35import Data.Packed.Internal.Signatures
36import Foreign.Marshal.Alloc(free)
37import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
38import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
39import Foreign.Ptr(Ptr)
40import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
41import Foreign.C.String
42import Foreign.C.Types
43import Data.Complex
44import Control.Monad(when)
45import System.IO.Unsafe(unsafePerformIO)
46
47#if __GLASGOW_HASKELL__ >= 605
48import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
49#else
50import Foreign.ForeignPtr (mallocForeignPtrBytes)
51#endif
52
53import GHC.Base
54#if __GLASGOW_HASKELL__ < 612
55import GHC.IOBase hiding (liftIO)
56#endif
57
58import qualified Data.Vector.Storable as Vector
59import Data.Vector.Storable(Vector,
60 fromList,
61 unsafeToForeignPtr,
62 unsafeFromForeignPtr,
63 unsafeWith)
64
65
66-- | Number of elements
67dim :: (Storable t) => Vector t -> Int
68dim = Vector.length
69
70
71-- C-Haskell vector adapter
72-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
73vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
74vec x f = unsafeWith x $ \p -> do
75 let v g = do
76 g (fi $ dim x) p
77 f v
78{-# INLINE vec #-}
79
80
81-- allocates memory for a new vector
82createVector :: Storable a => Int -> IO (Vector a)
83createVector n = do
84 when (n < 0) $ error ("trying to createVector of negative dim: "++show n)
85 fp <- doMalloc undefined
86 return $ unsafeFromForeignPtr fp 0 n
87 where
88 --
89 -- Use the much cheaper Haskell heap allocated storage
90 -- for foreign pointer space we control
91 --
92 doMalloc :: Storable b => b -> IO (ForeignPtr b)
93 doMalloc dummy = do
94#if __GLASGOW_HASKELL__ >= 605
95 mallocPlainForeignPtrBytes (n * sizeOf dummy)
96#else
97 mallocForeignPtrBytes (n * sizeOf dummy)
98#endif
99
100{- | creates a Vector from a list:
101
102@> fromList [2,3,5,7]
1034 |> [2.0,3.0,5.0,7.0]@
104
105-}
106
107safeRead v = inlinePerformIO . unsafeWith v
108{-# INLINE safeRead #-}
109
110inlinePerformIO :: IO a -> a
111inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
112{-# INLINE inlinePerformIO #-}
113
114{- | extracts the Vector elements to a list
115
116>>> toList (linspace 5 (1,10))
117[1.0,3.25,5.5,7.75,10.0]
118
119-}
120toList :: Storable a => Vector a -> [a]
121toList v = safeRead v $ peekArray (dim v)
122
123{- | Create a vector from a list of elements and explicit dimension. The input
124 list is explicitly truncated if it is too long, so it may safely
125 be used, for instance, with infinite lists.
126
127>>> 5 |> [1..]
128fromList [1.0,2.0,3.0,4.0,5.0]
129
130-}
131(|>) :: (Storable a) => Int -> [a] -> Vector a
132infixl 9 |>
133n |> l = if length l' == n
134 then fromList l'
135 else error "list too short for |>"
136 where l' = take n l
137
138
139-- | access to Vector elements without range checking
140at' :: Storable a => Vector a -> Int -> a
141at' v n = safeRead v $ flip peekElemOff n
142{-# INLINE at' #-}
143
144--
145-- turn off bounds checking with -funsafe at configure time.
146-- ghc will optimise away the salways true case at compile time.
147--
148#if defined(UNSAFE)
149safe :: Bool
150safe = False
151#else
152safe = True
153#endif
154
155-- | access to Vector elements with range checking.
156at :: Storable a => Vector a -> Int -> a
157at v n
158 | safe = if n >= 0 && n < dim v
159 then at' v n
160 else error "vector index out of range"
161 | otherwise = at' v n
162{-# INLINE at #-}
163
164{- | takes a number of consecutive elements from a Vector
165
166>>> subVector 2 3 (fromList [1..10])
167fromList [3.0,4.0,5.0]
168
169-}
170subVector :: Storable t => Int -- ^ index of the starting element
171 -> Int -- ^ number of elements to extract
172 -> Vector t -- ^ source
173 -> Vector t -- ^ result
174subVector = Vector.slice
175
176
177{- | Reads a vector position:
178
179>>> fromList [0..9] @> 7
1807.0
181
182-}
183(@>) :: Storable t => Vector t -> Int -> t
184infixl 9 @>
185(@>) = at
186
187
188{- | concatenate a list of vectors
189
190>>> vjoin [fromList [1..5::Double], konst 1 3]
191fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
192
193-}
194vjoin :: Storable t => [Vector t] -> Vector t
195vjoin [] = fromList []
196vjoin [v] = v
197vjoin as = unsafePerformIO $ do
198 let tot = sum (map dim as)
199 r <- createVector tot
200 unsafeWith r $ \ptr ->
201 joiner as tot ptr
202 return r
203 where joiner [] _ _ = return ()
204 joiner (v:cs) _ p = do
205 let n = dim v
206 unsafeWith v $ \pb -> copyArray p pb n
207 joiner cs 0 (advancePtr p n)
208
209
210{- | Extract consecutive subvectors of the given sizes.
211
212>>> takesV [3,4] (linspace 10 (1,10::Double))
213[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
214
215-}
216takesV :: Storable t => [Int] -> Vector t -> [Vector t]
217takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w)
218 | otherwise = go ms w
219 where go [] _ = []
220 go (n:ns) v = subVector 0 n v
221 : go ns (subVector n (dim v - n) v)
222
223---------------------------------------------------------------
224
225-- | transforms a complex vector into a real vector with alternating real and imaginary parts
226asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
227asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
228 where (fp,i,n) = unsafeToForeignPtr v
229
230-- | transforms a real vector into a complex vector with alternating real and imaginary parts
231asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
232asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
233 where (fp,i,n) = unsafeToForeignPtr v
234
235---------------------------------------------------------------
236
237float2DoubleV :: Vector Float -> Vector Double
238float2DoubleV v = unsafePerformIO $ do
239 r <- createVector (dim v)
240 app2 c_float2double vec v vec r "float2double"
241 return r
242
243double2FloatV :: Vector Double -> Vector Float
244double2FloatV v = unsafePerformIO $ do
245 r <- createVector (dim v)
246 app2 c_double2float vec v vec r "double2float2"
247 return r
248
249
250foreign import ccall unsafe "float2double" c_float2double:: TFV
251foreign import ccall unsafe "double2float" c_double2float:: TVF
252
253---------------------------------------------------------------
254
255stepF :: Vector Float -> Vector Float
256stepF v = unsafePerformIO $ do
257 r <- createVector (dim v)
258 app2 c_stepF vec v vec r "stepF"
259 return r
260
261stepD :: Vector Double -> Vector Double
262stepD v = unsafePerformIO $ do
263 r <- createVector (dim v)
264 app2 c_stepD vec v vec r "stepD"
265 return r
266
267foreign import ccall unsafe "stepF" c_stepF :: TFF
268foreign import ccall unsafe "stepD" c_stepD :: TVV
269
270---------------------------------------------------------------
271
272condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
273condF x y l e g = unsafePerformIO $ do
274 r <- createVector (dim x)
275 app6 c_condF vec x vec y vec l vec e vec g vec r "condF"
276 return r
277
278condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
279condD x y l e g = unsafePerformIO $ do
280 r <- createVector (dim x)
281 app6 c_condD vec x vec y vec l vec e vec g vec r "condD"
282 return r
283
284foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF
285foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV
286
287--------------------------------------------------------------------------------
288
289conjugateAux fun x = unsafePerformIO $ do
290 v <- createVector (dim x)
291 app2 fun vec x vec v "conjugateAux"
292 return v
293
294conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
295conjugateQ = conjugateAux c_conjugateQ
296foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV
297
298conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
299conjugateC = conjugateAux c_conjugateC
300foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV
301
302--------------------------------------------------------------------------------
303
304cloneVector :: Storable t => Vector t -> IO (Vector t)
305cloneVector v = do
306 let n = dim v
307 r <- createVector n
308 let f _ s _ d = copyArray d s n >> return 0
309 app2 f vec v vec r "cloneVector"
310 return r
311
312------------------------------------------------------------------
313
314-- | map on Vectors
315mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
316mapVector f v = unsafePerformIO $ do
317 w <- createVector (dim v)
318 unsafeWith v $ \p ->
319 unsafeWith w $ \q -> do
320 let go (-1) = return ()
321 go !k = do x <- peekElemOff p k
322 pokeElemOff q k (f x)
323 go (k-1)
324 go (dim v -1)
325 return w
326{-# INLINE mapVector #-}
327
328-- | zipWith for Vectors
329zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
330zipVectorWith f u v = unsafePerformIO $ do
331 let n = min (dim u) (dim v)
332 w <- createVector n
333 unsafeWith u $ \pu ->
334 unsafeWith v $ \pv ->
335 unsafeWith w $ \pw -> do
336 let go (-1) = return ()
337 go !k = do x <- peekElemOff pu k
338 y <- peekElemOff pv k
339 pokeElemOff pw k (f x y)
340 go (k-1)
341 go (n -1)
342 return w
343{-# INLINE zipVectorWith #-}
344
345-- | unzipWith for Vectors
346unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
347 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
348unzipVectorWith f u = unsafePerformIO $ do
349 let n = dim u
350 v <- createVector n
351 w <- createVector n
352 unsafeWith u $ \pu ->
353 unsafeWith v $ \pv ->
354 unsafeWith w $ \pw -> do
355 let go (-1) = return ()
356 go !k = do z <- peekElemOff pu k
357 let (x,y) = f z
358 pokeElemOff pv k x
359 pokeElemOff pw k y
360 go (k-1)
361 go (n-1)
362 return (v,w)
363{-# INLINE unzipVectorWith #-}
364
365foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
366foldVector f x v = unsafePerformIO $
367 unsafeWith v $ \p -> do
368 let go (-1) s = return s
369 go !k !s = do y <- peekElemOff p k
370 go (k-1::Int) (f y s)
371 go (dim v -1) x
372{-# INLINE foldVector #-}
373
374-- the zero-indexed index is passed to the folding function
375foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
376foldVectorWithIndex f x v = unsafePerformIO $
377 unsafeWith v $ \p -> do
378 let go (-1) s = return s
379 go !k !s = do y <- peekElemOff p k
380 go (k-1::Int) (f k y s)
381 go (dim v -1) x
382{-# INLINE foldVectorWithIndex #-}
383
384foldLoop f s0 d = go (d - 1) s0
385 where
386 go 0 s = f (0::Int) s
387 go !j !s = go (j - 1) (f j s)
388
389foldVectorG f s0 v = foldLoop g s0 (dim v)
390 where g !k !s = f k (at' v) s
391 {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
392{-# INLINE foldVectorG #-}
393
394-------------------------------------------------------------------
395
396-- | monadic map over Vectors
397-- the monad @m@ must be strict
398mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
399mapVectorM f v = do
400 w <- return $! unsafePerformIO $! createVector (dim v)
401 mapVectorM' w 0 (dim v -1)
402 return w
403 where mapVectorM' w' !k !t
404 | k == t = do
405 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
406 y <- f x
407 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
408 | otherwise = do
409 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
410 y <- f x
411 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
412 mapVectorM' w' (k+1) t
413{-# INLINE mapVectorM #-}
414
415-- | monadic map over Vectors
416mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
417mapVectorM_ f v = do
418 mapVectorM' 0 (dim v -1)
419 where mapVectorM' !k !t
420 | k == t = do
421 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
422 f x
423 | otherwise = do
424 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
425 _ <- f x
426 mapVectorM' (k+1) t
427{-# INLINE mapVectorM_ #-}
428
429-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
430-- the monad @m@ must be strict
431mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
432mapVectorWithIndexM f v = do
433 w <- return $! unsafePerformIO $! createVector (dim v)
434 mapVectorM' w 0 (dim v -1)
435 return w
436 where mapVectorM' w' !k !t
437 | k == t = do
438 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
439 y <- f k x
440 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
441 | otherwise = do
442 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
443 y <- f k x
444 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
445 mapVectorM' w' (k+1) t
446{-# INLINE mapVectorWithIndexM #-}
447
448-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
449mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
450mapVectorWithIndexM_ f v = do
451 mapVectorM' 0 (dim v -1)
452 where mapVectorM' !k !t
453 | k == t = do
454 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
455 f k x
456 | otherwise = do
457 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
458 _ <- f k x
459 mapVectorM' (k+1) t
460{-# INLINE mapVectorWithIndexM_ #-}
461
462
463mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
464--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
465mapVectorWithIndex f v = unsafePerformIO $ do
466 w <- createVector (dim v)
467 unsafeWith v $ \p ->
468 unsafeWith w $ \q -> do
469 let go (-1) = return ()
470 go !k = do x <- peekElemOff p k
471 pokeElemOff q k (f k x)
472 go (k-1)
473 go (dim v -1)
474 return w
475{-# INLINE mapVectorWithIndex #-}
476
477-------------------------------------------------------------------
478
479
480-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
481fscanfVector :: FilePath -> Int -> IO (Vector Double)
482fscanfVector filename n = do
483 charname <- newCString filename
484 res <- createVector n
485 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
486 free charname
487 return res
488
489foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
490
491-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
492fprintfVector :: FilePath -> String -> Vector Double -> IO ()
493fprintfVector filename fmt v = do
494 charname <- newCString filename
495 charfmt <- newCString fmt
496 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
497 free charname
498 free charfmt
499
500foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
501
502-- | Loads a vector from a binary file (the number of elements must be known in advance).
503freadVector :: FilePath -> Int -> IO (Vector Double)
504freadVector filename n = do
505 charname <- newCString filename
506 res <- createVector n
507 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
508 free charname
509 return res
510
511foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
512
513-- | Saves the elements of a vector to a binary file.
514fwriteVector :: FilePath -> Vector Double -> IO ()
515fwriteVector filename v = do
516 charname <- newCString filename
517 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
518 free charname
519
520foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
521