4 files changed, 135 insertions, 0 deletions

diff --git a/examples/devel/ej1/functions.c b/examples/devel/ej1/functions.c
new file mode 100644
index 0000000..02a4cdd
--- /dev/null
+++ b/examples/devel/ej1/functions.c
@@ -0,0 +1,35 @@
+/* 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
new file mode 100644
index 0000000..a88f74b
--- /dev/null
+++ b/examples/devel/ej1/wrappers.hs
@@ -0,0 +1,44 @@
+{-# 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
new file mode 100644
index 0000000..4dcd377
--- /dev/null
+++ b/examples/devel/ej2/functions.c
@@ -0,0 +1,24 @@
+/*  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
new file mode 100644
index 0000000..1c02a24
--- /dev/null
+++ b/examples/devel/ej2/wrappers.hs
@@ -0,0 +1,32 @@
+{-# 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