1 files changed, 90 insertions, 0 deletions

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;
+}