common trans

3 files changed, 123 insertions, 0 deletions

diff --git a/lib/Data/Packed/Internal.hs b/lib/Data/Packed/Internal.hs
index c8ad8d7..5e19e58 100644
--- a/lib/Data/Packed/Internal.hs
+++ b/lib/Data/Packed/Internal.hs
@@ -67,6 +67,15 @@ infixl 0 //
 vec :: Vector a -> (Int -> Ptr b -> t) -> t
 vec v f = f (dim v) (castPtr $ ptr v)
 
+mata :: Matrix a -> (Int-> Int -> Ptr b -> t) -> t
+mata m f = f (rows m) (cols m) (castPtr $ ptr (mat m))
+
+pd2pc :: Ptr Double -> Ptr (Complex (Double))
+pd2pc = castPtr
+
+pc2pd :: Ptr (Complex (Double)) -> Ptr Double
+pc2pd = castPtr
+
 check msg ls f = do
     err <- f
     when (err/=0) (error msg)
@@ -97,6 +106,9 @@ at :: Storable a => Vector a -> Int -> a
 at v n | n >= 0 && n < dim v = at' v n
        | otherwise          = error "vector index out of range"
 
+dsv v = sizeOf (v `at` 0)
+dsm m = (dsv.mat) m
+
 constant :: Storable a => Int -> a -> Vector a
 constant n x = unsafePerformIO $ do
     v <- createVector n
@@ -118,3 +130,18 @@ reshape n v = M { rows = dim v `div` n
                 , trMode = NoTrans
                 , isCOrder = True
                 }
+
+createMatrix r c = do
+    p <- createVector (r*c)
+    return (reshape c p)
+
+type CMat s = Int -> Int -> Ptr Double -> s
+type CVec s = Int -> Ptr Double -> s
+
+foreign import ccall safe "aux.h trans" ctrans :: Int -> CMat (CMat (IO Int))
+
+trans :: Storable a => Matrix a -> Matrix a
+trans m = unsafePerformIO $ do
+    r <- createMatrix (cols m) (rows m)
+    ctrans (dsm m) // mata m // mata r // check "trans" [mat m]
+    return r
diff --git a/lib/Data/Packed/aux.c b/lib/Data/Packed/aux.c
new file mode 100644
index 0000000..d772d90
--- /dev/null
+++ b/lib/Data/Packed/aux.c
@@ -0,0 +1,90 @@
+#include "aux.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_fft_complex.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_deriv.h>
+#include <gsl/gsl_poly.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <string.h>
+#include <stdio.h>
+
+#define MACRO(B) do {B} while (0)
+#define ERROR(CODE) MACRO(return CODE;)
+#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("GSL Wrapper "M": "); size_t t0 = time(NULL);
+#define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
+#else
+#define DEBUGMSG(M)
+#define OK return 0;
+#endif
+
+
+#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
+
+#ifdef DBG
+#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGMAT(MSG,X)
+#endif
+
+#ifdef DBG
+#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGVEC(MSG,X)
+#endif
+
+
+#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
+#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
+#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array(A##p,A##n)
+#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array(A##p,A##r,A##c)
+#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
+#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
+#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array(A##p,A##n)
+#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array(A##p,A##r,A##c)
+
+#define V(a) (&a.vector)
+#define M(a) (&a.matrix)
+
+#define GCVEC(A) int A##n, gsl_complex*A##p
+#define KGCVEC(A) int A##n, const gsl_complex*A##p
+
+#define BAD_SIZE 1000
+#define BAD_CODE 1001
+#define MEM      1002
+#define BAD_FILE 1003
+#define BAD_TYPE 1004
+
+
+int trans(int size,KMAT(x),MAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("trans");
+    if(size==8) {
+        DEBUGMSG("trans double");
+        KDMVIEW(x);
+        DMVIEW(t);
+        int res = gsl_matrix_transpose_memcpy(M(t),M(x));
+        CHECK(res,res);
+        OK
+    } else if (size==16) {
+        DEBUGMSG("trans complex double");
+        KCMVIEW(x);
+        CMVIEW(t);
+        int res = gsl_matrix_complex_transpose_memcpy(M(t),M(x));
+        CHECK(res,res);
+        OK
+    }
+    return BAD_TYPE;
+}
diff --git a/lib/Data/Packed/aux.h b/lib/Data/Packed/aux.h
new file mode 100644
index 0000000..c51234a
--- /dev/null
+++ b/lib/Data/Packed/aux.h
@@ -0,0 +1,6 @@
+#define VEC(A) int A##n, double*A##p
+#define MAT(A) int A##r, int A##c, double* A##p
+#define KVEC(A) int A##n, const double*A##p
+#define KMAT(A) int A##r, int A##c, const double* A##p
+
+int trans(int size, KMAT(x),MAT(t));