summaryrefslogtreecommitdiff
path: root/lib/Data/Packed/Internal
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Data/Packed/Internal')
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs40
-rw-r--r--lib/Data/Packed/Internal/Vector.hs18
2 files changed, 34 insertions, 24 deletions
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 33c9324..8bd6fb2 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -16,7 +16,7 @@
16-- #hide 16-- #hide
17 17
18module Data.Packed.Internal.Matrix( 18module Data.Packed.Internal.Matrix(
19 Matrix(..), 19 Matrix(..), rows, cols,
20 MatrixOrder(..), orderOf, 20 MatrixOrder(..), orderOf,
21 createMatrix, withMatrix, mat, 21 createMatrix, withMatrix, mat,
22 cmat, fmat, 22 cmat, fmat,
@@ -77,17 +77,23 @@ import Foreign.C.String
77data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) 77data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
78 78
79-- | Matrix representation suitable for GSL and LAPACK computations. 79-- | Matrix representation suitable for GSL and LAPACK computations.
80data Matrix t = MC { rows :: {-# UNPACK #-} !Int 80data Matrix t = MC { irows :: {-# UNPACK #-} !Int
81 , cols :: {-# UNPACK #-} !Int 81 , icols :: {-# UNPACK #-} !Int
82 , cdat :: {-# UNPACK #-} !(Vector t) } 82 , cdat :: {-# UNPACK #-} !(Vector t) }
83 83
84 | MF { rows :: {-# UNPACK #-} !Int 84 | MF { irows :: {-# UNPACK #-} !Int
85 , cols :: {-# UNPACK #-} !Int 85 , icols :: {-# UNPACK #-} !Int
86 , fdat :: {-# UNPACK #-} !(Vector t) } 86 , fdat :: {-# UNPACK #-} !(Vector t) }
87 87
88-- MC: preferred by C, fdat may require a transposition 88-- MC: preferred by C, fdat may require a transposition
89-- MF: preferred by LAPACK, cdat may require a transposition 89-- MF: preferred by LAPACK, cdat may require a transposition
90 90
91rows :: Matrix t -> Int
92rows = irows
93
94cols :: Matrix t -> Int
95cols = icols
96
91xdat MC {cdat = d } = d 97xdat MC {cdat = d } = d
92xdat MF {fdat = d } = d 98xdat MF {fdat = d } = d
93 99
@@ -97,16 +103,16 @@ orderOf MC{} = RowMajor
97 103
98-- | Matrix transpose. 104-- | Matrix transpose.
99trans :: Matrix t -> Matrix t 105trans :: Matrix t -> Matrix t
100trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d } 106trans MC {irows = r, icols = c, cdat = d } = MF {irows = c, icols = r, fdat = d }
101trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d } 107trans MF {irows = r, icols = c, fdat = d } = MC {irows = c, icols = r, cdat = d }
102 108
103cmat :: (Element t) => Matrix t -> Matrix t 109cmat :: (Element t) => Matrix t -> Matrix t
104cmat m@MC{} = m 110cmat m@MC{} = m
105cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c} 111cmat MF {irows = r, icols = c, fdat = d } = MC {irows = r, icols = c, cdat = transdata r d c}
106 112
107fmat :: (Element t) => Matrix t -> Matrix t 113fmat :: (Element t) => Matrix t -> Matrix t
108fmat m@MF{} = m 114fmat m@MF{} = m
109fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} 115fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r}
110 116
111-- C-Haskell matrix adapter 117-- C-Haskell matrix adapter
112mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r 118mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
@@ -170,13 +176,13 @@ infixl 9 @@>
170-- | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range" 176-- | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
171-- | otherwise = cdat m `at` (i*c+j) 177-- | otherwise = cdat m `at` (i*c+j)
172 178
173MC {rows = r, cols = c, cdat = v} @@> (i,j) 179MC {irows = r, icols = c, cdat = v} @@> (i,j)
174 | safe = if i<0 || i>=r || j<0 || j>=c 180 | safe = if i<0 || i>=r || j<0 || j>=c
175 then error "matrix indexing out of range" 181 then error "matrix indexing out of range"
176 else v `at` (i*c+j) 182 else v `at` (i*c+j)
177 | otherwise = v `at` (i*c+j) 183 | otherwise = v `at` (i*c+j)
178 184
179MF {rows = r, cols = c, fdat = v} @@> (i,j) 185MF {irows = r, icols = c, fdat = v} @@> (i,j)
180 | safe = if i<0 || i>=r || j<0 || j>=c 186 | safe = if i<0 || i>=r || j<0 || j>=c
181 then error "matrix indexing out of range" 187 then error "matrix indexing out of range"
182 else v `at` (j*r+i) 188 else v `at` (j*r+i)
@@ -184,18 +190,18 @@ MF {rows = r, cols = c, fdat = v} @@> (i,j)
184{-# INLINE (@@>) #-} 190{-# INLINE (@@>) #-}
185 191
186-- Unsafe matrix access without range checking 192-- Unsafe matrix access without range checking
187atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j) 193atM' MC {icols = c, cdat = v} i j = v `at'` (i*c+j)
188atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i) 194atM' MF {irows = r, fdat = v} i j = v `at'` (j*r+i)
189{-# INLINE atM' #-} 195{-# INLINE atM' #-}
190 196
191------------------------------------------------------------------ 197------------------------------------------------------------------
192 198
193matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v } 199matrixFromVector RowMajor c v = MC { irows = r, icols = c, cdat = v }
194 where (d,m) = dim v `divMod` c 200 where (d,m) = dim v `divMod` c
195 r | m==0 = d 201 r | m==0 = d
196 | otherwise = error "matrixFromVector" 202 | otherwise = error "matrixFromVector"
197 203
198matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v } 204matrixFromVector ColumnMajor c v = MF { irows = r, icols = c, fdat = v }
199 where (d,m) = dim v `divMod` c 205 where (d,m) = dim v `divMod` c
200 r | m==0 = d 206 r | m==0 = d
201 | otherwise = error "matrixFromVector" 207 | otherwise = error "matrixFromVector"
@@ -223,8 +229,8 @@ singleton x = reshape 1 (fromList [x])
223 229
224-- | application of a vector function on the flattened matrix elements 230-- | application of a vector function on the flattened matrix elements
225liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b 231liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
226liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d) 232liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor c (f d)
227liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) 233liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)
228 234
229-- | application of a vector function on the flattened matrices elements 235-- | application of a vector function on the flattened matrices elements
230liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t 236liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index bc623de..66acd87 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -15,7 +15,7 @@
15-- #hide 15-- #hide
16 16
17module Data.Packed.Internal.Vector ( 17module Data.Packed.Internal.Vector (
18 Vector(..), 18 Vector(..), dim,
19 fromList, toList, (|>), 19 fromList, toList, (|>),
20 join, (@>), safe, at, at', subVector, 20 join, (@>), safe, at, at', subVector,
21 mapVector, zipVector, 21 mapVector, zipVector,
@@ -47,10 +47,14 @@ import GHC.IOBase
47 47
48-- | A one-dimensional array of objects stored in a contiguous memory block. 48-- | A one-dimensional array of objects stored in a contiguous memory block.
49data Vector t = 49data Vector t =
50 V { dim :: {-# UNPACK #-} !Int -- ^ number of elements 50 V { idim :: {-# UNPACK #-} !Int -- ^ number of elements
51 , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block 51 , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block
52 } 52 }
53 53
54-- | Number of elements
55dim :: Vector t -> Int
56dim = idim
57
54-- C-Haskell vector adapter 58-- C-Haskell vector adapter
55vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r 59vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
56vec = withVector 60vec = withVector
@@ -157,7 +161,7 @@ subVector :: Storable t => Int -- ^ index of the starting element
157 -> Int -- ^ number of elements to extract 161 -> Int -- ^ number of elements to extract
158 -> Vector t -- ^ source 162 -> Vector t -- ^ source
159 -> Vector t -- ^ result 163 -> Vector t -- ^ result
160subVector k l (v@V {dim=n}) 164subVector k l (v@V {idim=n})
161 | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" 165 | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
162 | otherwise = unsafePerformIO $ do 166 | otherwise = unsafePerformIO $ do
163 r <- createVector l 167 r <- createVector l
@@ -192,23 +196,23 @@ join as = unsafePerformIO $ do
192 joiner as tot ptr 196 joiner as tot ptr
193 return r 197 return r
194 where joiner [] _ _ = return () 198 where joiner [] _ _ = return ()
195 joiner (V {dim = n, fptr = b} : cs) _ p = do 199 joiner (V {idim = n, fptr = b} : cs) _ p = do
196 withForeignPtr b $ \pb -> copyArray p pb n 200 withForeignPtr b $ \pb -> copyArray p pb n
197 joiner cs 0 (advancePtr p n) 201 joiner cs 0 (advancePtr p n)
198 202
199 203
200-- | transforms a complex vector into a real vector with alternating real and imaginary parts 204-- | transforms a complex vector into a real vector with alternating real and imaginary parts
201asReal :: Vector (Complex Double) -> Vector Double 205asReal :: Vector (Complex Double) -> Vector Double
202asReal v = V { dim = 2*dim v, fptr = castForeignPtr (fptr v) } 206asReal v = V { idim = 2*dim v, fptr = castForeignPtr (fptr v) }
203 207
204-- | transforms a real vector into a complex vector with alternating real and imaginary parts 208-- | transforms a real vector into a complex vector with alternating real and imaginary parts
205asComplex :: Vector Double -> Vector (Complex Double) 209asComplex :: Vector Double -> Vector (Complex Double)
206asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v) } 210asComplex v = V { idim = dim v `div` 2, fptr = castForeignPtr (fptr v) }
207 211
208---------------------------------------------------------------- 212----------------------------------------------------------------
209 213
210cloneVector :: Storable t => Vector t -> IO (Vector t) 214cloneVector :: Storable t => Vector t -> IO (Vector t)
211cloneVector (v@V {dim=n}) = do 215cloneVector (v@V {idim=n}) = do
212 r <- createVector n 216 r <- createVector n
213 let f _ s _ d = copyArray d s n >> return 0 217 let f _ s _ d = copyArray d s n >> return 0
214 app2 f vec v vec r "cloneVector" 218 app2 f vec v vec r "cloneVector"