3 files changed, 951 insertions, 0 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
+