add-vector-float

11 files changed, 251 insertions, 9 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 003e8ee..7b3b305 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -257,6 +257,10 @@ class (Storable a, Floating a) => Element a where
     constantD  :: a -> Int -> Vector a
     constantD = constant'
 
+instance Element Float where
+    transdata  = transdataAux ctransF
+    constantD  = constantAux cconstantF
+
 instance Element Double where
     transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
@@ -308,6 +312,7 @@ transdataAux fun c1 d c2 =
         r2 = dim d `div` c2
         noneed = r1 == 1 || c1 == 1
 
+foreign import ccall "transF" ctransF :: TFMFM
 foreign import ccall "transR" ctransR :: TMM
 foreign import ccall "transC" ctransC :: TCMCM
 ----------------------------------------------------------------------
@@ -329,6 +334,10 @@ constantAux fun x n = unsafePerformIO $ do
     free px
     return v
 
+constantF :: Float -> Int -> Vector Float
+constantF = constantAux cconstantF
+foreign import ccall "constantF" cconstantF :: Ptr Float -> TF
+
 constantR :: Double -> Int -> Vector Double
 constantR = constantAux cconstantR
 foreign import ccall "constantR" cconstantR :: Ptr Double -> TV
@@ -368,15 +377,15 @@ subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
 --------------------------------------------------------------------------
 
 -- | obtains the complex conjugate of a complex vector
-conjV :: Vector (Complex Double) -> Vector (Complex Double)
+conjV :: (Storable a, RealFloat a) => Vector (Complex a) -> Vector (Complex a)
 conjV = mapVector conjugate
 
 -- | creates a complex vector from vectors with real and imaginary parts
-toComplexV :: (Vector Double, Vector Double) ->  Vector (Complex Double)
+toComplexV :: Element a => (Vector a, Vector a) ->  Vector (Complex a)
 toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
 
 -- | the inverse of 'toComplex'
-fromComplexV :: Vector (Complex Double) -> (Vector Double, Vector Double)
+fromComplexV :: Element a => Vector (Complex a) -> (Vector a, Vector a)
 fromComplexV z = (r,i) where
     [r,i] = toColumns $ reshape 2 $ asReal z
 
diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs
index d3ce121..4d984e3 100644
--- a/lib/Data/Packed/Internal/Signatures.hs
+++ b/lib/Data/Packed/Internal/Signatures.hs
@@ -18,11 +18,17 @@ import Foreign
 import Data.Complex
 import Foreign.C.Types
 
+type PF = Ptr Float                             --
 type PD = Ptr Double                            --
 type PC = Ptr (Complex Double)                  --
+type TF = CInt -> PF -> IO CInt                 --
+type TFF = CInt -> PF -> TF                     --
+type TFFF = CInt -> PF -> TFF                   --
 type TV = CInt -> PD -> IO CInt                 --
 type TVV = CInt -> PD -> TV                     --
 type TVVV = CInt -> PD -> TVV                   --
+type TFM = CInt -> CInt -> PF -> IO CInt        --
+type TFMFM =  CInt -> CInt -> PF -> TFM         --
 type TM = CInt -> CInt -> PD -> IO CInt         --
 type TMM =  CInt -> CInt -> PD -> TM            --
 type TVMM = CInt -> PD -> TMM                   --
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index a6868d9..d97a86d 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -264,13 +264,13 @@ takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (
 ---------------------------------------------------------------
 
 -- | transforms a complex vector into a real vector with alternating real and imaginary parts 
-asReal :: Vector (Complex Double) -> Vector Double
+asReal :: Vector (Complex a) -> Vector a
 --asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr =  castForeignPtr (fptr v) }
 asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
     where (fp,i,n) = unsafeToForeignPtr v
 
 -- | transforms a real vector into a complex vector with alternating real and imaginary parts
-asComplex :: Vector Double -> Vector (Complex Double)
+asComplex :: Vector a -> Vector (Complex a)
 --asComplex v = V { ioff = ioff v `div` 2, idim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }
 asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
     where (fp,i,n) = unsafeToForeignPtr v
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 0c21b97..c6d8a90 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -455,6 +455,14 @@ class (Element e) => Container c e where
     real        :: c Double -> c e
     complex     :: c e -> c (Complex Double)
 
