summaryrefslogtreecommitdiff
path: root/examples/devel
diff options
context:
space:
mode:
Diffstat (limited to 'examples/devel')
-rw-r--r--examples/devel/ej1/functions.c35
-rw-r--r--examples/devel/ej1/wrappers.hs44
-rw-r--r--examples/devel/ej2/functions.c24
-rw-r--r--examples/devel/ej2/wrappers.hs32
-rw-r--r--examples/devel/example/functions.c22
-rw-r--r--examples/devel/example/wrappers.hs45
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 @@
1/* assuming row order */
2
3typedef struct { double r, i; } doublecomplex;
4
5#define DVEC(A) int A##n, double*A##p
6#define CVEC(A) int A##n, doublecomplex*A##p
7#define DMAT(A) int A##r, int A##c, double*A##p
8#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
9
10#define AT(M,row,col) (M##p[(row)*M##c + (col)])
11
12/*-----------------------------------------------------*/
13
14int c_scale_vector(double s, DVEC(x), DVEC(y)) {
15 int k;
16 for (k=0; k<=yn; k++) {
17 yp[k] = s*xp[k];
18 }
19 return 0;
20}
21
22/*-----------------------------------------------------*/
23
24int c_diag(DMAT(m),DVEC(y),DMAT(z)) {
25 int i,j;
26 for (j=0; j<yn; j++) {
27 yp[j] = AT(m,j,j);
28 }
29 for (i=0; i<mr; i++) {
30 for (j=0; j<mc; j++) {
31 AT(z,i,j) = i==j?yp[i]:0;
32 }
33 }
34 return 0;
35}
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 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3-- $ ghc -O2 --make wrappers.hs functions.c
4
5import Numeric.LinearAlgebra
6import Data.Packed.Development
7import Foreign(Ptr,unsafePerformIO)
8import Foreign.C.Types(CInt)
9
10-----------------------------------------------------
11
12main = do
13 print $ myScale 3.0 (fromList [1..10])
14 print $ myDiag $ (3><5) [1..]
15
16-----------------------------------------------------
17
18foreign import ccall unsafe "c_scale_vector"
19 cScaleVector :: Double -- scale
20 -> CInt -> Ptr Double -- argument
21 -> CInt -> Ptr Double -- result
22 -> IO CInt -- exit code
23
24myScale s x = unsafePerformIO $ do
25 y <- createVector (dim x)
26 app2 (cScaleVector s) vec x vec y "cScaleVector"
27 return y
28
29-----------------------------------------------------
30-- forcing row order
31
32foreign import ccall unsafe "c_diag"
33 cDiag :: CInt -> CInt -> Ptr Double -- argument
34 -> CInt -> Ptr Double -- result1
35 -> CInt -> CInt -> Ptr Double -- result2
36 -> IO CInt -- exit code
37
38myDiag m = unsafePerformIO $ do
39 y <- createVector (min r c)
40 z <- createMatrix RowMajor r c
41 app3 cDiag mat (cmat m) vec y mat z "cDiag"
42 return (y,z)
43 where r = rows m
44 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 @@
1/* general element order */
2
3typedef struct { double r, i; } doublecomplex;
4
5#define DVEC(A) int A##n, double*A##p
6#define CVEC(A) int A##n, doublecomplex*A##p
7#define DMAT(A) int A##r, int A##c, double*A##p
8#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
9
10#define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
11
12int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) {
13 int i,j,sr,sc;
14 if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
15 for (j=0; j<yn; j++) {
16 yp[j] = AT(m,j,j);
17 }
18 for (i=0; i<mr; i++) {
19 for (j=0; j<mc; j++) {
20 AT(z,i,j) = i==j?yp[i]:0;
21 }
22 }
23 return 0;
24}
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 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3-- $ ghc -O2 --make wrappers.hs functions.c
4
5import Numeric.LinearAlgebra
6import Data.Packed.Development
7import Foreign(Ptr,unsafePerformIO)
8import Foreign.C.Types(CInt)
9
10-----------------------------------------------------
11
12main = do
13 print $ myDiag $ (3><5) [1..]
14
15-----------------------------------------------------
16-- arbitrary data order
17
18foreign import ccall unsafe "c_diag"
19 cDiag :: CInt -- matrix order
20 -> CInt -> CInt -> Ptr Double -- argument
21 -> CInt -> Ptr Double -- result1
22 -> CInt -> CInt -> Ptr Double -- result2
23 -> IO CInt -- exit code
24
25myDiag m = unsafePerformIO $ do
26 y <- createVector (min r c)
27 z <- createMatrix (orderOf m) r c
28 app3 (cDiag o) mat m vec y mat z "cDiag"
29 return (y,z)
30 where r = rows m
31 c = cols m
32 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 @@
1
2typedef struct { double r, i; } doublecomplex;
3
4#define VEC(T,A) int A##n, T* A##p
5#define MAT(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
6
7#define AT(m,i,j) (m##p[(i)*m##Xr + (j)*m##Xc])
8#define TRAV(m,i,j) int i,j; for (i=0;i<m##r;i++) for (j=0;j<m##c;j++)
9
10
11int c_diag(MAT(double,m), VEC(double,y), MAT(double,z)) {
12 int k;
13 for (k=0; k<yn; k++) {
14 yp[k] = AT(m,k,k);
15 }
16 { TRAV(z,i,j) {
17 AT(z,i,j) = i==j?yp[i]:0;
18 }
19 }
20 return 0;
21}
22
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 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2{-# LANGUAGE TypeOperators #-}
3{-# LANGUAGE GADTs #-}
4
5{-
6 $ ghc -O2 wrappers.hs functions.c
7 $ ./wrappers
8-}
9
10import Numeric.LinearAlgebra
11import Numeric.LinearAlgebra.Devel
12import System.IO.Unsafe(unsafePerformIO)
13import Foreign.C.Types(CInt(..))
14import Foreign.Ptr(Ptr)
15
16
17infixl 1 #
18a # b = apply a b
19{-# INLINE (#) #-}
20
21infixr 5 :>, ::>
22type (:>) t r = CInt -> Ptr t -> r
23type (::>) t r = CInt -> CInt -> CInt -> CInt -> Ptr t -> r
24type Ok = IO CInt
25
26-----------------------------------------------------
27
28x = (3><5) [1..]
29
30main = do
31 print x
32 print $ myDiag x
33 print $ myDiag (tr x)
34
35-----------------------------------------------------
36foreign import ccall unsafe "c_diag" cDiag :: Double ::> Double :> Double ::> Ok
37
38myDiag m = unsafePerformIO $ do
39 y <- createVector (min r c)
40 z <- createMatrix RowMajor r c
41 cDiag # m # y # z #| "cDiag"
42 return (y,z)
43 where
44 (r,c) = size m
45
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 02:56:40 +0000