diff options
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 40 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 18 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Instances.hs | 2 |
3 files changed, 35 insertions, 25 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 | ||
18 | module Data.Packed.Internal.Matrix( | 18 | module 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 | |||
77 | data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) | 77 | data 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. |
80 | data Matrix t = MC { rows :: {-# UNPACK #-} !Int | 80 | data 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 | ||
91 | rows :: Matrix t -> Int | ||
92 | rows = irows | ||
93 | |||
94 | cols :: Matrix t -> Int | ||
95 | cols = icols | ||
96 | |||
91 | xdat MC {cdat = d } = d | 97 | xdat MC {cdat = d } = d |
92 | xdat MF {fdat = d } = d | 98 | xdat MF {fdat = d } = d |
93 | 99 | ||
@@ -97,16 +103,16 @@ orderOf MC{} = RowMajor | |||
97 | 103 | ||
98 | -- | Matrix transpose. | 104 | -- | Matrix transpose. |
99 | trans :: Matrix t -> Matrix t | 105 | trans :: Matrix t -> Matrix t |
100 | trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d } | 106 | trans MC {irows = r, icols = c, cdat = d } = MF {irows = c, icols = r, fdat = d } |
101 | trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d } | 107 | trans MF {irows = r, icols = c, fdat = d } = MC {irows = c, icols = r, cdat = d } |
102 | 108 | ||
103 | cmat :: (Element t) => Matrix t -> Matrix t | 109 | cmat :: (Element t) => Matrix t -> Matrix t |
104 | cmat m@MC{} = m | 110 | cmat m@MC{} = m |
105 | cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c} | 111 | cmat MF {irows = r, icols = c, fdat = d } = MC {irows = r, icols = c, cdat = transdata r d c} |
106 | 112 | ||
107 | fmat :: (Element t) => Matrix t -> Matrix t | 113 | fmat :: (Element t) => Matrix t -> Matrix t |
108 | fmat m@MF{} = m | 114 | fmat m@MF{} = m |
109 | fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} | 115 | fmat 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 |
112 | mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r | 118 | mat :: 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 | ||
173 | MC {rows = r, cols = c, cdat = v} @@> (i,j) | 179 | MC {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 | ||
179 | MF {rows = r, cols = c, fdat = v} @@> (i,j) | 185 | MF {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 |
187 | atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j) | 193 | atM' MC {icols = c, cdat = v} i j = v `at'` (i*c+j) |
188 | atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i) | 194 | atM' MF {irows = r, fdat = v} i j = v `at'` (j*r+i) |
189 | {-# INLINE atM' #-} | 195 | {-# INLINE atM' #-} |
190 | 196 | ||
191 | ------------------------------------------------------------------ | 197 | ------------------------------------------------------------------ |
192 | 198 | ||
193 | matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v } | 199 | matrixFromVector 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 | ||
198 | matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v } | 204 | matrixFromVector 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 |
225 | liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b | 231 | liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b |
226 | liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d) | 232 | liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor c (f d) |
227 | liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) | 233 | liftMatrix 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 |
230 | liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t | 236 | liftMatrix2 :: (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 | ||
17 | module Data.Packed.Internal.Vector ( | 17 | module 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. |
49 | data Vector t = | 49 | data 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 | ||
55 | dim :: Vector t -> Int | ||
56 | dim = idim | ||
57 | |||
54 | -- C-Haskell vector adapter | 58 | -- C-Haskell vector adapter |
55 | vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r | 59 | vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r |
56 | vec = withVector | 60 | vec = 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 |
160 | subVector k l (v@V {dim=n}) | 164 | subVector 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 |
201 | asReal :: Vector (Complex Double) -> Vector Double | 205 | asReal :: Vector (Complex Double) -> Vector Double |
202 | asReal v = V { dim = 2*dim v, fptr = castForeignPtr (fptr v) } | 206 | asReal 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 |
205 | asComplex :: Vector Double -> Vector (Complex Double) | 209 | asComplex :: Vector Double -> Vector (Complex Double) |
206 | asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v) } | 210 | asComplex v = V { idim = dim v `div` 2, fptr = castForeignPtr (fptr v) } |
207 | 211 | ||
208 | ---------------------------------------------------------------- | 212 | ---------------------------------------------------------------- |
209 | 213 | ||
210 | cloneVector :: Storable t => Vector t -> IO (Vector t) | 214 | cloneVector :: Storable t => Vector t -> IO (Vector t) |
211 | cloneVector (v@V {dim=n}) = do | 215 | cloneVector (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" |
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs index ffc4f17..1f8b5a0 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs | |||
@@ -203,7 +203,7 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati | |||
203 | --------------------------------------------------------------- | 203 | --------------------------------------------------------------- |
204 | 204 | ||
205 | instance (Storable a, Num (Vector a)) => Monoid (Vector a) where | 205 | instance (Storable a, Num (Vector a)) => Monoid (Vector a) where |
206 | mempty = 0 { dim = 0 } | 206 | mempty = 0 { idim = 0 } |
207 | mappend a b = mconcat [a,b] | 207 | mappend a b = mconcat [a,b] |
208 | mconcat = j . filter ((>0).dim) | 208 | mconcat = j . filter ((>0).dim) |
209 | where j [] = mempty | 209 | where j [] = mempty |