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 @@
+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<yn; j++) {
+        yp[j] = AT(m,j,j);
+    }
+    for (i=0; i<mr; i++) {
+        for (j=0; j<mc; j++) {
+            AT(z,i,j) = i==j?yp[i]:0;
+        }
+    }
+    return 0;
+}
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 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- $ ghc -O2 --make wrappers.hs functions.c
+
+import Numeric.LinearAlgebra
+import Data.Packed.Development
+import Foreign(Ptr,unsafePerformIO)
+import Foreign.C.Types(CInt)
+
+-----------------------------------------------------
+
+main = do
+    print $ myScale 3.0 (fromList [1..10])
+    print $ myDiag $ (3><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 <aruiz@um.es>
+-- 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)