From 0f9575462eb37a7c9985583ca33ae315f6e6431d Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 14 May 2014 19:46:44 +0200 Subject: vectorized operations in base --- packages/hmatrix/hmatrix.cabal | 3 +- packages/hmatrix/src/Numeric/GSL/Vector.hs | 34 -- packages/hmatrix/src/Numeric/GSL/gsl-aux.c | 596 ------------------------ packages/hmatrix/src/Numeric/GSL/gsl-vector.c | 642 ++++++++++++++++++++++++++ 4 files changed, 644 insertions(+), 631 deletions(-) create mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-vector.c (limited to 'packages/hmatrix') 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 + +#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 +#include +#include +#include +#include +#include + +#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