diff options
author | Alberto Ruiz <aruiz@um.es> | 2009-05-13 09:39:25 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2009-05-13 09:39:25 +0000 |
commit | 6809103cf34a9345f8cb60a0ec3a8f55dd18d5ef (patch) | |
tree | c100100bfd26021c6130be16013e493db78cea16 | |
parent | edabccf9d9564573359ebf795f4ce9cce7b15816 (diff) |
added Development module
-rw-r--r-- | examples/devel/functions.c | 34 | ||||
-rw-r--r-- | examples/devel/wrappers.hs | 45 | ||||
-rw-r--r-- | hmatrix.cabal | 26 | ||||
-rw-r--r-- | lib/Data/Packed/Development.hs | 27 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Common.hs | 38 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 9 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 5 |
7 files changed, 178 insertions, 6 deletions
diff --git a/examples/devel/functions.c b/examples/devel/functions.c new file mode 100644 index 0000000..6885a0e --- /dev/null +++ b/examples/devel/functions.c | |||
@@ -0,0 +1,34 @@ | |||
1 | typedef struct { double r, i; } doublecomplex; | ||
2 | |||
3 | #define DVEC(A) int A##n, double*A##p | ||
4 | #define CVEC(A) int A##n, doublecomplex*A##p | ||
5 | #define DMAT(A) int A##r, int A##c, double*A##p | ||
6 | #define CMAT(A) int A##r, int A##c, doublecomplex*A##p | ||
7 | |||
8 | #define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) | ||
9 | |||
10 | /*-----------------------------------------------------*/ | ||
11 | |||
12 | int c_scale_vector(double s, DVEC(x), DVEC(y)) { | ||
13 | int k; | ||
14 | for (k=0; k<=yn; k++) { | ||
15 | yp[k] = s*xp[k]; | ||
16 | } | ||
17 | return 0; | ||
18 | } | ||
19 | |||
20 | /*-----------------------------------------------------*/ | ||
21 | |||
22 | int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) { | ||
23 | int i,j,sr,sc; | ||
24 | if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} | ||
25 | for (j=0; j<yn; j++) { | ||
26 | yp[j] = AT(m,j,j); | ||
27 | } | ||
28 | for (i=0; i<mr; i++) { | ||
29 | for (j=0; j<mc; j++) { | ||
30 | AT(z,i,j) = i==j?yp[i]:0; | ||
31 | } | ||
32 | } | ||
33 | return 0; | ||
34 | } | ||
diff --git a/examples/devel/wrappers.hs b/examples/devel/wrappers.hs new file mode 100644 index 0000000..f9e258a --- /dev/null +++ b/examples/devel/wrappers.hs | |||
@@ -0,0 +1,45 @@ | |||
1 | {-# LANGUAGE ForeignFunctionInterface #-} | ||
2 | |||
3 | -- $ ghc -O2 --make wrappers.hs functions.c | ||
4 | |||
5 | import Numeric.LinearAlgebra | ||
6 | import Data.Packed.Development | ||
7 | import Foreign(Ptr,unsafePerformIO) | ||
8 | import Foreign.C.Types(CInt) | ||
9 | |||
10 | ----------------------------------------------------- | ||
11 | |||
12 | main = do | ||
13 | print $ myScale 3.0 (fromList [1..10]) | ||
14 | print $ myDiag $ (3><5) [1..] | ||
15 | |||
16 | ----------------------------------------------------- | ||
17 | |||
18 | foreign import ccall "c_scale_vector" | ||
19 | cScaleVector :: Double -- scale | ||
20 | -> CInt -> Ptr Double -- argument | ||
21 | -> CInt -> Ptr Double -- result | ||
22 | -> IO CInt -- exit code | ||
23 | |||
24 | myScale s x = unsafePerformIO $ do | ||
25 | y <- createVector (dim x) | ||
26 | app2 (cScaleVector s) vec x vec y "cScaleVector" | ||
27 | return y | ||
28 | |||
29 | ----------------------------------------------------- | ||
30 | |||
31 | foreign import ccall "c_diag" | ||
32 | cDiag :: CInt -- matrix order | ||
33 | -> CInt -> CInt -> Ptr Double -- argument | ||
34 | -> CInt -> Ptr Double -- result1 | ||
35 | -> CInt -> CInt -> Ptr Double -- result2 | ||
36 | -> IO CInt -- exit code | ||
37 | |||
38 | myDiag m = unsafePerformIO $ do | ||
39 | y <- createVector (min r c) | ||
40 | z <- createMatrix (orderOf m) r c | ||
41 | app3 (cDiag o) mat m vec y mat z "cDiag" | ||
42 | return (y,z) | ||
43 | where r = rows m | ||
44 | c = cols m | ||
45 | o = if orderOf m == RowMajor then 1 else 0 | ||
diff --git a/hmatrix.cabal b/hmatrix.cabal index 8235753..e185d8a 100644 --- a/hmatrix.cabal +++ b/hmatrix.cabal | |||
@@ -11,7 +11,7 @@ Description: This library provides a purely functional interface to basic | |||
11 | and other numerical computations, internally implemented using | 11 | and other numerical computations, internally implemented using |
12 | GSL, BLAS and LAPACK. | 12 | GSL, BLAS and LAPACK. |
13 | Category: Math | 13 | Category: Math |
14 | tested-with: GHC ==6.10.2 | 14 | tested-with: GHC ==6.10.3 |
15 | 15 | ||
16 | cabal-version: >=1.2 | 16 | cabal-version: >=1.2 |
17 | build-type: Custom | 17 | build-type: Custom |
@@ -20,6 +20,23 @@ extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h | |||
20 | extra-source-files: configure configure.hs README INSTALL | 20 | extra-source-files: configure configure.hs README INSTALL |
21 | extra-tmp-files: hmatrix.buildinfo | 21 | extra-tmp-files: hmatrix.buildinfo |
22 | 22 | ||
23 | extra-source-files: examples/tests.hs | ||
24 | examples/deriv.hs | ||
25 | examples/integrate.hs | ||
26 | examples/minimize.hs | ||
27 | examples/pca1.hs | ||
28 | examples/pca2.hs | ||
29 | examples/pinv.hs | ||
30 | examples/lie.hs | ||
31 | examples/kalman.hs | ||
32 | examples/plot.hs | ||
33 | examples/latexmat.hs | ||
34 | examples/inplace.hs | ||
35 | examples/benchmarks.hs | ||
36 | examples/error.hs | ||
37 | examples/devel/wrappers.hs | ||
38 | examples/devel/functions.hs | ||
39 | |||
23 | flag splitBase | 40 | flag splitBase |
24 | description: Choose the new smaller, split-up base package. | 41 | description: Choose the new smaller, split-up base package. |
25 | 42 | ||
@@ -90,9 +107,10 @@ library | |||
90 | Numeric.LinearAlgebra.Interface, | 107 | Numeric.LinearAlgebra.Interface, |
91 | Numeric.LinearAlgebra.Algorithms, | 108 | Numeric.LinearAlgebra.Algorithms, |
92 | Graphics.Plot, | 109 | Graphics.Plot, |
93 | Numeric.LinearAlgebra.Tests | 110 | Numeric.LinearAlgebra.Tests, |
94 | Data.Packed.Convert | 111 | Data.Packed.Convert, |
95 | Data.Packed.ST | 112 | Data.Packed.ST, |
113 | Data.Packed.Development | ||
96 | other-modules: Data.Packed.Internal, | 114 | other-modules: Data.Packed.Internal, |
97 | Data.Packed.Internal.Common, | 115 | Data.Packed.Internal.Common, |
98 | Data.Packed.Internal.Vector, | 116 | Data.Packed.Internal.Vector, |
diff --git a/lib/Data/Packed/Development.hs b/lib/Data/Packed/Development.hs new file mode 100644 index 0000000..cbb2f81 --- /dev/null +++ b/lib/Data/Packed/Development.hs | |||
@@ -0,0 +1,27 @@ | |||
1 | |||
2 | ----------------------------------------------------------------------------- | ||
3 | -- | | ||
4 | -- Module : Data.Packed.Development | ||
5 | -- Copyright : (c) Alberto Ruiz 2009 | ||
6 | -- License : GPL | ||
7 | -- | ||
8 | -- Maintainer : Alberto Ruiz <aruiz@um.es> | ||
9 | -- Stability : provisional | ||
10 | -- Portability : portable | ||
11 | -- | ||
12 | -- The implementation of the 'Vector' and 'Matrix' types is not exposed, | ||
13 | -- but the library can be easily extended with additional foreign functions | ||
14 | -- using the tools in this module. Illustrative usage examples can be found | ||
15 | -- in the @examples\/devel@ folder included in the package. | ||
16 | -- | ||
17 | ----------------------------------------------------------------------------- | ||
18 | |||
19 | module Data.Packed.Development ( | ||
20 | createVector, createMatrix, | ||
21 | Adapt, | ||
22 | vec, mat, | ||
23 | app1, app2, app3, app4, | ||
24 | MatrixOrder(..), orderOf, cmat, fmat, | ||
25 | ) where | ||
26 | |||
27 | import Data.Packed.Internal \ No newline at end of file | ||
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs index 2310f5f..879dade 100644 --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs | |||
@@ -60,6 +60,44 @@ ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 | |||
60 | ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ \a1 -> ww2 w2 o2 w3 o3 (f a1) | 60 | ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ \a1 -> ww2 w2 o2 w3 o3 (f a1) |
61 | ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ \a1 -> ww3 w2 o2 w3 o3 w4 o4 (f a1) | 61 | ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ \a1 -> ww3 w2 o2 w3 o3 w4 o4 (f a1) |
62 | 62 | ||
63 | type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() | ||
64 | |||
65 | app1 :: f | ||
66 | -> Adapt f t (IO CInt) | ||
67 | -> t | ||
68 | -> String | ||
69 | -> IO() | ||
70 | |||
71 | app2 :: f | ||
72 | -> Adapt f t1 r | ||
73 | -> t1 | ||
74 | -> Adapt r t2 (IO CInt) | ||
75 | -> t2 | ||
76 | -> String | ||
77 | -> IO() | ||
78 | |||
79 | app3 :: f | ||
80 | -> Adapt f t1 r1 | ||
81 | -> t1 | ||
82 | -> Adapt r1 t2 r2 | ||
83 | -> t2 | ||
84 | -> Adapt r2 t3 (IO CInt) | ||
85 | -> t3 | ||
86 | -> String | ||
87 | -> IO() | ||
88 | |||
89 | app4 :: f | ||
90 | -> Adapt f t1 r1 | ||
91 | -> t1 | ||
92 | -> Adapt r1 t2 r2 | ||
93 | -> t2 | ||
94 | -> Adapt r2 t3 r3 | ||
95 | -> t3 | ||
96 | -> Adapt r3 t4 (IO CInt) | ||
97 | -> t4 | ||
98 | -> String | ||
99 | -> IO() | ||
100 | |||
63 | app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s | 101 | app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s |
64 | app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s | 102 | app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s |
65 | app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $ | 103 | app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $ |
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 13ffc34..8a074a6 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -74,6 +74,7 @@ data Matrix t = MC { rows :: {-# UNPACK #-} !Int | |||
74 | xdat MC {cdat = d } = d | 74 | xdat MC {cdat = d } = d |
75 | xdat MF {fdat = d } = d | 75 | xdat MF {fdat = d } = d |
76 | 76 | ||
77 | orderOf :: Matrix t -> MatrixOrder | ||
77 | orderOf MF{} = ColumnMajor | 78 | orderOf MF{} = ColumnMajor |
78 | orderOf MC{} = RowMajor | 79 | orderOf MC{} = RowMajor |
79 | 80 | ||
@@ -82,12 +83,16 @@ trans :: Matrix t -> Matrix t | |||
82 | trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d } | 83 | trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d } |
83 | trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d } | 84 | trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d } |
84 | 85 | ||
86 | cmat :: (Element t) => Matrix t -> Matrix t | ||
85 | cmat m@MC{} = m | 87 | cmat m@MC{} = m |
86 | cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c} | 88 | cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c} |
87 | 89 | ||
90 | fmat :: (Element t) => Matrix t -> Matrix t | ||
88 | fmat m@MF{} = m | 91 | fmat m@MF{} = m |
89 | fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} | 92 | fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} |
90 | 93 | ||
94 | -- C-Haskell matrix adapter | ||
95 | mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r | ||
91 | mat = withMatrix | 96 | mat = withMatrix |
92 | 97 | ||
93 | withMatrix a f = | 98 | withMatrix a f = |
@@ -156,7 +161,7 @@ MF {rows = r, cols = c, fdat = v} @@> (i,j) | |||
156 | | otherwise = v `at` (j*r+i) | 161 | | otherwise = v `at` (j*r+i) |
157 | {-# INLINE (@@>) #-} | 162 | {-# INLINE (@@>) #-} |
158 | 163 | ||
159 | -- | Unsafe matrix access without range checking | 164 | -- Unsafe matrix access without range checking |
160 | atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j) | 165 | atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j) |
161 | atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i) | 166 | atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i) |
162 | {-# INLINE atM' #-} | 167 | {-# INLINE atM' #-} |
@@ -173,6 +178,8 @@ matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v } | |||
173 | r | m==0 = d | 178 | r | m==0 = d |
174 | | otherwise = error "matrixFromVector" | 179 | | otherwise = error "matrixFromVector" |
175 | 180 | ||
181 | -- allocates memory for a new matrix | ||
182 | createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a) | ||
176 | createMatrix order r c = do | 183 | createMatrix order r c = do |
177 | p <- createVector (r*c) | 184 | p <- createVector (r*c) |
178 | return (matrixFromVector order c p) | 185 | return (matrixFromVector order c p) |
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs index fc8a6be..1b572a5 100644 --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs | |||
@@ -18,6 +18,7 @@ module Data.Packed.Internal.Vector where | |||
18 | 18 | ||
19 | import Data.Packed.Internal.Common | 19 | import Data.Packed.Internal.Common |
20 | import Foreign | 20 | import Foreign |
21 | import Foreign.C.Types(CInt) | ||
21 | import Complex | 22 | import Complex |
22 | import Control.Monad(when) | 23 | import Control.Monad(when) |
23 | 24 | ||
@@ -36,6 +37,8 @@ data Vector t = | |||
36 | , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block | 37 | , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block |
37 | } | 38 | } |
38 | 39 | ||
40 | -- C-Haskell vector adapter | ||
41 | vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r | ||
39 | vec = withVector | 42 | vec = withVector |
40 | 43 | ||
41 | withVector (V n fp) f = withForeignPtr fp $ \p -> do | 44 | withVector (V n fp) f = withForeignPtr fp $ \p -> do |
@@ -43,7 +46,7 @@ withVector (V n fp) f = withForeignPtr fp $ \p -> do | |||
43 | g (fi n) p | 46 | g (fi n) p |
44 | f v | 47 | f v |
45 | 48 | ||
46 | -- | allocates memory for a new vector | 49 | -- allocates memory for a new vector |
47 | createVector :: Storable a => Int -> IO (Vector a) | 50 | createVector :: Storable a => Int -> IO (Vector a) |
48 | createVector n = do | 51 | createVector n = do |
49 | when (n <= 0) $ error ("trying to createVector of dim "++show n) | 52 | when (n <= 0) $ error ("trying to createVector of dim "++show n) |