added Vectors typeclass and refactored

2 files changed, 168 insertions, 4 deletions

diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index d09323b..97a0f9c 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -14,6 +14,8 @@
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.Vector (
+    sumF, sumR, sumQ, sumC,
+    dotF, dotR, dotQ, dotC,
     FunCodeS(..), toScalarR, toScalarF,
     FunCodeV(..), vectorMapR, vectorMapC, vectorMapF,
     FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF,
@@ -76,6 +78,74 @@ data FunCodeS = Norm2
 
 ------------------------------------------------------------------
 
+-- | 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 safe "gsl-aux.h sumF" c_sumF :: TFF
+foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV
+foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV
+foreign import ccall safe "gsl-aux.h sumC" c_sumC :: 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 safe "gsl-aux.h dotF" c_dotF :: TFFF
+foreign import ccall safe "gsl-aux.h dotR" c_dotR :: TVVV
+foreign import ccall safe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV
+foreign import ccall safe "gsl-aux.h dotC" c_dotC :: TCVCVCV
+
+------------------------------------------------------------------
+
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
     app2 (fun (fromei code)) vec v vec r "toScalarAux"
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 6bb16f0..fe33766 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -76,12 +76,12 @@
 
 #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_float_complex_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
-#define QMVIEW(A) gsl_matrix_float_complex_view A = gsl_matrix_float_complex_view_array((float*)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_float_complex_const_view A = gsl_vector_float_complex_const_view_array((float*)A##p,A##n)
-#define KQMVIEW(A) gsl_matrix_float_complex_const_view A = gsl_matrix_float_complex_const_view_array((float*)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)
@@ -103,6 +103,100 @@ 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 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");