+instance Container Vector Float where
+    toComplex = toComplexV
+    fromComplex = fromComplexV
+    comp v = toComplex (v,constant 0 (dim v))
+    conj = conjV
+    real = mapVector realToFrac
+    complex = (mapVector (\(r :+ i) -> (realToFrac r :+ realToFrac i))) . comp
+
 instance Container Vector Double where
     toComplex = toComplexV
     fromComplex = fromComplexV
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index f6b3fc6..66aa71d 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -19,6 +19,7 @@ module Data.Packed.Vector (
     subVector, takesV, join,
     constant, linspace,
     vecdisp,
+    vectorFMax, vectorFMin, vectorFMaxIndex, vectorFMinIndex,
     vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,
     mapVector, zipVector,
     fscanfVector, fprintfVector, freadVector, fwriteVector,
@@ -63,6 +64,18 @@ linspace n (a,b) = add a $ scale s  $ fromList [0 .. fromIntegral n-1]
           add   = vectorMapValR AddConstant
           s = (b-a)/fromIntegral (n-1)
 
+vectorFMax :: Vector Float -> Float
+vectorFMax = toScalarF Max
+
+vectorFMin :: Vector Float -> Float
+vectorFMin = toScalarF Min
+
+vectorFMaxIndex :: Vector Float -> Int
+vectorFMaxIndex = round . toScalarF MaxIdx
+
+vectorFMinIndex :: Vector Float -> Int
+vectorFMinIndex = round . toScalarF MinIdx
+
 vectorMax :: Vector Double -> Double
 vectorMax = toScalarR Max
 
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 28c3b06..d09323b 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -14,10 +14,10 @@
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.Vector (
-    FunCodeS(..), toScalarR,
-    FunCodeV(..), vectorMapR, vectorMapC,
-    FunCodeSV(..), vectorMapValR, vectorMapValC,
-    FunCodeVV(..), vectorZipR, vectorZipC,
+    FunCodeS(..), toScalarR, toScalarF,
+    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF,
+    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF,
+    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF,
     RandDist(..), randomVector
 ) where
 
@@ -106,6 +106,12 @@ toScalarR oper =  toScalarAux c_toScalarR (fromei oper)
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF
+
 ------------------------------------------------------------------
 
 -- | map of real vectors with given function
@@ -120,6 +126,12 @@ vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper)
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF
+
 -------------------------------------------------------------------
 
 -- | map of real vectors with given function
@@ -134,6 +146,12 @@ vectorMapValC = vectorMapValAux c_vectorMapValC
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF
+
 -------------------------------------------------------------------
 
 -- | elementwise operation on real vectors
@@ -148,6 +166,12 @@ vectorZipC = vectorZipAux c_vectorZipC
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF
+
 -----------------------------------------------------------------------
 
 data RandDist = Uniform  -- ^ uniform distribution in [0,1)
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index f5a95f4..6bb16f0 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -10,6 +10,16 @@
 #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>
@@ -64,12 +74,24 @@
 #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_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 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 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
@@ -99,6 +121,24 @@ int toScalarR(int code, KRVEC(x), RVEC(r)) {
     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
+}
+
 
 inline double sign(double x) {
     if(x>0) {
@@ -110,6 +150,16 @@ inline double sign(double x) {
     }
 }
 
+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);
@@ -159,6 +209,32 @@ int mapR(int code, KRVEC(x), RVEC(r)) {
     }
 }
 
+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;
@@ -210,6 +286,22 @@ int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
     }
 }
 
+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;
@@ -255,6 +347,28 @@ int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
 }
 
 
+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;
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index 1992db0..bba89c8 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -94,6 +94,14 @@ instance Linear Vector a => Eq (Vector a) where
 
 #endif
 
+instance Num (Vector Float) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapF Sign
+    abs = vectorMapF Abs
+    fromInteger = fromList . return . fromInteger
+
 instance Num (Vector Double) where
     (+) = adaptScalar addConstant add (flip addConstant)
     negate = scale (-1)
