updated devel example

6 files changed, 67 insertions, 135 deletions

diff --git a/examples/devel/ej1/functions.c b/examples/devel/ej1/functions.c
deleted file mode 100644
index 02a4cdd..0000000
--- a/examples/devel/ej1/functions.c
+++ /dev/null
@@ -1,35 +0,0 @@
-/* assuming row order */
-
-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,row,col) (M##p[(row)*M##c + (col)])
-
-/*-----------------------------------------------------*/
-
-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(DMAT(m),DVEC(y),DMAT(z)) {
-    int i,j;
-    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/ej1/wrappers.hs b/examples/devel/ej1/wrappers.hs
deleted file mode 100644
index a88f74b..0000000
--- a/examples/devel/ej1/wrappers.hs
+++ /dev/null
@@ -1,44 +0,0 @@
-{-# 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 unsafe "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
-
------------------------------------------------------
--- forcing row order
-
-foreign import ccall unsafe "c_diag"
-    cDiag :: 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 RowMajor r c
-    app3 cDiag mat (cmat m) vec y mat z "cDiag"
-    return (y,z)
-  where r = rows m
-        c = cols m
diff --git a/examples/devel/ej2/functions.c b/examples/devel/ej2/functions.c
deleted file mode 100644
index 4dcd377..0000000
--- a/examples/devel/ej2/functions.c
+++ /dev/null
@@ -1,24 +0,0 @@
-/*  general element order   */
-
-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_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/ej2/wrappers.hs b/examples/devel/ej2/wrappers.hs
deleted file mode 100644
index 1c02a24..0000000
--- a/examples/devel/ej2/wrappers.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# 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 $ myDiag $ (3><5) [1..]
-
------------------------------------------------------
--- arbitrary data order
-
-foreign import ccall unsafe "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/examples/devel/example/functions.c b/examples/devel/example/functions.c
new file mode 100644
index 0000000..67d3270
--- /dev/null
+++ b/examples/devel/example/functions.c
@@ -0,0 +1,22 @@
+
+typedef struct { double r, i; } doublecomplex;
+
+#define VEC(T,A) int A##n, T* A##p
+#define MAT(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
+
+#define AT(m,i,j) (m##p[(i)*m##Xr + (j)*m##Xc])
+#define TRAV(m,i,j) int i,j; for (i=0;i<m##r;i++) for (j=0;j<m##c;j++)
+
+
+int c_diag(MAT(double,m), VEC(double,y), MAT(double,z)) {
+    int k;
+    for (k=0; k<yn; k++) {
+        yp[k] = AT(m,k,k);
+    }
+    { TRAV(z,i,j) {
+        AT(z,i,j) = i==j?yp[i]:0;
+        }
+    }
+    return 0;
+}
+
diff --git a/examples/devel/example/wrappers.hs b/examples/devel/example/wrappers.hs
new file mode 100644
index 0000000..f4e0f0b
--- /dev/null
+++ b/examples/devel/example/wrappers.hs
@@ -0,0 +1,45 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GADTs #-}
+
+{-
+    $ ghc -O2 wrappers.hs functions.c
+    $ ./wrappers
+-}
+
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Devel
+import System.IO.Unsafe(unsafePerformIO)
+import Foreign.C.Types(CInt(..))
+import Foreign.Ptr(Ptr)
+
+
+infixl 1 #
+a # b = apply a b
+{-# INLINE (#) #-}
+
+infixr 5 :>, ::>
+type (:>)  t r = CInt -> Ptr t -> r
+type (::>) t r =  CInt -> CInt -> CInt -> CInt -> Ptr t -> r
+type Ok = IO CInt
+
+-----------------------------------------------------
+
+x = (3><5) [1..]
+
+main = do
+    print   x
+    print $ myDiag x
+    print $ myDiag (tr x)
+
+-----------------------------------------------------
+foreign import ccall unsafe "c_diag" cDiag :: Double ::> Double :> Double ::> Ok
+
+myDiag m = unsafePerformIO $ do
+    y <- createVector (min r c)
+    z <- createMatrix RowMajor r c
+    cDiag # m # y # z #| "cDiag"
+    return (y,z)
+  where
+    (r,c) = size m
+