add vector Product to Vectors

3 files changed, 139 insertions, 50 deletions

diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 14ba0ff..e8df27c 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -15,6 +15,7 @@
 
 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,
@@ -111,6 +112,39 @@ 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
 
+-- | 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 safe "gsl-aux.h prodF" c_prodF :: TFF
+foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV
+foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV
+foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV
+
 -- | dot product
 dotF :: Vector Float -> Vector Float -> Float
 dotF x y = unsafePerformIO $ do
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 689989d..3e9e2cb 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -153,6 +153,56 @@ int sumC(KCVEC(x),CVEC(r)) {
     OK
 }
 
+int prodF(KFVEC(x),FVEC(r)) {
+    DEBUGMSG("prodF");
+    REQUIRES(rn==1,BAD_SIZE);
+    int i;
+    float res = 0;
+    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 = 0;
+    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;
+    res.dat[0] = 0;
+    res.dat[1] = 0;
+    for (i = 0; i < xn; i++) {
+      res.dat[0] = 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];
+    }
+    rp[0] = res;
+    OK
+}
+    
+int prodC(KCVEC(x),CVEC(r)) {
+    DEBUGMSG("prodC");
+    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] = 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];
+    }
+    rp[0] = res;
+    OK
+}
+
 int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
     DEBUGMSG("dotF");
     REQUIRES(xn==yn,BAD_SIZE); 
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 4d7f608..778b976 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -39,65 +39,70 @@ import Control.Monad(ap)
 -- | basic Vector functions
 class Num e => Vectors a e where
     -- the C functions sumX are twice as fast as using foldVector
-    vectorSum :: a e -> e
-    absSum    :: a e -> e
-    vectorMin :: a e -> e
-    vectorMax :: a e -> e
-    minIdx    :: a e -> Int
-    maxIdx    :: a e -> Int
-    dot       :: a e -> a e -> e
-    norm1     :: a e -> e
-    norm2     :: a e -> e
-    normInf   :: a e -> e
+    vectorSum  :: a e -> e
+    vectorProd :: a e -> e
+    absSum     :: a e -> e
+    vectorMin  :: a e -> e
+    vectorMax  :: a e -> e
+    minIdx     :: a e -> Int
+    maxIdx     :: a e -> Int
+    dot        :: a e -> a e -> e
+    norm1      :: a e -> e
+    norm2      :: a e -> e
+    normInf    :: a e -> e
 
 
 instance Vectors Vector Float where
-    vectorSum = sumF
-    norm2     = toScalarF Norm2
-    absSum    = toScalarF AbsSum
-    vectorMin = toScalarF Min
-    vectorMax = toScalarF Max
-    minIdx    = round . toScalarF MinIdx
-    maxIdx    = round . toScalarF MaxIdx
-    dot       = dotF
-    norm1     = toScalarF AbsSum
-    normInf   = vectorMax . vectorMapF Abs
+    vectorSum  = sumF
+    vectorProd = prodF
+    norm2      = toScalarF Norm2
+    absSum     = toScalarF AbsSum
+    vectorMin  = toScalarF Min
+    vectorMax  = toScalarF Max
+    minIdx     = round . toScalarF MinIdx
+    maxIdx     = round . toScalarF MaxIdx
+    dot        = dotF
+    norm1      = toScalarF AbsSum
+    normInf    = vectorMax . vectorMapF Abs
 
 instance Vectors Vector Double where
-    vectorSum = sumR
-    norm2     = toScalarR Norm2
-    absSum    = toScalarR AbsSum
-    vectorMin = toScalarR Min
-    vectorMax = toScalarR Max
-    minIdx    = round . toScalarR MinIdx
-    maxIdx    = round . toScalarR MaxIdx
-    dot       = dotR
-    norm1     = toScalarR AbsSum
-    normInf   = vectorMax . vectorMapR Abs
+    vectorSum  = sumR
+    vectorProd = prodR
+    norm2      = toScalarR Norm2
+    absSum     = toScalarR AbsSum
+    vectorMin  = toScalarR Min
+    vectorMax  = toScalarR Max
+    minIdx     = round . toScalarR MinIdx
+    maxIdx     = round . toScalarR MaxIdx
+    dot        = dotR
+    norm1      = toScalarR AbsSum
+    normInf    = vectorMax . vectorMapR Abs
 
 instance Vectors Vector (Complex Float) where
-    vectorSum = sumQ
-    norm2     = (:+ 0) . toScalarQ Norm2
-    absSum    = (:+ 0) . toScalarQ AbsSum
-    vectorMin = ap (@>) minIdx
-    vectorMax = ap (@>) maxIdx
-    minIdx    = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIdx    = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    dot       = dotQ
-    norm1     = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
-    normInf   = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
+    vectorSum  = sumQ
+    vectorProd = prodQ
+    norm2      = (:+ 0) . toScalarQ Norm2
+    absSum     = (:+ 0) . toScalarQ AbsSum
+    vectorMin  = ap (@>) minIdx
+    vectorMax  = ap (@>) maxIdx
+    minIdx     = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIdx     = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    dot        = dotQ
+    norm1      = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
+    normInf    = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
 
 instance Vectors Vector (Complex Double) where
-    vectorSum = sumC
-    norm2     = (:+ 0) . toScalarC Norm2
-    absSum    = (:+ 0) . toScalarC AbsSum
-    vectorMin = ap (@>) minIdx
-    vectorMax = ap (@>) maxIdx
-    minIdx    = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIdx    = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    dot       = dotC
-    norm1     = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
-    normInf   = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
+    vectorSum  = sumC
+    vectorProd = prodC
+    norm2      = (:+ 0) . toScalarC Norm2
+    absSum     = (:+ 0) . toScalarC AbsSum
+    vectorMin  = ap (@>) minIdx
+    vectorMax  = ap (@>) maxIdx
+    minIdx     = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIdx     = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    dot        = dotC
+    norm1      = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
+    normInf    = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
 
 ----------------------------------------------------