@@ -138,6 +146,27 @@ instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional
 
 ---------------------------------------------------------
 
+instance Floating (Vector Float) where
+    sin   = vectorMapF Sin
+    cos   = vectorMapF Cos
+    tan   = vectorMapF Tan
+    asin  = vectorMapF ASin
+    acos  = vectorMapF ACos
+    atan  = vectorMapF ATan
+    sinh  = vectorMapF Sinh
+    cosh  = vectorMapF Cosh
+    tanh  = vectorMapF Tanh
+    asinh = vectorMapF ASinh
+    acosh = vectorMapF ACosh
+    atanh = vectorMapF ATanh
+    exp   = vectorMapF Exp
+    log   = vectorMapF Log
+    sqrt  = vectorMapF Sqrt
+    (**)  = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))
+    pi    = fromList [pi]
+
+-------------------------------------------------------------
+
 instance Floating (Vector Double) where
     sin   = vectorMapR Sin
     cos   = vectorMapR Cos
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index fd840e3..b9c2572 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -1039,6 +1039,18 @@ int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
 
 //////////////////// transpose /////////////////////////
 
+int transF(KFMAT(x),FMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transF");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        tp[i*tc+j] = xp[j*xc+i];
+        }
+    }
+    OK
+}
+
 int transR(KDMAT(x),DMAT(t)) {
     REQUIRES(xr==tc && xc==tr,BAD_SIZE);
     DEBUGMSG("transR");
@@ -1065,6 +1077,16 @@ int transC(KCMAT(x),CMAT(t)) {
 
 //////////////////// constant /////////////////////////
 
+int constantF(float * pval, FVEC(r)) {
+    DEBUGMSG("constantF")
+    int k;
+    double val = *pval;
+    for(k=0;k<rn;k++) {
+        rp[k]=val;
+    }
+    OK
+}
+
 int constantR(double * pval, DVEC(r)) {
     DEBUGMSG("constantR")
     int k;
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index adf096e..415a6ab 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -40,13 +40,17 @@ typedef short ftnlen;
 
 /********************************************************/
 
+#define FVEC(A) int A##n, float*A##p
 #define DVEC(A) int A##n, double*A##p
 #define CVEC(A) int A##n, double*A##p
+#define FMAT(A) int A##r, int A##c, float* A##p
 #define DMAT(A) int A##r, int A##c, double* A##p
 #define CMAT(A) int A##r, int A##c, double* A##p
 
+#define KFVEC(A) int A##n, const float*A##p
 #define KDVEC(A) int A##n, const double*A##p
 #define KCVEC(A) int A##n, const double*A##p
+#define KFMAT(A) int A##r, int A##c, const float* A##p
 #define KDMAT(A) int A##r, int A##c, const double* A##p
 #define KCMAT(A) int A##r, int A##c, const double* A##p
 
@@ -55,9 +59,11 @@ typedef short ftnlen;
 int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r));
 int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r));
 
+int transF(KFMAT(x),FMAT(t));
 int transR(KDMAT(x),DMAT(t));
 int transC(KCMAT(x),CMAT(t));
 
+int constantF(float * pval, FVEC(r));
 int constantR(double * pval, DVEC(r));
 int constantC(doublecomplex* pval, CVEC(r));
 
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index c802712..481d72a 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -42,6 +42,17 @@ class (Container c e) => Linear c e where
     equal       :: c e -> c e -> Bool
 
 
+instance Linear Vector Float where
+    scale = vectorMapValF Scale
+    scaleRecip = vectorMapValF Recip
+    addConstant = vectorMapValF AddConstant
+    add = vectorZipF Add
+    sub = vectorZipF Sub
+    mul = vectorZipF Mul
+    divide = vectorZipF Div
+    equal u v = dim u == dim v && vectorFMax (vectorMapF Abs (sub u v)) == 0.0
+    scalar x = fromList [x]
+
 instance Linear Vector Double where
     scale = vectorMapValR Scale
     scaleRecip = vectorMapValR Recip