vectorized operations in base

7 files changed, 1595 insertions, 631 deletions

diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index bd356df..cf84f30 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -37,6 +37,7 @@ library
                         Data.Packed.Development,
                         Numeric.Conversion
                         Numeric.LinearAlgebra.LAPACK
+                        Numeric.Vectorized
                         
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
@@ -45,6 +46,7 @@ library
                         Data.Packed.Internal.Matrix
 
     C-sources:          src/C/lapack-aux.c
+                        src/C/vector-aux.c
                         
 
     extensions:         ForeignFunctionInterface,
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
new file mode 100644
index 0000000..07ccf76
--- /dev/null
+++ b/packages/base/src/C/vector-aux.c
@@ -0,0 +1,676 @@
+#include "lapack-aux.h"
+
+#define V(x) x##n,x##p
+
+#include <string.h>
+#include <math.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 OK return 0;
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
+#else
+#define DEBUGMSG(M)
+#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 BAD_SIZE 2000
+#define BAD_CODE 2001
+#define MEM      2002
+#define BAD_FILE 2003
+
+
+int sumF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("sumF");
+    printf("hello!\n");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int sumR(KDVEC(x),DVEC(r)) {
+    DEBUGMSG("sumR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+
+
+int sumQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("sumQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    complex res;
+    res.r = 0;
+    res.i = 0;
+    for (i = 0; i < xn; i++) {
+      res.r += xp[i].r;
+      res.i += xp[i].i;
+    }
+    rp[0] = res;
+    OK
+}
+    
+int sumC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("sumC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    doublecomplex res;
+    res.r = 0;
+    res.i = 0;
+    for (i = 0; i < xn; i++)  {
+      res.r += xp[i].r;
+      res.i += xp[i].i;
+    }
+    rp[0] = res;
+    OK
+}
+
+
+int prodF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("prodF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int prodR(KDVEC(x),DVEC(r)) {
+    DEBUGMSG("prodR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+
+
+int prodQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("prodQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    complex res;
+    float temp;
+    res.r = 1;
+    res.i = 0;
+    for (i = 0; i < xn; i++) {
+      temp  = res.r * xp[i].r - res.i * xp[i].i;
+      res.i = res.r * xp[i].i + res.i * xp[i].r;
+      res.r = temp;
+    }
+    rp[0] = res;
+    OK
+}
+    
+int prodC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("prodC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    doublecomplex res;
+    double temp;
+    res.r = 1;
+    res.i = 0;
+    for (i = 0; i < xn; i++)  {
+      temp  = res.r * xp[i].r - res.i * xp[i].i;
+      res.i = res.r * xp[i].i + res.i * xp[i].r;
+      res.r = temp;
+    }
+    rp[0] = res;
+    OK
+}
+
+    
+double dnrm2_(integer*, const double*, integer*);
+double dasum_(integer*, const double*, integer*);
+
+double vector_max(KDVEC(x)) {
+    double r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+double vector_min(KDVEC(x)) {
+    double r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+double vector_max_index(KDVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[0]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+double vector_min_index(KDVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[0]) {
+            r = k;
+        }
+    }
+    return r;
+}
+   
+int toScalarR(int code, KDVEC(x), DVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarR");
+    double res;
+    integer one = 1;
+    integer n = rn;
+    switch(code) {
+        case 0: { res = dnrm2_(&n,xp,&one); break; }
+        case 1: { res = dasum_(&n,xp,&one);  break; }
+        case 2: { res = vector_max_index(V(x));  break; }
+        case 3: { res = vector_max(V(x));  break; }
+        case 4: { res = vector_min_index(V(x)); break; }
+        case 5: { res = vector_min(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+float snrm2_(integer*, const float*, integer*);
+float sasum_(integer*, const float*, integer*);
+
+float vector_max_f(KFVEC(x)) {
+    float r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+float vector_min_f(KFVEC(x)) {
+    float r = xp[0];
+    int k;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<r) {
+            r = xp[k];
+        }
+    }
+    return r;
+}
+
+float vector_max_index_f(KFVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]>xp[0]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+float vector_min_index_f(KFVEC(x)) {
+    int k, r = 0;
+    for (k = 1; k<xn; k++) {
+        if(xp[k]<xp[0]) {
+            r = k;
+        }
+    }
+    return r;
+}
+
+
+int toScalarF(int code, KFVEC(x), FVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarF");
+    float res;
+    integer one = 1;
+    integer n = rn;
+    switch(code) {
+        case 0: { res = snrm2_(&n,xp,&one); break; }
+        case 1: { res = sasum_(&n,xp,&one);  break; }
+        case 2: { res = vector_max_index_f(V(x));  break; }
+        case 3: { res = vector_max_f(V(x));  break; }
+        case 4: { res = vector_min_index_f(V(x)); break; }
+        case 5: { res = vector_min_f(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+double dznrm2_(integer*, const doublecomplex*, integer*);
+double dzasum_(integer*, const doublecomplex*, integer*);
+
+int toScalarC(int code, KCVEC(x), DVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarC");
+    double res;
+    integer one = 1;
+    integer n = rn;
+    switch(code) {
+        case 0: { res = dznrm2_(&n,xp,&one); break; }
+        case 1: { res = dzasum_(&n,xp,&one);  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+double scnrm2_(integer*, const complex*, integer*);
+double scasum_(integer*, const complex*, integer*);
+
+int toScalarQ(int code, KQVEC(x), FVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarQ");
+    float res;
+    integer one = 1;
+    integer n = rn;
+    switch(code) {
+        case 0: { res = scnrm2_(&n,xp,&one); break; }
+        case 1: { res = scasum_(&n,xp,&one);  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+inline double sign(double x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+inline float float_sign(float x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+
+#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
+#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
+int mapR(int code, KDVEC(x), DVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapR");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan) /* atan2 mediante vectorZip */
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+//        OP(10,gsl_asinh)
+//        OP(11,gsl_acosh)
+//        OP(12,gsl_atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapF(int code, KFVEC(x), FVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapF");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan) /* atan2 mediante vectorZip */
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+//        OP(10,gsl_asinh)
+//        OP(11,gsl_acosh)
+//        OP(12,gsl_atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+inline double abs_complex(doublecomplex z) {
+    return sqrt(z.r*z.r + z.i*z.i);
+}
+
+inline doublecomplex complex_abs_complex(doublecomplex z) {
+    doublecomplex r;
+    r.r = abs_complex(z);
+    r.i = 0;
+    return r;
+}
+
+inline doublecomplex complex_signum_complex(doublecomplex z) {
+    doublecomplex r;
+    double mag;
+    if (z.r == 0 && z.i == 0) {
+        r.r = 0;
+        r.i = 0;
+    } else {
+        mag = abs_complex(z);
+        r.r = z.r/mag;
+        r.i = z.i/mag;
+    }
+    return r;
+}
+
+
+
+int mapC(int code, KCVEC(x), CVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapC");
+    switch (code) {
+/*
+        OP(0,gsl_complex_sin)
+        OP(1,gsl_complex_cos)
+        OP(2,gsl_complex_tan)
+*/
+        OP(3,complex_abs_complex)
+/*
+        OP(4,gsl_complex_arcsin)
+        OP(5,gsl_complex_arccos)
+        OP(6,gsl_complex_arctan)
+        OP(7,gsl_complex_sinh)
+        OP(8,gsl_complex_cosh)
+        OP(9,gsl_complex_tanh)
+        OP(10,gsl_complex_arcsinh)
+        OP(11,gsl_complex_arccosh)
+        OP(12,gsl_complex_arctanh)
+        OP(13,gsl_complex_exp)
+        OP(14,gsl_complex_log)
+*/
+        OP(15,complex_signum_complex)
+        
+//      OP(16,gsl_complex_sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+inline complex complex_f_math_fun(doublecomplex (*cf)(doublecomplex), complex a)
+{
+  doublecomplex c;
+  doublecomplex r;
+
+  complex float_r;
+
+  c.r = a.r;
+  c.i = a.i;
+
+  r = (*cf)(c);
+
+  float_r.r = r.r;
+  float_r.i = r.i;
+
+  return float_r;
+}
+
+inline complex complex_f_math_op(doublecomplex (*cf)(doublecomplex,doublecomplex),
+                                 complex a,complex b)
+{
+  doublecomplex c1,c2,r;
+
+  complex float_r;
+
+  c1.r = a.r;
+  c1.i = a.i;
+
+  c2.r = b.r;
+  c2.i = b.i;
+
+  r = (*cf)(c1,c2);
+
+  float_r.r = r.r;
+  float_r.i = r.i;
+
+  return float_r;
+}
+
+#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_fun(&F,xp[k]); OK }
+#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_f_math_op(&F,A,B); OK }
+int mapQ(int code, KQVEC(x), QVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapQ");
+    switch (code) {
+/*
+        OPC(0,gsl_complex_sin)
+        OPC(1,gsl_complex_cos)
+        OPC(2,gsl_complex_tan)
+*/
+        OPC(3,complex_abs_complex)
+/*
+        OPC(4,gsl_complex_arcsin)
+        OPC(5,gsl_complex_arccos)
+        OPC(6,gsl_complex_arctan)
+        OPC(7,gsl_complex_sinh)
+        OPC(8,gsl_complex_cosh)
+        OPC(9,gsl_complex_tanh)
+        OPC(10,gsl_complex_arcsinh)
+        OPC(11,gsl_complex_arccosh)
+        OPC(12,gsl_complex_arctanh)
+        OPC(13,gsl_complex_exp)
+        OPC(14,gsl_complex_log)
+*/
+        OPC(15,complex_signum_complex)
+        
+//      OPC(16,gsl_complex_sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapValR(int code, double* pval, KDVEC(x), DVEC(r)) {
+    int k;
+    double val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValR");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
+    int k;
+    float val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValF");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+inline doublecomplex complex_add(doublecomplex a, doublecomplex b) {
+    doublecomplex r;
+    r.r = a.r+b.r;
+    r.i = a.i+b.i;
+    return r;
+}
+
+
+int mapValC(int code, doublecomplex* pval, KCVEC(x), CVEC(r)) {
+    int k;
+    doublecomplex val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValC");
+    switch (code) {
+//        OPV(0,gsl_complex_mul(val,xp[k]))
+//        OPV(1,gsl_complex_div(val,xp[k]))
+        OPV(2,complex_add(val,xp[k]))
+//        OPV(3,gsl_complex_sub(val,xp[k]))
+//        OPV(4,gsl_complex_pow(val,xp[k]))
+//        OPV(5,gsl_complex_pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapValQ(int code, complex* pval, KQVEC(x), QVEC(r)) {
+    int k;
+    complex val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValQ");
+    switch (code) {
+//        OPCA(0,gsl_complex_mul,val,xp[k])
+//        OPCA(1,gsl_complex_div,val,xp[k])
+        OPCA(2,complex_add,val,xp[k])
+//        OPCA(3,gsl_complex_sub,val,xp[k])
+//        OPCA(4,gsl_complex_pow,val,xp[k])
+//        OPCA(5,gsl_complex_pow,xp[k],val)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
+#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
+#define OPZO(C,msg,O) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = ap[k] O bp[k]; OK }
+
+int zipR(int code, KDVEC(a), KDVEC(b), DVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipR Add",+)
+        OPZO(1,"zipR Sub",-)
+        OPZO(2,"zipR Mul",*)
+        OPZO(3,"zipR Div",/)
+        OPZE(4,"zipR Pow",  pow)
+        OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
+REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZO(0,"zipR Add",+)
+        OPZO(1,"zipR Sub",-)
+        OPZO(2,"zipR Mul",*)
+        OPZO(3,"zipR Div",/)
+        OPZE(4,"zipR Pow",  pow)
+        OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+
+int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZE(0,"zipC Add",complex_add)
+//        OPZE(1,"zipC Sub",gsl_complex_sub)
+//        OPZE(2,"zipC Mul",gsl_complex_mul)
+//        OPZE(3,"zipC Div",gsl_complex_div)
+//        OPZE(4,"zipC Pow",gsl_complex_pow)
+//        //OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+
+
+
+#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_f_math_op(&E,ap[k],bp[k]); OK }
+
+int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPCZE(0,"zipQ Add",complex_add)
+//        OPCZE(1,"zipQ Sub",gsl_complex_sub)
+//        OPCZE(2,"zipQ Mul",gsl_complex_mul)
+//        OPCZE(3,"zipQ Div",gsl_complex_div)
+//        OPCZE(4,"zipQ Pow",gsl_complex_pow)
+//        //OPZE(5,"zipR ATan2",atan2)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
diff --git a/packages/base/src/Numeric/Vectorized.hs b/packages/base/src/Numeric/Vectorized.hs
new file mode 100644
index 0000000..3814579
--- /dev/null
+++ b/packages/base/src/Numeric/Vectorized.hs
@@ -0,0 +1,273 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Vectorized
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Low level interface to vector operations.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Vectorized (
+    sumF, sumR, sumQ, sumC,
+    prodF, prodR, prodQ, prodC,
+    FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
+    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
+    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
+    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ
+) where
+
+import Data.Packed.Internal.Common
+import Data.Packed.Internal.Signatures
+import Data.Packed.Internal.Vector
+
+import Data.Complex
+import Foreign.Marshal.Alloc(free)
+import Foreign.Marshal.Array(newArray)
+import Foreign.Ptr(Ptr)
+import Foreign.C.Types
+import System.IO.Unsafe(unsafePerformIO)
+
+fromei x = fromIntegral (fromEnum x) :: CInt
+
+data FunCodeV = Sin
+              | Cos
+              | Tan
+              | Abs
+              | ASin
+              | ACos
+              | ATan
+              | Sinh
+              | Cosh
+              | Tanh
+              | ASinh
+              | ACosh
+              | ATanh
+              | Exp
+              | Log
+              | Sign
+              | Sqrt
+              deriving Enum
+
+data FunCodeSV = Scale
+               | Recip
+               | AddConstant
+               | Negate
+               | PowSV
+               | PowVS
+               deriving Enum
+
+data FunCodeVV = Add
+               | Sub
+               | Mul
+               | Div
+               | Pow
+               | ATan2
+               deriving Enum
+
+data FunCodeS = Norm2
+              | AbsSum
+              | MaxIdx
+              | Max
+              | MinIdx
+              | Min
+              deriving Enum
+
+------------------------------------------------------------------
+
+-- | sum of elements
+sumF :: Vector Float -> Float
+sumF x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_sumF vec x vec r "sumF"
+           return $ r @> 0
+
+-- | sum of elements
+sumR :: Vector Double -> Double
+sumR x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_sumR vec x vec r "sumR"
+           return $ r @> 0
+
+-- | sum of elements
+sumQ :: Vector (Complex Float) -> Complex Float
+sumQ x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_sumQ vec x vec r "sumQ"
+           return $ r @> 0
+
+-- | sum of elements
+sumC :: Vector (Complex Double) -> Complex Double
+sumC x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_sumC vec x vec r "sumC"
+           return $ r @> 0
+
+foreign import ccall unsafe "sumF" c_sumF :: TFF
+foreign import ccall unsafe "sumR" c_sumR :: TVV
+foreign import ccall unsafe "sumQ" c_sumQ :: TQVQV
+foreign import ccall unsafe "sumC" c_sumC :: TCVCV
+
+-- | product of elements
+prodF :: Vector Float -> Float
+prodF x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodF vec x vec r "prodF"
+           return $ r @> 0
+
+-- | product of elements
+prodR :: Vector Double -> Double
+prodR x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodR vec x vec r "prodR"
+           return $ r @> 0
+
+-- | product of elements
+prodQ :: Vector (Complex Float) -> Complex Float
+prodQ x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodQ vec x vec r "prodQ"
+           return $ r @> 0
+
+-- | product of elements
+prodC :: Vector (Complex Double) -> Complex Double
+prodC x = unsafePerformIO $ do
+           r <- createVector 1
+           app2 c_prodC vec x vec r "prodC"
+           return $ r @> 0
+
+foreign import ccall unsafe "prodF" c_prodF :: TFF
+foreign import ccall unsafe "prodR" c_prodR :: TVV
+foreign import ccall unsafe "prodQ" c_prodQ :: TQVQV
+foreign import ccall unsafe "prodC" c_prodC :: TCVCV
+
+------------------------------------------------------------------
+
+toScalarAux fun code v = unsafePerformIO $ do
+    r <- createVector 1
+    app2 (fun (fromei code)) vec v vec r "toScalarAux"
+    return (r `at` 0)
+
+vectorMapAux fun code v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 (fun (fromei code)) vec v vec r "vectorMapAux"
+    return r
+
+vectorMapValAux fun code val v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    pval <- newArray [val]
+    app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux"
+    free pval
+    return r
+
+vectorZipAux fun code u v = unsafePerformIO $ do
+    r <- createVector (dim u)
+    app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux"
+    return r
+
+---------------------------------------------------------------------
+
+-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
+toScalarR :: FunCodeS -> Vector Double -> Double
+toScalarR oper =  toScalarAux c_toScalarR (fromei oper)
+
+foreign import ccall unsafe "toScalarR" c_toScalarR :: CInt -> TVV
+
+-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
+toScalarF :: FunCodeS -> Vector Float -> Float
+toScalarF oper =  toScalarAux c_toScalarF (fromei oper)
+
+foreign import ccall unsafe "toScalarF" c_toScalarF :: CInt -> TFF
+
+-- | obtains different functions of a vector: only norm1, norm2
+toScalarC :: FunCodeS -> Vector (Complex Double) -> Double
+toScalarC oper =  toScalarAux c_toScalarC (fromei oper)
+
+foreign import ccall unsafe "toScalarC" c_toScalarC :: CInt -> TCVV
+
+-- | obtains different functions of a vector: only norm1, norm2
+toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float
+toScalarQ oper =  toScalarAux c_toScalarQ (fromei oper)
+
+foreign import ccall unsafe "toScalarQ" c_toScalarQ :: CInt -> TQVF
+
+------------------------------------------------------------------
+
+-- | map of real vectors with given function
+vectorMapR :: FunCodeV -> Vector Double -> Vector Double
+vectorMapR = vectorMapAux c_vectorMapR
+
+foreign import ccall unsafe "mapR" c_vectorMapR :: CInt -> TVV
+
+-- | map of complex vectors with given function
+vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
+vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper)
+
+foreign import ccall unsafe "mapC" c_vectorMapC :: CInt -> TCVCV
+
+-- | map of real vectors with given function
+vectorMapF :: FunCodeV -> Vector Float -> Vector Float
+vectorMapF = vectorMapAux c_vectorMapF
+
+foreign import ccall unsafe "mapF" c_vectorMapF :: CInt -> TFF
+
+-- | map of real vectors with given function
+vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float)
+vectorMapQ = vectorMapAux c_vectorMapQ
+
+foreign import ccall unsafe "mapQ" c_vectorMapQ :: CInt -> TQVQV
+
+-------------------------------------------------------------------
+
+-- | map of real vectors with given function
+vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
+vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper)
+
+foreign import ccall unsafe "mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV
+
+-- | map of complex vectors with given function
+vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
+vectorMapValC = vectorMapValAux c_vectorMapValC
+
+foreign import ccall unsafe "mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV
+
+-- | map of real vectors with given function
+vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float
+vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)
+
+foreign import ccall unsafe "mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF
+
+-- | map of complex vectors with given function
+vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)
+vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)
+
+foreign import ccall unsafe "mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV
+
+-------------------------------------------------------------------
+
+-- | elementwise operation on real vectors
+vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
+vectorZipR = vectorZipAux c_vectorZipR
+
+foreign import ccall unsafe "zipR" c_vectorZipR :: CInt -> TVVV
+
+-- | elementwise operation on complex vectors
+vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
+vectorZipC = vectorZipAux c_vectorZipC
+
+foreign import ccall unsafe "zipC" c_vectorZipC :: CInt -> TCVCVCV
+
+-- | elementwise operation on real vectors
+vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float
+vectorZipF = vectorZipAux c_vectorZipF
+
+foreign import ccall unsafe "zipF" c_vectorZipF :: CInt -> TFFF
+
+-- | elementwise operation on complex vectors
+vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)
+vectorZipQ = vectorZipAux c_vectorZipQ
+
+foreign import ccall unsafe "zipQ" c_vectorZipQ :: CInt -> TQVQVQV
+
diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal
index 8496663..eaf4262 100644
--- a/packages/hmatrix/hmatrix.cabal
+++ b/packages/hmatrix/hmatrix.cabal
@@ -109,7 +109,8 @@ library
                         Numeric.LinearAlgebra.Util.Convolution
                         
 
-    C-sources:          src/Numeric/GSL/gsl-aux.c
+    C-sources:          src/Numeric/GSL/gsl-vector.c
+                        src/Numeric/GSL/gsl-aux.c
 
     cc-options:         -O4 -msse2 -Wall
 
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
index 29c8bb7..3591289 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs
@@ -13,7 +13,6 @@
 module Numeric.GSL.Vector (
     sumF, sumR, sumQ, sumC,
     prodF, prodR, prodQ, prodC,
-    dotF, dotR, dotQ, dotC,
     FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
     FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
     FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
@@ -148,39 +147,6 @@ foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV
 foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV
 foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV
 
--- | dot product
-dotF :: Vector Float -> Vector Float -> Float
-dotF x y = unsafePerformIO $ do
-           r <- createVector 1
-           app3 c_dotF vec x vec y vec r "dotF"
-           return $ r @> 0
-
--- | dot product
-dotR :: Vector Double -> Vector Double -> Double
-dotR x y = unsafePerformIO $ do
-           r <- createVector 1
-           app3 c_dotR vec x vec y vec r "dotR"
-           return $ r @> 0
-
--- | dot product
-dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float
-dotQ x y = unsafePerformIO $ do
-           r <- createVector 1
-           app3 c_dotQ vec x vec y vec r "dotQ"
-           return $ r @> 0
-
--- | dot product
-dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double
-dotC x y = unsafePerformIO $ do
-           r <- createVector 1
-           app3 c_dotC vec x vec y vec r "dotC"
-           return $ r @> 0
-
-foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF
-foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV
-foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV
-foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV
-
 ------------------------------------------------------------------
 
 toScalarAux fun code v = unsafePerformIO $ do
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
index 410d157..5da94ca 100644
--- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
@@ -104,602 +104,6 @@ void no_abort_on_error() {
 }
 
 
-int sumF(KFVEC(x),FVEC(r)) {
-    DEBUGMSG("sumF");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    float res = 0;
-    for (i = 0; i < xn; i++) res += xp[i];
-    rp[0] = res;
-    OK
-}
-    
-int sumR(KRVEC(x),RVEC(r)) {
-    DEBUGMSG("sumR");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    double res = 0;
-    for (i = 0; i < xn; i++) res += xp[i];
-    rp[0] = res;
-    OK
-}
-    
-int sumQ(KQVEC(x),QVEC(r)) {
-    DEBUGMSG("sumQ");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    gsl_complex_float res;
-    res.dat[0] = 0;
-    res.dat[1] = 0;
-    for (i = 0; i < xn; i++) {
-      res.dat[0] += xp[i].dat[0];
-      res.dat[1] += xp[i].dat[1];
-    }
-    rp[0] = res;
-    OK
-}
-    
-int sumC(KCVEC(x),CVEC(r)) {
-    DEBUGMSG("sumC");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    gsl_complex res;
-    res.dat[0] = 0;
-    res.dat[1] = 0;
-    for (i = 0; i < xn; i++)  {
-      res.dat[0] += xp[i].dat[0];
-      res.dat[1] += xp[i].dat[1];
-    }
-    rp[0] = res;
-    OK
-}
-
-int prodF(KFVEC(x),FVEC(r)) {
-    DEBUGMSG("prodF");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    float res = 1;
-    for (i = 0; i < xn; i++) res *= xp[i];
-    rp[0] = res;
-    OK
-}
-    
-int prodR(KRVEC(x),RVEC(r)) {
-    DEBUGMSG("prodR");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    double res = 1;
-    for (i = 0; i < xn; i++) res *= xp[i];
-    rp[0] = res;
-    OK
-}
-    
-int prodQ(KQVEC(x),QVEC(r)) {
-    DEBUGMSG("prodQ");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    gsl_complex_float res;
-    float temp;
-    res.dat[0] = 1;
-    res.dat[1] = 0;
-    for (i = 0; i < xn; i++) {
-      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
-      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
-      res.dat[0] = temp;
-    }
-    rp[0] = res;
-    OK
-}
-    
-int prodC(KCVEC(x),CVEC(r)) {
-    DEBUGMSG("prodC");
-    REQUIRES(rn==1,BAD_SIZE);
-    int i;
-    gsl_complex res;
-    double temp;
-    res.dat[0] = 1;
-    res.dat[1] = 0;
-    for (i = 0; i < xn; i++)  {
-      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
-      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
-      res.dat[0] = temp;
-    }
-    rp[0] = res;
-    OK
-}
-
-int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
-    DEBUGMSG("dotF");
-    REQUIRES(xn==yn,BAD_SIZE); 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("dotF");
-    KFVVIEW(x);
-    KFVVIEW(y);
-    gsl_blas_sdot(V(x),V(y),rp);
-    OK
-}
-    
-int dotR(KRVEC(x), KRVEC(y), RVEC(r)) {
-    DEBUGMSG("dotR");
-    REQUIRES(xn==yn,BAD_SIZE); 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("dotR");
-    KDVVIEW(x);
-    KDVVIEW(y);
-    gsl_blas_ddot(V(x),V(y),rp);
-    OK
-}
-    
-int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) {
-    DEBUGMSG("dotQ");
-    REQUIRES(xn==yn,BAD_SIZE); 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("dotQ");
-    KQVVIEW(x);
-    KQVVIEW(y);
-    gsl_blas_cdotu(V(x),V(y),rp);
-    OK
-}
-    
-int dotC(KCVEC(x), KCVEC(y), CVEC(r)) {
-    DEBUGMSG("dotC");
-    REQUIRES(xn==yn,BAD_SIZE); 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("dotC");
-    KCVVIEW(x);
-    KCVVIEW(y);
-    gsl_blas_zdotu(V(x),V(y),rp);
-    OK
-}
-    
-int toScalarR(int code, KRVEC(x), RVEC(r)) { 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("toScalarR");
-    KDVVIEW(x);
-    double res;
-    switch(code) {
-        case 0: { res = gsl_blas_dnrm2(V(x)); break; } 
-        case 1: { res = gsl_blas_dasum(V(x));  break; }
-        case 2: { res = gsl_vector_max_index(V(x));  break; }
-        case 3: { res = gsl_vector_max(V(x));  break; }
-        case 4: { res = gsl_vector_min_index(V(x)); break; }
-        case 5: { res = gsl_vector_min(V(x)); break; }
-        default: ERROR(BAD_CODE);
-    }
-    rp[0] = res;
-    OK
-}
-
-int toScalarF(int code, KFVEC(x), FVEC(r)) { 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("toScalarF");
-    KFVVIEW(x);
-    float res;
-    switch(code) {
-        case 0: { res = gsl_blas_snrm2(V(x)); break; } 
-        case 1: { res = gsl_blas_sasum(V(x));  break; }
-        case 2: { res = gsl_vector_float_max_index(V(x));  break; }
-        case 3: { res = gsl_vector_float_max(V(x));  break; }
-        case 4: { res = gsl_vector_float_min_index(V(x)); break; }
-        case 5: { res = gsl_vector_float_min(V(x)); break; }
-        default: ERROR(BAD_CODE);
-    }
-    rp[0] = res;
-    OK
-}
-
-
-int toScalarC(int code, KCVEC(x), RVEC(r)) { 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("toScalarC");
-    KCVVIEW(x);
-    double res;
-    switch(code) {
-        case 0: { res = gsl_blas_dznrm2(V(x)); break; } 
-        case 1: { res = gsl_blas_dzasum(V(x));  break; }
-        default: ERROR(BAD_CODE);
-    }
-    rp[0] = res;
-    OK
-}
-
-int toScalarQ(int code, KQVEC(x), FVEC(r)) { 
-    REQUIRES(rn==1,BAD_SIZE);
-    DEBUGMSG("toScalarQ");
-    KQVVIEW(x);
-    float res;
-    switch(code) {
-        case 0: { res = gsl_blas_scnrm2(V(x)); break; } 
-        case 1: { res = gsl_blas_scasum(V(x));  break; }
-        default: ERROR(BAD_CODE);
-    }
-    rp[0] = res;
-    OK
-}
-
-
-inline double sign(double x) {
-    if(x>0) {
-        return +1.0;
-    } else if (x<0) {
-        return -1.0;
-    } else {
-        return 0.0;
-    }
-}
-
-inline float float_sign(float x) {
-    if(x>0) {
-        return +1.0;
-    } else if (x<0) {
-        return -1.0;
-    } else {
-        return 0.0;
-    }
-}
-
-inline gsl_complex complex_abs(gsl_complex z) {
-    gsl_complex r;
-    r.dat[0] = gsl_complex_abs(z);
-    r.dat[1] = 0;
-    return r;
-}
-
-inline gsl_complex complex_signum(gsl_complex z) {
-    gsl_complex r;
-    double mag;
-    if (z.dat[0] == 0 && z.dat[1] == 0) {
-        r.dat[0] = 0;
-        r.dat[1] = 0;
-    } else {
-        mag = gsl_complex_abs(z);
-        r.dat[0] = z.dat[0]/mag;
-        r.dat[1] = z.dat[1]/mag;
-    }
-    return r;
-}
-
-#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
-#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
-int mapR(int code, KRVEC(x), RVEC(r)) {
-    int k;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapR");
-    switch (code) {
-        OP(0,sin)
-        OP(1,cos)
-        OP(2,tan)
-        OP(3,fabs)
-        OP(4,asin)
-        OP(5,acos)
-        OP(6,atan) /* atan2 mediante vectorZip */
-        OP(7,sinh)
-        OP(8,cosh)
-        OP(9,tanh)
-        OP(10,gsl_asinh)
-        OP(11,gsl_acosh)
-        OP(12,gsl_atanh)
-        OP(13,exp)
-        OP(14,log)
-        OP(15,sign)
-        OP(16,sqrt)
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapF(int code, KFVEC(x), FVEC(r)) {
-    int k;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapF");
-    switch (code) {
-        OP(0,sin)
-        OP(1,cos)
-        OP(2,tan)
-        OP(3,fabs)
-        OP(4,asin)
-        OP(5,acos)
-        OP(6,atan) /* atan2 mediante vectorZip */
-        OP(7,sinh)
-        OP(8,cosh)
-        OP(9,tanh)
-        OP(10,gsl_asinh)
-        OP(11,gsl_acosh)
-        OP(12,gsl_atanh)
-        OP(13,exp)
-        OP(14,log)
-        OP(15,sign)
-        OP(16,sqrt)
-        default: ERROR(BAD_CODE);
-    }
-}
-
-
-int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
-    int k;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapC");
-    switch (code) {
-        OP(0,gsl_complex_sin)
-        OP(1,gsl_complex_cos)
-        OP(2,gsl_complex_tan)
-        OP(3,complex_abs)
-        OP(4,gsl_complex_arcsin)
-        OP(5,gsl_complex_arccos)
-        OP(6,gsl_complex_arctan)
-        OP(7,gsl_complex_sinh)
-        OP(8,gsl_complex_cosh)
-        OP(9,gsl_complex_tanh)
-        OP(10,gsl_complex_arcsinh)
-        OP(11,gsl_complex_arccosh)
-        OP(12,gsl_complex_arctanh)
-        OP(13,gsl_complex_exp)
-        OP(14,gsl_complex_log)
-        OP(15,complex_signum)
-        OP(16,gsl_complex_sqrt)
-
-        // gsl_complex_arg
-        // gsl_complex_abs
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapC(int code, KCVEC(x), CVEC(r)) {
-    return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
-}
-
-
-gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
-{
-  gsl_complex c;
-  gsl_complex r;
-
-  gsl_complex_float float_r;
-
-  c.dat[0] = a.dat[0];
-  c.dat[1] = a.dat[1];
-
-  r = (*cf)(c);
-
-  float_r.dat[0] = r.dat[0];
-  float_r.dat[1] = r.dat[1];
-
-  return float_r;
-}
-
-gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), 
-                                        gsl_complex_float a,gsl_complex_float b)
-{
-  gsl_complex c1;
-  gsl_complex c2;
-  gsl_complex r;
-
-  gsl_complex_float float_r;
-
-  c1.dat[0] = a.dat[0];
-  c1.dat[1] = a.dat[1];
-
-  c2.dat[0] = b.dat[0];
-  c2.dat[1] = b.dat[1];
-
-  r = (*cf)(c1,c2);
-
-  float_r.dat[0] = r.dat[0];
-  float_r.dat[1] = r.dat[1];
-
-  return float_r;
-}
-
-#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
-#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
-int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
-    int k;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapQ");
-    switch (code) {
-        OPC(0,gsl_complex_sin)
-        OPC(1,gsl_complex_cos)
-        OPC(2,gsl_complex_tan)
-        OPC(3,complex_abs)
-        OPC(4,gsl_complex_arcsin)
-        OPC(5,gsl_complex_arccos)
-        OPC(6,gsl_complex_arctan)
-        OPC(7,gsl_complex_sinh)
-        OPC(8,gsl_complex_cosh)
-        OPC(9,gsl_complex_tanh)
-        OPC(10,gsl_complex_arcsinh)
-        OPC(11,gsl_complex_arccosh)
-        OPC(12,gsl_complex_arctanh)
-        OPC(13,gsl_complex_exp)
-        OPC(14,gsl_complex_log)
-        OPC(15,complex_signum)
-        OPC(16,gsl_complex_sqrt)
-
-        // gsl_complex_arg
-        // gsl_complex_abs
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapQ(int code, KQVEC(x), QVEC(r)) {
-    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
-}
-
-
-int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
-    int k;
-    double val = *pval;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapValR");
-    switch (code) {
-        OPV(0,val*xp[k])
-        OPV(1,val/xp[k])
-        OPV(2,val+xp[k])
-        OPV(3,val-xp[k])
-        OPV(4,pow(val,xp[k]))
-        OPV(5,pow(xp[k],val))
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
-    int k;
-    float val = *pval;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapValF");
-    switch (code) {
-        OPV(0,val*xp[k])
-        OPV(1,val/xp[k])
-        OPV(2,val+xp[k])
-        OPV(3,val-xp[k])
-        OPV(4,pow(val,xp[k]))
-        OPV(5,pow(xp[k],val))
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
-    int k;
-    gsl_complex val = *pval;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapValC");
-    switch (code) {
-        OPV(0,gsl_complex_mul(val,xp[k]))
-        OPV(1,gsl_complex_div(val,xp[k]))
-        OPV(2,gsl_complex_add(val,xp[k]))
-        OPV(3,gsl_complex_sub(val,xp[k]))
-        OPV(4,gsl_complex_pow(val,xp[k]))
-        OPV(5,gsl_complex_pow(xp[k],val))
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
-    return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
-}
-
-
-int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
-    int k;
-    gsl_complex_float val = *pval;
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("mapValQ");
-    switch (code) {
-        OPCA(0,gsl_complex_mul,val,xp[k])
-        OPCA(1,gsl_complex_div,val,xp[k])
-        OPCA(2,gsl_complex_add,val,xp[k])
-        OPCA(3,gsl_complex_sub,val,xp[k])
-        OPCA(4,gsl_complex_pow,val,xp[k])
-        OPCA(5,gsl_complex_pow,xp[k],val)
-        default: ERROR(BAD_CODE);
-    }
-}
-
-int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
-    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
-}
-
-
-#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
-#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
-int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
-    REQUIRES(an == bn && an == rn, BAD_SIZE);
-    int k;
-    switch(code) {
-        OPZE(4,"zipR Pow",pow)
-        OPZE(5,"zipR ATan2",atan2)
-    }
-    KDVVIEW(a);
-    KDVVIEW(b);
-    DVVIEW(r);
-    gsl_vector_memcpy(V(r),V(a));
-    int res;
-    switch(code) {
-        OPZV(0,"zipR Add",gsl_vector_add)
-        OPZV(1,"zipR Sub",gsl_vector_sub)
-        OPZV(2,"zipR Mul",gsl_vector_mul)
-        OPZV(3,"zipR Div",gsl_vector_div)
-        default: ERROR(BAD_CODE);
-    }
-}
-
-
-int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
-    REQUIRES(an == bn && an == rn, BAD_SIZE);
-    int k;
-    switch(code) {
-        OPZE(4,"zipF Pow",pow)
-        OPZE(5,"zipF ATan2",atan2)
-    }
-    KFVVIEW(a);
-    KFVVIEW(b);
-    FVVIEW(r);
-    gsl_vector_float_memcpy(V(r),V(a));
-    int res;
-    switch(code) {
-        OPZV(0,"zipF Add",gsl_vector_float_add)
-        OPZV(1,"zipF Sub",gsl_vector_float_sub)
-        OPZV(2,"zipF Mul",gsl_vector_float_mul)
-        OPZV(3,"zipF Div",gsl_vector_float_div)
-        default: ERROR(BAD_CODE);
-    }
-}
-
-
-int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
-    REQUIRES(an == bn && an == rn, BAD_SIZE);
-    int k;
-    switch(code) {
-        OPZE(0,"zipC Add",gsl_complex_add)
-        OPZE(1,"zipC Sub",gsl_complex_sub)
-        OPZE(2,"zipC Mul",gsl_complex_mul)
-        OPZE(3,"zipC Div",gsl_complex_div)
-        OPZE(4,"zipC Pow",gsl_complex_pow)
-        //OPZE(5,"zipR ATan2",atan2)
-    }
-    //KCVVIEW(a);
-    //KCVVIEW(b);
-    //CVVIEW(r);
-    //gsl_vector_memcpy(V(r),V(a));
-    //int res;
-    switch(code) {
-        default: ERROR(BAD_CODE);
-    }
-}
-
-
-int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
-    return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
-}
-
-
-#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
-int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
-    REQUIRES(an == bn && an == rn, BAD_SIZE);
-    int k;
-    switch(code) {
-        OPCZE(0,"zipQ Add",gsl_complex_add)
-        OPCZE(1,"zipQ Sub",gsl_complex_sub)
-        OPCZE(2,"zipQ Mul",gsl_complex_mul)
-        OPCZE(3,"zipQ Div",gsl_complex_div)
-        OPCZE(4,"zipQ Pow",gsl_complex_pow)
-        //OPZE(5,"zipR ATan2",atan2)
-    }
-    //KCVVIEW(a);
-    //KCVVIEW(b);
-    //CVVIEW(r);
-    //gsl_vector_memcpy(V(r),V(a));
-    //int res;
-    switch(code) {
-        default: ERROR(BAD_CODE);
-    }
-}
-
-
-int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
-    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
-}
-
-
 
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-vector.c b/packages/hmatrix/src/Numeric/GSL/gsl-vector.c
new file mode 100644
index 0000000..40a086a
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-vector.c
@@ -0,0 +1,642 @@
+#include <gsl/gsl_complex.h>
+
+#define RVEC(A) int A##n, double*A##p
+#define RMAT(A) int A##r, int A##c, double* A##p
+#define KRVEC(A) int A##n, const double*A##p
+#define KRMAT(A) int A##r, int A##c, const double* A##p
+
+#define CVEC(A) int A##n, gsl_complex*A##p
+#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
+#define KCVEC(A) int A##n, const gsl_complex*A##p
+#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
+
+#define FVEC(A) int A##n, float*A##p
+#define FMAT(A) int A##r, int A##c, float* A##p
+#define KFVEC(A) int A##n, const float*A##p
+#define KFMAT(A) int A##r, int A##c, const float* A##p
+
+#define QVEC(A) int A##n, gsl_complex_float*A##p
+#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
+#define KQVEC(A) int A##n, const gsl_complex_float*A##p
+#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.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 OK return 0;
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
+#else
+#define DEBUGMSG(M)
+#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((double*)A##p,A##n)
+#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)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((double*)A##p,A##n)
+#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
+
+#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
+#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
+#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
+#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
+#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
+#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
+#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
+#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)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 GQVEC(A) int A##n, gsl_complex_float*A##p
+#define KGQVEC(A) int A##n, const gsl_complex_float*A##p
+
+#define BAD_SIZE 2000
+#define BAD_CODE 2001
+#define MEM      2002
+#define BAD_FILE 2003
+
+
+int sumF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("sumF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int sumR(KRVEC(x),RVEC(r)) {
+    DEBUGMSG("sumR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 0;
+    for (i = 0; i < xn; i++) res += xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int sumQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("sumQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex_float res;
+    res.dat[0] = 0;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++) {
+      res.dat[0] += xp[i].dat[0];
+      res.dat[1] += xp[i].dat[1];
+    }
+    rp[0] = res;
+    OK
+}
+    
+int sumC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("sumC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex res;
+    res.dat[0] = 0;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++)  {
+      res.dat[0] += xp[i].dat[0];
+      res.dat[1] += xp[i].dat[1];
+    }
+    rp[0] = res;
+    OK
+}
+
+int prodF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("prodF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int prodR(KRVEC(x),RVEC(r)) {
+    DEBUGMSG("prodR");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    double res = 1;
+    for (i = 0; i < xn; i++) res *= xp[i];
+    rp[0] = res;
+    OK
+}
+    
+int prodQ(KQVEC(x),QVEC(r)) {
+    DEBUGMSG("prodQ");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex_float res;
+    float temp;
+    res.dat[0] = 1;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++) {
+      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
+      res.dat[0] = temp;
+    }
+    rp[0] = res;
+    OK
+}
+    
+int prodC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("prodC");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    gsl_complex res;
+    double temp;
+    res.dat[0] = 1;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++)  {
+      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
+      res.dat[0] = temp;
+    }
+    rp[0] = res;
+    OK
+}
+
+    
+int toScalarR(int code, KRVEC(x), RVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarR");
+    KDVVIEW(x);
+    double res;
+    switch(code) {
+        case 0: { res = gsl_blas_dnrm2(V(x)); break; } 
+        case 1: { res = gsl_blas_dasum(V(x));  break; }
+        case 2: { res = gsl_vector_max_index(V(x));  break; }
+        case 3: { res = gsl_vector_max(V(x));  break; }
+        case 4: { res = gsl_vector_min_index(V(x)); break; }
+        case 5: { res = gsl_vector_min(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+int toScalarF(int code, KFVEC(x), FVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarF");
+    KFVVIEW(x);
+    float res;
+    switch(code) {
+        case 0: { res = gsl_blas_snrm2(V(x)); break; } 
+        case 1: { res = gsl_blas_sasum(V(x));  break; }
+        case 2: { res = gsl_vector_float_max_index(V(x));  break; }
+        case 3: { res = gsl_vector_float_max(V(x));  break; }
+        case 4: { res = gsl_vector_float_min_index(V(x)); break; }
+        case 5: { res = gsl_vector_float_min(V(x)); break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+int toScalarC(int code, KCVEC(x), RVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarC");
+    KCVVIEW(x);
+    double res;
+    switch(code) {
+        case 0: { res = gsl_blas_dznrm2(V(x)); break; } 
+        case 1: { res = gsl_blas_dzasum(V(x));  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+int toScalarQ(int code, KQVEC(x), FVEC(r)) { 
+    REQUIRES(rn==1,BAD_SIZE);
+    DEBUGMSG("toScalarQ");
+    KQVVIEW(x);
+    float res;
+    switch(code) {
+        case 0: { res = gsl_blas_scnrm2(V(x)); break; } 
+        case 1: { res = gsl_blas_scasum(V(x));  break; }
+        default: ERROR(BAD_CODE);
+    }
+    rp[0] = res;
+    OK
+}
+
+
+inline double sign(double x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+inline float float_sign(float x) {
+    if(x>0) {
+        return +1.0;
+    } else if (x<0) {
+        return -1.0;
+    } else {
+        return 0.0;
+    }
+}
+
+inline gsl_complex complex_abs(gsl_complex z) {
+    gsl_complex r;
+    r.dat[0] = gsl_complex_abs(z);
+    r.dat[1] = 0;
+    return r;
+}
+
+inline gsl_complex complex_signum(gsl_complex z) {
+    gsl_complex r;
+    double mag;
+    if (z.dat[0] == 0 && z.dat[1] == 0) {
+        r.dat[0] = 0;
+        r.dat[1] = 0;
+    } else {
+        mag = gsl_complex_abs(z);
+        r.dat[0] = z.dat[0]/mag;
+        r.dat[1] = z.dat[1]/mag;
+    }
+    return r;
+}
+
+#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
+#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
+int mapR(int code, KRVEC(x), RVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapR");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan) /* atan2 mediante vectorZip */
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+        OP(10,gsl_asinh)
+        OP(11,gsl_acosh)
+        OP(12,gsl_atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapF(int code, KFVEC(x), FVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapF");
+    switch (code) {
+        OP(0,sin)
+        OP(1,cos)
+        OP(2,tan)
+        OP(3,fabs)
+        OP(4,asin)
+        OP(5,acos)
+        OP(6,atan) /* atan2 mediante vectorZip */
+        OP(7,sinh)
+        OP(8,cosh)
+        OP(9,tanh)
+        OP(10,gsl_asinh)
+        OP(11,gsl_acosh)
+        OP(12,gsl_atanh)
+        OP(13,exp)
+        OP(14,log)
+        OP(15,sign)
+        OP(16,sqrt)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapC");
+    switch (code) {
+        OP(0,gsl_complex_sin)
+        OP(1,gsl_complex_cos)
+        OP(2,gsl_complex_tan)
+        OP(3,complex_abs)
+        OP(4,gsl_complex_arcsin)
+        OP(5,gsl_complex_arccos)
+        OP(6,gsl_complex_arctan)
+        OP(7,gsl_complex_sinh)
+        OP(8,gsl_complex_cosh)
+        OP(9,gsl_complex_tanh)
+        OP(10,gsl_complex_arcsinh)
+        OP(11,gsl_complex_arccosh)
+        OP(12,gsl_complex_arctanh)
+        OP(13,gsl_complex_exp)
+        OP(14,gsl_complex_log)
+        OP(15,complex_signum)
+        OP(16,gsl_complex_sqrt)
+
+        // gsl_complex_arg
+        // gsl_complex_abs
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapC(int code, KCVEC(x), CVEC(r)) {
+    return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
+}
+
+
+gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
+{
+  gsl_complex c;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c.dat[0] = a.dat[0];
+  c.dat[1] = a.dat[1];
+
+  r = (*cf)(c);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), 
+                                        gsl_complex_float a,gsl_complex_float b)
+{
+  gsl_complex c1;
+  gsl_complex c2;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c1.dat[0] = a.dat[0];
+  c1.dat[1] = a.dat[1];
+
+  c2.dat[0] = b.dat[0];
+  c2.dat[1] = b.dat[1];
+
+  r = (*cf)(c1,c2);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
+#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
+int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
+    int k;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapQ");
+    switch (code) {
+        OPC(0,gsl_complex_sin)
+        OPC(1,gsl_complex_cos)
+        OPC(2,gsl_complex_tan)
+        OPC(3,complex_abs)
+        OPC(4,gsl_complex_arcsin)
+        OPC(5,gsl_complex_arccos)
+        OPC(6,gsl_complex_arctan)
+        OPC(7,gsl_complex_sinh)
+        OPC(8,gsl_complex_cosh)
+        OPC(9,gsl_complex_tanh)
+        OPC(10,gsl_complex_arcsinh)
+        OPC(11,gsl_complex_arccosh)
+        OPC(12,gsl_complex_arctanh)
+        OPC(13,gsl_complex_exp)
+        OPC(14,gsl_complex_log)
+        OPC(15,complex_signum)
+        OPC(16,gsl_complex_sqrt)
+
+        // gsl_complex_arg
+        // gsl_complex_abs
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapQ(int code, KQVEC(x), QVEC(r)) {
+    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
+int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
+    int k;
+    double val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValR");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
+    int k;
+    float val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValF");
+    switch (code) {
+        OPV(0,val*xp[k])
+        OPV(1,val/xp[k])
+        OPV(2,val+xp[k])
+        OPV(3,val-xp[k])
+        OPV(4,pow(val,xp[k]))
+        OPV(5,pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
+    int k;
+    gsl_complex val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValC");
+    switch (code) {
+        OPV(0,gsl_complex_mul(val,xp[k]))
+        OPV(1,gsl_complex_div(val,xp[k]))
+        OPV(2,gsl_complex_add(val,xp[k]))
+        OPV(3,gsl_complex_sub(val,xp[k]))
+        OPV(4,gsl_complex_pow(val,xp[k]))
+        OPV(5,gsl_complex_pow(xp[k],val))
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
+    return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
+}
+
+
+int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
+    int k;
+    gsl_complex_float val = *pval;
+    REQUIRES(xn == rn,BAD_SIZE);
+    DEBUGMSG("mapValQ");
+    switch (code) {
+        OPCA(0,gsl_complex_mul,val,xp[k])
+        OPCA(1,gsl_complex_div,val,xp[k])
+        OPCA(2,gsl_complex_add,val,xp[k])
+        OPCA(3,gsl_complex_sub,val,xp[k])
+        OPCA(4,gsl_complex_pow,val,xp[k])
+        OPCA(5,gsl_complex_pow,xp[k],val)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
+    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
+#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
+#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
+int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZE(4,"zipR Pow",pow)
+        OPZE(5,"zipR ATan2",atan2)
+    }
+    KDVVIEW(a);
+    KDVVIEW(b);
+    DVVIEW(r);
+    gsl_vector_memcpy(V(r),V(a));
+    int res;
+    switch(code) {
+        OPZV(0,"zipR Add",gsl_vector_add)
+        OPZV(1,"zipR Sub",gsl_vector_sub)
+        OPZV(2,"zipR Mul",gsl_vector_mul)
+        OPZV(3,"zipR Div",gsl_vector_div)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZE(4,"zipF Pow",pow)
+        OPZE(5,"zipF ATan2",atan2)
+    }
+    KFVVIEW(a);
+    KFVVIEW(b);
+    FVVIEW(r);
+    gsl_vector_float_memcpy(V(r),V(a));
+    int res;
+    switch(code) {
+        OPZV(0,"zipF Add",gsl_vector_float_add)
+        OPZV(1,"zipF Sub",gsl_vector_float_sub)
+        OPZV(2,"zipF Mul",gsl_vector_float_mul)
+        OPZV(3,"zipF Div",gsl_vector_float_div)
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPZE(0,"zipC Add",gsl_complex_add)
+        OPZE(1,"zipC Sub",gsl_complex_sub)
+        OPZE(2,"zipC Mul",gsl_complex_mul)
+        OPZE(3,"zipC Div",gsl_complex_div)
+        OPZE(4,"zipC Pow",gsl_complex_pow)
+        //OPZE(5,"zipR ATan2",atan2)
+    }
+    //KCVVIEW(a);
+    //KCVVIEW(b);
+    //CVVIEW(r);
+    //gsl_vector_memcpy(V(r),V(a));
+    //int res;
+    switch(code) {
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
+    return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
+}
+
+
+#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
+int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPCZE(0,"zipQ Add",gsl_complex_add)
+        OPCZE(1,"zipQ Sub",gsl_complex_sub)
+        OPCZE(2,"zipQ Mul",gsl_complex_mul)
+        OPCZE(3,"zipQ Div",gsl_complex_div)
+        OPCZE(4,"zipQ Pow",gsl_complex_pow)
+        //OPZE(5,"zipR ATan2",atan2)
+    }
+    //KCVVIEW(a);
+    //KCVVIEW(b);
+    //CVVIEW(r);
+    //gsl_vector_memcpy(V(r),V(a));
+    //int res;
+    switch(code) {
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
+    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
+}
+