From 6809103cf34a9345f8cb60a0ec3a8f55dd18d5ef Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 13 May 2009 09:39:25 +0000 Subject: added Development module --- examples/devel/functions.c | 34 ++++++++++++++++++++++++++++ examples/devel/wrappers.hs | 45 ++++++++++++++++++++++++++++++++++++++ hmatrix.cabal | 26 ++++++++++++++++++---- lib/Data/Packed/Development.hs | 27 +++++++++++++++++++++++ lib/Data/Packed/Internal/Common.hs | 38 ++++++++++++++++++++++++++++++++ lib/Data/Packed/Internal/Matrix.hs | 9 +++++++- lib/Data/Packed/Internal/Vector.hs | 5 ++++- 7 files changed, 178 insertions(+), 6 deletions(-) create mode 100644 examples/devel/functions.c create mode 100644 examples/devel/wrappers.hs create mode 100644 lib/Data/Packed/Development.hs 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 @@ +typedef struct { double r, i; } doublecomplex; + +#define DVEC(A) int A##n, double*A##p +#define CVEC(A) int A##n, doublecomplex*A##p +#define DMAT(A) int A##r, int A##c, double*A##p +#define CMAT(A) int A##r, int A##c, doublecomplex*A##p + +#define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) + +/*-----------------------------------------------------*/ + +int c_scale_vector(double s, DVEC(x), DVEC(y)) { + int k; + for (k=0; k<=yn; k++) { + yp[k] = s*xp[k]; + } + return 0; +} + +/*-----------------------------------------------------*/ + +int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) { + int i,j,sr,sc; + if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} + for (j=0; j<5) [1..] + +----------------------------------------------------- + +foreign import ccall "c_scale_vector" + cScaleVector :: Double -- scale + -> CInt -> Ptr Double -- argument + -> CInt -> Ptr Double -- result + -> IO CInt -- exit code + +myScale s x = unsafePerformIO $ do + y <- createVector (dim x) + app2 (cScaleVector s) vec x vec y "cScaleVector" + return y + +----------------------------------------------------- + +foreign import ccall "c_diag" + cDiag :: CInt -- matrix order + -> CInt -> CInt -> Ptr Double -- argument + -> CInt -> Ptr Double -- result1 + -> CInt -> CInt -> Ptr Double -- result2 + -> IO CInt -- exit code + +myDiag m = unsafePerformIO $ do + y <- createVector (min r c) + z <- createMatrix (orderOf m) r c + app3 (cDiag o) mat m vec y mat z "cDiag" + return (y,z) + where r = rows m + c = cols m + 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 and other numerical computations, internally implemented using GSL, BLAS and LAPACK. Category: Math -tested-with: GHC ==6.10.2 +tested-with: GHC ==6.10.3 cabal-version: >=1.2 build-type: Custom @@ -20,6 +20,23 @@ extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h extra-source-files: configure configure.hs README INSTALL extra-tmp-files: hmatrix.buildinfo +extra-source-files: examples/tests.hs + examples/deriv.hs + examples/integrate.hs + examples/minimize.hs + examples/pca1.hs + examples/pca2.hs + examples/pinv.hs + examples/lie.hs + examples/kalman.hs + examples/plot.hs + examples/latexmat.hs + examples/inplace.hs + examples/benchmarks.hs + examples/error.hs + examples/devel/wrappers.hs + examples/devel/functions.hs + flag splitBase description: Choose the new smaller, split-up base package. @@ -90,9 +107,10 @@ library Numeric.LinearAlgebra.Interface, Numeric.LinearAlgebra.Algorithms, Graphics.Plot, - Numeric.LinearAlgebra.Tests - Data.Packed.Convert - Data.Packed.ST + Numeric.LinearAlgebra.Tests, + Data.Packed.Convert, + Data.Packed.ST, + Data.Packed.Development other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, 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 @@ + +----------------------------------------------------------------------------- +-- | +-- Module : Data.Packed.Development +-- Copyright : (c) Alberto Ruiz 2009 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- The implementation of the 'Vector' and 'Matrix' types is not exposed, +-- but the library can be easily extended with additional foreign functions +-- using the tools in this module. Illustrative usage examples can be found +-- in the @examples\/devel@ folder included in the package. +-- +----------------------------------------------------------------------------- + +module Data.Packed.Development ( + createVector, createMatrix, + Adapt, + vec, mat, + app1, app2, app3, app4, + MatrixOrder(..), orderOf, cmat, fmat, +) where + +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 ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ \a1 -> ww2 w2 o2 w3 o3 (f a1) ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ \a1 -> ww3 w2 o2 w3 o3 w4 o4 (f a1) +type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() + +app1 :: f + -> Adapt f t (IO CInt) + -> t + -> String + -> IO() + +app2 :: f + -> Adapt f t1 r + -> t1 + -> Adapt r t2 (IO CInt) + -> t2 + -> String + -> IO() + +app3 :: f + -> Adapt f t1 r1 + -> t1 + -> Adapt r1 t2 r2 + -> t2 + -> Adapt r2 t3 (IO CInt) + -> t3 + -> String + -> IO() + +app4 :: f + -> Adapt f t1 r1 + -> t1 + -> Adapt r1 t2 r2 + -> t2 + -> Adapt r2 t3 r3 + -> t3 + -> Adapt r3 t4 (IO CInt) + -> t4 + -> String + -> IO() + app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s 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 xdat MC {cdat = d } = d xdat MF {fdat = d } = d +orderOf :: Matrix t -> MatrixOrder orderOf MF{} = ColumnMajor orderOf MC{} = RowMajor @@ -82,12 +83,16 @@ trans :: Matrix t -> Matrix t trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d } trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d } +cmat :: (Element t) => Matrix t -> Matrix t cmat m@MC{} = m cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c} +fmat :: (Element t) => Matrix t -> Matrix t fmat m@MF{} = m fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} +-- C-Haskell matrix adapter +mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r mat = withMatrix withMatrix a f = @@ -156,7 +161,7 @@ MF {rows = r, cols = c, fdat = v} @@> (i,j) | otherwise = v `at` (j*r+i) {-# INLINE (@@>) #-} --- | Unsafe matrix access without range checking +-- Unsafe matrix access without range checking atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j) atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i) {-# INLINE atM' #-} @@ -173,6 +178,8 @@ matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v } r | m==0 = d | otherwise = error "matrixFromVector" +-- allocates memory for a new matrix +createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a) createMatrix order r c = do p <- createVector (r*c) 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 import Data.Packed.Internal.Common import Foreign +import Foreign.C.Types(CInt) import Complex import Control.Monad(when) @@ -36,6 +37,8 @@ data Vector t = , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block } +-- C-Haskell vector adapter +vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r vec = withVector withVector (V n fp) f = withForeignPtr fp $ \p -> do @@ -43,7 +46,7 @@ withVector (V n fp) f = withForeignPtr fp $ \p -> do g (fi n) p f v --- | allocates memory for a new vector +-- allocates memory for a new vector createVector :: Storable a => Int -> IO (Vector a) createVector n = do when (n <= 0) $ error ("trying to createVector of dim "++show n) -- cgit v1.2.3