diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-09-12 19:09:47 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-09-12 19:09:47 +0000 |
commit | 0ff13d993b880739295de343bca62f06fac5ca0c (patch) | |
tree | 252a51b4314c19c04a9eda87973eeaae63167a41 /lib/Data/Packed | |
parent | cd937c2be2900b8f13506d9ae7c731ad43d74e05 (diff) |
documentation
Diffstat (limited to 'lib/Data/Packed')
-rw-r--r-- | lib/Data/Packed/Internal.hs | 5 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Common.hs | 24 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 116 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 62 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/aux.c | 8 | ||||
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 62 | ||||
-rw-r--r-- | lib/Data/Packed/Plot.hs | 1 | ||||
-rw-r--r-- | lib/Data/Packed/Vector.hs | 35 |
8 files changed, 190 insertions, 123 deletions
diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs index 822bb1b..e5028d3 100644 --- a/lib/Data/Packed/Internal.hs +++ b/lib/Data/Packed/Internal.hs | |||
@@ -11,15 +11,14 @@ | |||
11 | -- Reexports all internal modules | 11 | -- Reexports all internal modules |
12 | -- | 12 | -- |
13 | ----------------------------------------------------------------------------- | 13 | ----------------------------------------------------------------------------- |
14 | -- #hide | ||
14 | 15 | ||
15 | module Data.Packed.Internal ( | 16 | module Data.Packed.Internal ( |
16 | module Data.Packed.Internal.Common, | 17 | module Data.Packed.Internal.Common, |
17 | module Data.Packed.Internal.Vector, | 18 | module Data.Packed.Internal.Vector, |
18 | module Data.Packed.Internal.Matrix, | 19 | module Data.Packed.Internal.Matrix, |
19 | -- module Data.Packed.Internal.Tensor | ||
20 | ) where | 20 | ) where |
21 | 21 | ||
22 | import Data.Packed.Internal.Common | 22 | import Data.Packed.Internal.Common |
23 | import Data.Packed.Internal.Vector | 23 | import Data.Packed.Internal.Vector |
24 | import Data.Packed.Internal.Matrix | 24 | import Data.Packed.Internal.Matrix \ No newline at end of file |
25 | --import Data.Packed.Internal.Tensor \ No newline at end of file | ||
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs index 5548285..218fb6b 100644 --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs | |||
@@ -9,9 +9,10 @@ | |||
9 | -- Stability : provisional | 9 | -- Stability : provisional |
10 | -- Portability : portable (uses FFI) | 10 | -- Portability : portable (uses FFI) |
11 | -- | 11 | -- |
12 | -- Common tools | 12 | -- Development utilities. |
13 | -- | 13 | -- |
14 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
15 | -- #hide | ||
15 | 16 | ||
16 | module Data.Packed.Internal.Common where | 17 | module Data.Packed.Internal.Common where |
17 | 18 | ||
@@ -33,11 +34,15 @@ instance (Storable a, RealFloat a) => Storable (Complex a) where -- | |||
33 | poke p (a :+ b) = pokeArray (castPtr p) [a,b] -- | 34 | poke p (a :+ b) = pokeArray (castPtr p) [a,b] -- |
34 | ---------------------------------------------------------------------- | 35 | ---------------------------------------------------------------------- |
35 | 36 | ||
37 | -- | @debug x = trace (show x) x@ | ||
38 | debug :: (Show a) => a -> a | ||
36 | debug x = trace (show x) x | 39 | debug x = trace (show x) x |
37 | 40 | ||
41 | -- | useful for expressions like @sortBy (compare \`on\` length)@ | ||
38 | on :: (a -> a -> b) -> (t -> a) -> t -> t -> b | 42 | on :: (a -> a -> b) -> (t -> a) -> t -> t -> b |
39 | on f g = \x y -> f (g x) (g y) | 43 | on f g = \x y -> f (g x) (g y) |
40 | 44 | ||
45 | -- | @partit 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@ | ||
41 | partit :: Int -> [a] -> [[a]] | 46 | partit :: Int -> [a] -> [[a]] |
42 | partit _ [] = [] | 47 | partit _ [] = [] |
43 | partit n l = take n l : partit n (drop n l) | 48 | partit n l = take n l : partit n (drop n l) |
@@ -50,19 +55,20 @@ common f = commonval . map f where | |||
50 | commonval [a] = Just a | 55 | commonval [a] = Just a |
51 | commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing | 56 | commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing |
52 | 57 | ||
58 | -- | postfix function application (@flip ($)@) | ||
53 | (//) :: x -> (x -> y) -> y | 59 | (//) :: x -> (x -> y) -> y |
54 | infixl 0 // | 60 | infixl 0 // |
55 | (//) = flip ($) | 61 | (//) = flip ($) |
56 | 62 | ||
57 | -- our codes should start from 1024 | 63 | -- GSL error codes are <= 1024 |
58 | 64 | -- | error codes for the auxiliary functions required by the wrappers | |
59 | errorCode :: Int -> String | 65 | errorCode :: Int -> String |
60 | errorCode 1000 = "bad size" | 66 | errorCode 2000 = "bad size" |
61 | errorCode 1001 = "bad function code" | 67 | errorCode 2001 = "bad function code" |
62 | errorCode 1002 = "memory problem" | 68 | errorCode 2002 = "memory problem" |
63 | errorCode 1003 = "bad file" | 69 | errorCode 2003 = "bad file" |
64 | errorCode 1004 = "singular" | 70 | errorCode 2004 = "singular" |
65 | errorCode 1005 = "didn't converge" | 71 | errorCode 2005 = "didn't converge" |
66 | errorCode n = "code "++show n | 72 | errorCode n = "code "++show n |
67 | 73 | ||
68 | {- | conversion of Haskell functions into function pointers that can be used in the C side | 74 | {- | conversion of Haskell functions into function pointers that can be used in the C side |
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 2db4838..f9dd9a9 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -12,7 +12,7 @@ | |||
12 | -- Internal matrix representation | 12 | -- Internal matrix representation |
13 | -- | 13 | -- |
14 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
15 | -- --#hide | 15 | -- #hide |
16 | 16 | ||
17 | module Data.Packed.Internal.Matrix where | 17 | module Data.Packed.Internal.Matrix where |
18 | 18 | ||
@@ -57,10 +57,14 @@ import Data.Maybe(fromJust) | |||
57 | 57 | ||
58 | data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) | 58 | data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) |
59 | 59 | ||
60 | -- | Matrix representation suitable for GSL and LAPACK computations. | ||
60 | data Matrix t = MC { rows :: Int, cols :: Int, cdat :: Vector t, fdat :: Vector t } | 61 | data Matrix t = MC { rows :: Int, cols :: Int, cdat :: Vector t, fdat :: Vector t } |
61 | | MF { rows :: Int, cols :: Int, fdat :: Vector t, cdat :: Vector t } | 62 | | MF { rows :: Int, cols :: Int, fdat :: Vector t, cdat :: Vector t } |
62 | 63 | ||
63 | -- transposition just changes the data order | 64 | -- MC: preferred by C, fdat may require a transposition |
65 | -- MF: preferred by LAPACK, cdat may require a transposition | ||
66 | |||
67 | -- | matrix transpose | ||
64 | trans :: Matrix t -> Matrix t | 68 | trans :: Matrix t -> Matrix t |
65 | trans MC {rows = r, cols = c, cdat = d, fdat = dt } = MF {rows = c, cols = r, fdat = d, cdat = dt } | 69 | trans MC {rows = r, cols = c, cdat = d, fdat = dt } = MF {rows = c, cols = r, fdat = d, cdat = dt } |
66 | trans MF {rows = r, cols = c, fdat = d, cdat = dt } = MC {rows = c, cols = r, cdat = d, fdat = dt } | 70 | trans MF {rows = r, cols = c, fdat = d, cdat = dt } = MC {rows = c, cols = r, cdat = d, fdat = dt } |
@@ -166,32 +170,29 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 | |||
166 | 170 | ||
167 | ---------------------------------------------------------------- | 171 | ---------------------------------------------------------------- |
168 | 172 | ||
169 | -- | element types for which optimized matrix computations are provided | 173 | -- | Optimized matrix computations are provided for elements in the Field class. |
170 | class Storable a => Field a where | 174 | class Storable a => Field a where |
171 | -- | @constant val n@ creates a vector with @n@ elements, all equal to @val@. | 175 | constantD :: a -> Int -> Vector a |
172 | constant :: a -> Int -> Vector a | ||
173 | transdata :: Int -> Vector a -> Int -> Vector a | 176 | transdata :: Int -> Vector a -> Int -> Vector a |
174 | multiplyD :: MatrixOrder -> Matrix a -> Matrix a -> Matrix a | 177 | multiplyD :: Matrix a -> Matrix a -> Matrix a |
175 | -- | extracts a submatrix froma a matrix | 178 | subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position |
176 | subMatrix :: (Int,Int) -- ^ (r0,c0) starting position | 179 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix |
177 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | 180 | -> Matrix a -> Matrix a |
178 | -> Matrix a -> Matrix a | 181 | diagD :: Vector a -> Matrix a |
179 | -- | creates a square matrix with the given diagonal | ||
180 | diag :: Vector a -> Matrix a | ||
181 | 182 | ||
182 | instance Field Double where | 183 | instance Field Double where |
183 | constant = constantR | 184 | constantD = constantR |
184 | transdata = transdataR | 185 | transdata = transdataR |
185 | multiplyD = multiplyR | 186 | multiplyD = multiplyR |
186 | subMatrix = subMatrixR | 187 | subMatrixD = subMatrixR |
187 | diag = diagR | 188 | diagD = diagR |
188 | 189 | ||
189 | instance Field (Complex Double) where | 190 | instance Field (Complex Double) where |
190 | constant = constantC | 191 | constantD = constantC |
191 | transdata = transdataC | 192 | transdata = transdataC |
192 | multiplyD = multiplyC | 193 | multiplyD = multiplyC |
193 | subMatrix = subMatrixC | 194 | subMatrixD = subMatrixC |
194 | diag = diagC | 195 | diagD = diagC |
195 | 196 | ||
196 | ------------------------------------------------------------------ | 197 | ------------------------------------------------------------------ |
197 | 198 | ||
@@ -209,6 +210,15 @@ dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unw | |||
209 | 210 | ||
210 | ------------------------------------------------------------------ | 211 | ------------------------------------------------------------------ |
211 | 212 | ||
213 | (>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a | ||
214 | r >|< c = f where | ||
215 | f l | dim v == r*c = matrixFromVector ColumnMajor c v | ||
216 | | otherwise = error $ "inconsistent list size = " | ||
217 | ++show (dim v) ++" in ("++show r++"><"++show c++")" | ||
218 | where v = fromList l | ||
219 | |||
220 | ------------------------------------------------------------------- | ||
221 | |||
212 | transdataR :: Int -> Vector Double -> Int -> Vector Double | 222 | transdataR :: Int -> Vector Double -> Int -> Vector Double |
213 | transdataR = transdataAux ctransR | 223 | transdataR = transdataAux ctransR |
214 | 224 | ||
@@ -237,10 +247,10 @@ foreign import ccall safe "aux.h transC" | |||
237 | gmatC MF {rows = r, cols = c, fdat = d} f = f 1 c r (ptr d) | 247 | gmatC MF {rows = r, cols = c, fdat = d} f = f 1 c r (ptr d) |
238 | gmatC MC {rows = r, cols = c, cdat = d} f = f 0 r c (ptr d) | 248 | gmatC MC {rows = r, cols = c, cdat = d} f = f 0 r c (ptr d) |
239 | 249 | ||
240 | multiplyAux fun order a b = unsafePerformIO $ do | 250 | multiplyAux fun a b = unsafePerformIO $ do |
241 | when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++ | 251 | when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++ |
242 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) | 252 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) |
243 | r <- createMatrix order (rows a) (cols b) | 253 | r <- createMatrix RowMajor (rows a) (cols b) |
244 | fun // gmatC a // gmatC b // mat dat r // check "multiplyAux" [dat a, dat b] | 254 | fun // gmatC a // gmatC b // mat dat r // check "multiplyAux" [dat a, dat b] |
245 | return r | 255 | return r |
246 | 256 | ||
@@ -258,33 +268,41 @@ foreign import ccall safe "aux.h multiplyC" | |||
258 | -> Int -> Int -> Ptr (Complex Double) | 268 | -> Int -> Int -> Ptr (Complex Double) |
259 | -> IO Int | 269 | -> IO Int |
260 | 270 | ||
261 | multiply :: (Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a | 271 | multiply' :: (Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a |
262 | multiply RowMajor a b = multiplyD RowMajor a b | 272 | multiply' RowMajor a b = multiplyD a b |
263 | multiply ColumnMajor a b = MF {rows = c, cols = r, fdat = d, cdat = dt } | 273 | multiply' ColumnMajor a b = trans $ multiplyD (trans b) (trans a) |
264 | where MC {rows = r, cols = c, cdat = d, fdat = dt } = multiplyD RowMajor (trans b) (trans a) | 274 | |
265 | -- FIXME using MatrixFromVector | 275 | |
276 | -- | matrix product | ||
277 | multiply :: (Field a) => Matrix a -> Matrix a -> Matrix a | ||
278 | multiply = multiplyD | ||
266 | 279 | ||
267 | ---------------------------------------------------------------------- | 280 | ---------------------------------------------------------------------- |
268 | 281 | ||
269 | -- | extraction of a submatrix of a real matrix | 282 | -- | extraction of a submatrix from a real matrix |
270 | subMatrixR :: (Int,Int) -- ^ (r0,c0) starting position | 283 | subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double |
271 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | ||
272 | -> Matrix Double -> Matrix Double | ||
273 | subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do | 284 | subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do |
274 | r <- createMatrix RowMajor rt ct | 285 | r <- createMatrix RowMajor rt ct |
275 | c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r] | 286 | c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r] |
276 | return r | 287 | return r |
277 | foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM | 288 | foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM |
278 | 289 | ||
279 | -- | extraction of a submatrix of a complex matrix | 290 | -- | extraction of a submatrix from a complex matrix |
280 | subMatrixC :: (Int,Int) -- ^ (r0,c0) starting position | 291 | subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double) |
281 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | ||
282 | -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
283 | subMatrixC (r0,c0) (rt,ct) x = | 292 | subMatrixC (r0,c0) (rt,ct) x = |
284 | reshape ct . asComplex . cdat . | 293 | reshape ct . asComplex . cdat . |
285 | subMatrixR (r0,2*c0) (rt,2*ct) . | 294 | subMatrixR (r0,2*c0) (rt,2*ct) . |
286 | reshape (2*cols x) . asReal . cdat $ x | 295 | reshape (2*cols x) . asReal . cdat $ x |
287 | 296 | ||
297 | -- | Extracts a submatrix from a matrix. | ||
298 | subMatrix :: Field a | ||
299 | => (Int,Int) -- ^ (r0,c0) starting position | ||
300 | -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix | ||
301 | -> Matrix a -- ^ input matrix | ||
302 | -> Matrix a -- ^ result | ||
303 | subMatrix = subMatrixD | ||
304 | |||
305 | |||
288 | --------------------------------------------------------------------- | 306 | --------------------------------------------------------------------- |
289 | 307 | ||
290 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do | 308 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do |
@@ -302,6 +320,10 @@ diagC :: Vector (Complex Double) -> Matrix (Complex Double) | |||
302 | diagC = diagAux c_diagC "diagC" | 320 | diagC = diagAux c_diagC "diagC" |
303 | foreign import ccall "aux.h diagC" c_diagC :: TCVCM | 321 | foreign import ccall "aux.h diagC" c_diagC :: TCVCM |
304 | 322 | ||
323 | -- | creates a square matrix with the given diagonal | ||
324 | diag :: Field a => Vector a -> Matrix a | ||
325 | diag = diagD | ||
326 | |||
305 | ------------------------------------------------------------------------ | 327 | ------------------------------------------------------------------------ |
306 | 328 | ||
307 | constantAux fun x n = unsafePerformIO $ do | 329 | constantAux fun x n = unsafePerformIO $ do |
@@ -321,6 +343,28 @@ constantC = constantAux cconstantC | |||
321 | foreign import ccall safe "aux.h constantC" | 343 | foreign import ccall safe "aux.h constantC" |
322 | cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int | 344 | cconstantC :: Ptr (Complex Double) -> TCV -- Complex Double :> IO Int |
323 | 345 | ||
346 | {- | creates a vector with a given number of equal components: | ||
347 | |||
348 | @> constant 2 7 | ||
349 | 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ | ||
350 | -} | ||
351 | constant :: Field a => a -> Int -> Vector a | ||
352 | constant = constantD | ||
353 | |||
354 | -------------------------------------------------------------------------- | ||
355 | |||
356 | -- | obtains the complex conjugate of a complex vector | ||
357 | conj :: Vector (Complex Double) -> Vector (Complex Double) | ||
358 | conj v = asComplex $ cdat $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1]) | ||
359 | |||
360 | -- | creates a complex vector from vectors with real and imaginary parts | ||
361 | toComplex :: (Vector Double, Vector Double) -> Vector (Complex Double) | ||
362 | toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i] | ||
363 | |||
364 | -- | converts a real vector into a complex representation (with zero imaginary parts) | ||
365 | comp :: Vector Double -> Vector (Complex Double) | ||
366 | comp v = toComplex (v,constant 0 (dim v)) | ||
367 | |||
324 | ------------------------------------------------------------------------- | 368 | ------------------------------------------------------------------------- |
325 | 369 | ||
326 | -- Generic definitions | 370 | -- Generic definitions |
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs index 0d9dc70..f0ef8b6 100644 --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs | |||
@@ -12,6 +12,7 @@ | |||
12 | -- Vector implementation | 12 | -- Vector implementation |
13 | -- | 13 | -- |
14 | ----------------------------------------------------------------------------- | 14 | ----------------------------------------------------------------------------- |
15 | -- #hide | ||
15 | 16 | ||
16 | module Data.Packed.Internal.Vector where | 17 | module Data.Packed.Internal.Vector where |
17 | 18 | ||
@@ -24,34 +25,39 @@ import Debug.Trace(trace) | |||
24 | import Foreign.C.String(peekCString) | 25 | import Foreign.C.String(peekCString) |
25 | import Foreign.C.Types | 26 | import Foreign.C.Types |
26 | 27 | ||
27 | 28 | -- | A one-dimensional array of objects stored in a contiguous memory block. | |
28 | data Vector t = V { dim :: Int | 29 | data Vector t = V { dim :: Int -- ^ number of elements |
29 | , fptr :: ForeignPtr t | 30 | , fptr :: ForeignPtr t -- ^ foreign pointer to the memory block |
30 | , ptr :: Ptr t | 31 | , ptr :: Ptr t -- ^ ordinary pointer to the actual starting address (usually the same) |
31 | } | 32 | } |
32 | 33 | ||
34 | -- | check the error code and touch foreign ptr of vector arguments (if any) | ||
33 | check :: String -> [Vector a] -> IO Int -> IO () | 35 | check :: String -> [Vector a] -> IO Int -> IO () |
34 | check msg ls f = do | 36 | check msg ls f = do |
35 | err <- f | 37 | err <- f |
36 | when (err/=0) $ if err > 999 -- FIXME, it should be 1024 | 38 | when (err/=0) $ if err > 1024 |
37 | then (error (msg++": "++errorCode err)) | 39 | then (error (msg++": "++errorCode err)) -- our errors |
38 | else do | 40 | else do -- GSL errors |
39 | ps <- gsl_strerror err | 41 | ps <- gsl_strerror err |
40 | s <- peekCString ps | 42 | s <- peekCString ps |
41 | error (msg++": "++s) | 43 | error (msg++": "++s) |
42 | mapM_ (touchForeignPtr . fptr) ls | 44 | mapM_ (touchForeignPtr . fptr) ls |
43 | return () | 45 | return () |
44 | 46 | ||
47 | -- | description of GSL error codes | ||
45 | foreign import ccall "aux.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar) | 48 | foreign import ccall "aux.h gsl_strerror" gsl_strerror :: Int -> IO (Ptr CChar) |
46 | 49 | ||
50 | -- | signature of foreign functions admitting C-style vectors | ||
47 | type Vc t s = Int -> Ptr t -> s | 51 | type Vc t s = Int -> Ptr t -> s |
48 | -- not yet admitted by my haddock version | 52 | -- not yet admitted by my haddock version |
49 | -- infixr 5 :> | 53 | -- infixr 5 :> |
50 | -- type t :> s = Vc t s | 54 | -- type t :> s = Vc t s |
51 | 55 | ||
56 | -- | adaptation of our vectors to be admitted by foreign functions: @f \/\/ vec v@ | ||
52 | vec :: Vector t -> (Vc t s) -> s | 57 | vec :: Vector t -> (Vc t s) -> s |
53 | vec v f = f (dim v) (ptr v) | 58 | vec v f = f (dim v) (ptr v) |
54 | 59 | ||
60 | -- | allocates memory for a new vector | ||
55 | createVector :: Storable a => Int -> IO (Vector a) | 61 | createVector :: Storable a => Int -> IO (Vector a) |
56 | createVector n = do | 62 | createVector n = do |
57 | when (n <= 0) $ error ("trying to createVector of dim "++show n) | 63 | when (n <= 0) $ error ("trying to createVector of dim "++show n) |
@@ -60,6 +66,12 @@ createVector n = do | |||
60 | --putStrLn ("\n---------> V"++show n) | 66 | --putStrLn ("\n---------> V"++show n) |
61 | return $ V n fp p | 67 | return $ V n fp p |
62 | 68 | ||
69 | {- | creates a Vector from a list: | ||
70 | |||
71 | @> fromList [2,3,5,7] | ||
72 | 4 |> [2.0,3.0,5.0,7.0]@ | ||
73 | |||
74 | -} | ||
63 | fromList :: Storable a => [a] -> Vector a | 75 | fromList :: Storable a => [a] -> Vector a |
64 | fromList l = unsafePerformIO $ do | 76 | fromList l = unsafePerformIO $ do |
65 | v <- createVector (length l) | 77 | v <- createVector (length l) |
@@ -67,14 +79,25 @@ fromList l = unsafePerformIO $ do | |||
67 | f // vec v // check "fromList" [] | 79 | f // vec v // check "fromList" [] |
68 | return v | 80 | return v |
69 | 81 | ||
82 | {- | extracts the Vector elements to a list | ||
83 | |||
84 | @> toList (linspace 5 (1,10)) | ||
85 | [1.0,3.25,5.5,7.75,10.0]@ | ||
86 | |||
87 | -} | ||
70 | toList :: Storable a => Vector a -> [a] | 88 | toList :: Storable a => Vector a -> [a] |
71 | toList v = unsafePerformIO $ peekArray (dim v) (ptr v) | 89 | toList v = unsafePerformIO $ peekArray (dim v) (ptr v) |
72 | 90 | ||
91 | -- | an alternative to 'fromList' with explicit dimension, used also in the instances for Show (Vector a). | ||
92 | (|>) :: (Storable a) => Int -> [a] -> Vector a | ||
93 | infixl 9 |> | ||
73 | n |> l = if length l == n then fromList l else error "|> with wrong size" | 94 | n |> l = if length l == n then fromList l else error "|> with wrong size" |
74 | 95 | ||
96 | -- | access to Vector elements without range checking | ||
75 | at' :: Storable a => Vector a -> Int -> a | 97 | at' :: Storable a => Vector a -> Int -> a |
76 | at' v n = unsafePerformIO $ peekElemOff (ptr v) n | 98 | at' v n = unsafePerformIO $ peekElemOff (ptr v) n |
77 | 99 | ||
100 | -- | access to Vector elements with range checking. | ||
78 | at :: Storable a => Vector a -> Int -> a | 101 | at :: Storable a => Vector a -> Int -> a |
79 | at v n | n >= 0 && n < dim v = at' v n | 102 | at v n | n >= 0 && n < dim v = at' v n |
80 | | otherwise = error "vector index out of range" | 103 | | otherwise = error "vector index out of range" |
@@ -82,9 +105,14 @@ at v n | n >= 0 && n < dim v = at' v n | |||
82 | instance (Show a, Storable a) => (Show (Vector a)) where | 105 | instance (Show a, Storable a) => (Show (Vector a)) where |
83 | show v = (show (dim v))++" |> " ++ show (toList v) | 106 | show v = (show (dim v))++" |> " ++ show (toList v) |
84 | 107 | ||
85 | -- | creates a Vector taking a number of consecutive toList from another Vector | 108 | {- | takes a number of consecutive elements from a Vector |
109 | |||
110 | @> subVector 2 3 (fromList [1..10]) | ||
111 | 3 |> [3.0,4.0,5.0]@ | ||
112 | |||
113 | -} | ||
86 | subVector :: Storable t => Int -- ^ index of the starting element | 114 | subVector :: Storable t => Int -- ^ index of the starting element |
87 | -> Int -- ^ number of toList to extract | 115 | -> Int -- ^ number of elements to extract |
88 | -> Vector t -- ^ source | 116 | -> Vector t -- ^ source |
89 | -> Vector t -- ^ result | 117 | -> Vector t -- ^ result |
90 | subVector k l (v@V {dim=n, ptr=p, fptr=fp}) | 118 | subVector k l (v@V {dim=n, ptr=p, fptr=fp}) |
@@ -100,13 +128,23 @@ subVector' k l (v@V {dim=n, ptr=p, fptr=fp}) | |||
100 | | otherwise = v {dim=l, ptr=advancePtr p k} | 128 | | otherwise = v {dim=l, ptr=advancePtr p k} |
101 | 129 | ||
102 | 130 | ||
103 | -- | Reads a vector position. | 131 | {- | Reads a vector position: |
132 | |||
133 | @> fromList [0..9] \@\> 7 | ||
134 | 7.0@ | ||
135 | |||
136 | -} | ||
104 | (@>) :: Storable t => Vector t -> Int -> t | 137 | (@>) :: Storable t => Vector t -> Int -> t |
105 | infixl 9 @> | 138 | infixl 9 @> |
106 | (@>) = at | 139 | (@>) = at |
107 | 140 | ||
108 | 141 | ||
109 | -- | creates a new Vector by joining a list of Vectors | 142 | {- | creates a new Vector by joining a list of Vectors |
143 | |||
144 | @> join [fromList [1..5], constant 1 3] | ||
145 | 8 |> [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]@ | ||
146 | |||
147 | -} | ||
110 | join :: Storable t => [Vector t] -> Vector t | 148 | join :: Storable t => [Vector t] -> Vector t |
111 | join [] = error "joining zero vectors" | 149 | join [] = error "joining zero vectors" |
112 | join as = unsafePerformIO $ do | 150 | join as = unsafePerformIO $ do |
@@ -131,9 +169,11 @@ asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v), ptr = ca | |||
131 | 169 | ||
132 | ---------------------------------------------------------------- | 170 | ---------------------------------------------------------------- |
133 | 171 | ||
172 | -- | map on Vectors | ||
134 | liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b | 173 | liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b |
135 | liftVector f = fromList . map f . toList | 174 | liftVector f = fromList . map f . toList |
136 | 175 | ||
176 | -- | zipWith for Vectors | ||
137 | liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c | 177 | liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c |
138 | liftVector2 f u v = fromList $ zipWith f (toList u) (toList v) | 178 | liftVector2 f u v = fromList $ zipWith f (toList u) (toList v) |
139 | 179 | ||
diff --git a/lib/Data/Packed/Internal/aux.c b/lib/Data/Packed/Internal/aux.c index fe611e2..3db3535 100644 --- a/lib/Data/Packed/Internal/aux.c +++ b/lib/Data/Packed/Internal/aux.c | |||
@@ -58,10 +58,10 @@ | |||
58 | #define GCVEC(A) int A##n, gsl_complex*A##p | 58 | #define GCVEC(A) int A##n, gsl_complex*A##p |
59 | #define KGCVEC(A) int A##n, const gsl_complex*A##p | 59 | #define KGCVEC(A) int A##n, const gsl_complex*A##p |
60 | 60 | ||
61 | #define BAD_SIZE 1000 | 61 | #define BAD_SIZE 2000 |
62 | #define BAD_CODE 1001 | 62 | #define BAD_CODE 2001 |
63 | #define MEM 1002 | 63 | #define MEM 2002 |
64 | #define BAD_FILE 1003 | 64 | #define BAD_FILE 2003 |
65 | 65 | ||
66 | int transR(KRMAT(x),RMAT(t)) { | 66 | int transR(KRMAT(x),RMAT(t)) { |
67 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | 67 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); |
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index 45aaaba..01e8133 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs | |||
@@ -13,19 +13,19 @@ | |||
13 | ----------------------------------------------------------------------------- | 13 | ----------------------------------------------------------------------------- |
14 | 14 | ||
15 | module Data.Packed.Matrix ( | 15 | module Data.Packed.Matrix ( |
16 | Matrix(rows,cols), | 16 | Field, |
17 | fromLists, toLists, (><), (>|<), (@@>), | 17 | Matrix,rows,cols, |
18 | (><), reshape, flatten, | ||
19 | fromLists, toLists, | ||
20 | (@@>), | ||
18 | trans, conjTrans, | 21 | trans, conjTrans, |
19 | reshape, flatten, asRow, asColumn, | 22 | asRow, asColumn, |
20 | fromRows, toRows, fromColumns, toColumns, fromBlocks, | 23 | fromRows, toRows, fromColumns, toColumns, |
21 | joinVert, joinHoriz, | 24 | fromBlocks, joinVert, joinHoriz, |
22 | flipud, fliprl, | 25 | flipud, fliprl, |
26 | subMatrix, takeRows, dropRows, takeColumns, dropColumns, | ||
27 | diag, takeDiag, diagRect, ident, | ||
23 | liftMatrix, liftMatrix2, | 28 | liftMatrix, liftMatrix2, |
24 | multiply, | ||
25 | outer, | ||
26 | subMatrix, | ||
27 | takeRows, dropRows, takeColumns, dropColumns, | ||
28 | diag, takeDiag, diagRect, ident | ||
29 | ) where | 29 | ) where |
30 | 30 | ||
31 | import Data.Packed.Internal | 31 | import Data.Packed.Internal |
@@ -83,6 +83,18 @@ takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols | |||
83 | ident :: (Num t, Field t) => Int -> Matrix t | 83 | ident :: (Num t, Field t) => Int -> Matrix t |
84 | ident n = diag (constant 1 n) | 84 | ident n = diag (constant 1 n) |
85 | 85 | ||
86 | ------------------------------------------------------------ | ||
87 | |||
88 | {- | An easy way to create a matrix: | ||
89 | |||
90 | @\> (2><3)[1..6] | ||
91 | (2><3) | ||
92 | [ 1.0, 2.0, 3.0 | ||
93 | , 4.0, 5.0, 6.0 ]@ | ||
94 | |||
95 | This is the format produced by the instances of Show (Matrix a), which | ||
96 | can also be used for input. | ||
97 | -} | ||
86 | (><) :: (Field a) => Int -> Int -> [a] -> Matrix a | 98 | (><) :: (Field a) => Int -> Int -> [a] -> Matrix a |
87 | r >< c = f where | 99 | r >< c = f where |
88 | f l | dim v == r*c = matrixFromVector RowMajor c v | 100 | f l | dim v == r*c = matrixFromVector RowMajor c v |
@@ -90,13 +102,6 @@ r >< c = f where | |||
90 | ++show (dim v) ++" in ("++show r++"><"++show c++")" | 102 | ++show (dim v) ++" in ("++show r++"><"++show c++")" |
91 | where v = fromList l | 103 | where v = fromList l |
92 | 104 | ||
93 | (>|<) :: (Field a) => Int -> Int -> [a] -> Matrix a | ||
94 | r >|< c = f where | ||
95 | f l | dim v == r*c = matrixFromVector ColumnMajor c v | ||
96 | | otherwise = error $ "inconsistent list size = " | ||
97 | ++show (dim v) ++" in ("++show r++"><"++show c++")" | ||
98 | where v = fromList l | ||
99 | |||
100 | ---------------------------------------------------------------- | 105 | ---------------------------------------------------------------- |
101 | 106 | ||
102 | -- | Creates a matrix with the first n rows of another matrix | 107 | -- | Creates a matrix with the first n rows of another matrix |
@@ -117,12 +122,19 @@ dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat | |||
117 | {- | Creates a vector by concatenation of rows | 122 | {- | Creates a vector by concatenation of rows |
118 | 123 | ||
119 | @\> flatten ('ident' 3) | 124 | @\> flatten ('ident' 3) |
120 | 9 # [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ | 125 | 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ |
121 | -} | 126 | -} |
122 | flatten :: Field t => Matrix t -> Vector t | 127 | flatten :: Field t => Matrix t -> Vector t |
123 | flatten = cdat | 128 | flatten = cdat |
124 | 129 | ||
125 | -- | Creates a 'Matrix' from a list of lists (considered as rows). | 130 | {- | Creates a 'Matrix' from a list of lists (considered as rows). |
131 | |||
132 | @\> fromLists [[1,2],[3,4],[5,6]] | ||
133 | (3><2) | ||
134 | [ 1.0, 2.0 | ||
135 | , 3.0, 4.0 | ||
136 | , 5.0, 6.0 ]@ | ||
137 | -} | ||
126 | fromLists :: Field t => [[t]] -> Matrix t | 138 | fromLists :: Field t => [[t]] -> Matrix t |
127 | fromLists = fromRows . map fromList | 139 | fromLists = fromRows . map fromList |
128 | 140 | ||
@@ -137,15 +149,3 @@ asColumn v = reshape 1 v | |||
137 | 149 | ||
138 | ------------------------------------------------ | 150 | ------------------------------------------------ |
139 | 151 | ||
140 | {- | Outer product of two vectors. | ||
141 | |||
142 | @\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] | ||
143 | (3><3) | ||
144 | [ 5.0, 2.0, 3.0 | ||
145 | , 10.0, 4.0, 6.0 | ||
146 | , 15.0, 6.0, 9.0 ]@ | ||
147 | -} | ||
148 | outer :: (Num t, Field t) => Vector t -> Vector t -> Matrix t | ||
149 | outer u v = multiply RowMajor r c | ||
150 | where r = matrixFromVector RowMajor 1 u | ||
151 | c = matrixFromVector RowMajor (dim v) v | ||
diff --git a/lib/Data/Packed/Plot.hs b/lib/Data/Packed/Plot.hs index a0a4aae..7d63426 100644 --- a/lib/Data/Packed/Plot.hs +++ b/lib/Data/Packed/Plot.hs | |||
@@ -26,6 +26,7 @@ module Data.Packed.Plot( | |||
26 | 26 | ||
27 | import Data.Packed.Vector | 27 | import Data.Packed.Vector |
28 | import Data.Packed.Matrix | 28 | import Data.Packed.Matrix |
29 | import LinearAlgebra.Linear(outer) | ||
29 | import GSL.Vector(FunCodeS(Max,Min),toScalarR) | 30 | import GSL.Vector(FunCodeS(Max,Min),toScalarR) |
30 | import Data.List(intersperse) | 31 | import Data.List(intersperse) |
31 | import System | 32 | import System |
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs index 867b77b..9818d58 100644 --- a/lib/Data/Packed/Vector.hs +++ b/lib/Data/Packed/Vector.hs | |||
@@ -13,33 +13,17 @@ | |||
13 | ----------------------------------------------------------------------------- | 13 | ----------------------------------------------------------------------------- |
14 | 14 | ||
15 | module Data.Packed.Vector ( | 15 | module Data.Packed.Vector ( |
16 | Vector(dim), Field, | 16 | Vector, |
17 | fromList, toList, | 17 | fromList, (|>), toList, |
18 | (@>), | 18 | dim, (@>), |
19 | subVector, join, | 19 | subVector, join, |
20 | constant, | 20 | constant, linspace, |
21 | toComplex, comp, | 21 | |
22 | conj, | ||
23 | dot, | ||
24 | linspace, | ||
25 | liftVector, liftVector2 | 22 | liftVector, liftVector2 |
26 | ) where | 23 | ) where |
27 | 24 | ||
28 | import Data.Packed.Internal | 25 | import Data.Packed.Internal |
29 | import Complex | 26 | import Complex |
30 | --import GSL.Vector | ||
31 | |||
32 | -- | creates a complex vector from vectors with real and imaginary parts | ||
33 | toComplex :: (Vector Double, Vector Double) -> Vector (Complex Double) | ||
34 | toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i] | ||
35 | |||
36 | -- | obtains the complex conjugate of a complex vector | ||
37 | conj :: Vector (Complex Double) -> Vector (Complex Double) | ||
38 | conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1]) | ||
39 | where mulC = multiply RowMajor | ||
40 | |||
41 | comp :: Vector Double -> Vector (Complex Double) | ||
42 | comp v = toComplex (v,constant 0 (dim v)) | ||
43 | 27 | ||
44 | {- | Creates a real vector containing a range of values: | 28 | {- | Creates a real vector containing a range of values: |
45 | 29 | ||
@@ -47,12 +31,5 @@ comp v = toComplex (v,constant 0 (dim v)) | |||
47 | 5 |> [-3.0,-0.5,2.0,4.5,7.0]@ | 31 | 5 |> [-3.0,-0.5,2.0,4.5,7.0]@ |
48 | -} | 32 | -} |
49 | linspace :: Int -> (Double, Double) -> Vector Double | 33 | linspace :: Int -> (Double, Double) -> Vector Double |
50 | linspace n (a,b) = fromList [a::Double,a+delta .. b] | 34 | linspace n (a,b) = fromList [a, a+delta .. b] |
51 | where delta = (b-a)/(fromIntegral n -1) | 35 | where delta = (b-a)/(fromIntegral n -1) |
52 | |||
53 | |||
54 | dot :: (Field t) => Vector t -> Vector t -> t | ||
55 | dot u v = dat (multiply RowMajor r c) `at` 0 | ||
56 | where r = matrixFromVector RowMajor (dim u) u | ||
57 | c = matrixFromVector RowMajor 1 v | ||
58 | |||