From 97e8a48be58fd53afccc7ae01ee6ec5805d5c1cd Mon Sep 17 00:00:00 2001
From: Vivian McPhail <haskell.vivian.mcphail@gmail.com>
Date: Thu, 8 Jul 2010 23:03:48 +0000
Subject: Linear and Floating (Complex Float)

---
 lib/Data/Packed/Internal/Matrix.hs            |  10 ++
 lib/Data/Packed/Internal/Signatures.hs        |   2 +
 lib/Data/Packed/Matrix.hs                     |   9 ++
 lib/Numeric/GSL/Vector.hs                     |  24 ++++-
 lib/Numeric/GSL/gsl-aux.c                     | 126 ++++++++++++++++++++++++++
 lib/Numeric/LinearAlgebra/Instances.hs        |  29 ++++++
 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c |  22 +++++
 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h |   6 ++
 lib/Numeric/LinearAlgebra/Linear.hs           |  11 +++
 9 files changed, 236 insertions(+), 3 deletions(-)

(limited to 'lib')

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 7b3b305..861c72a 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -265,6 +265,10 @@ instance Element Double where
     transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
 
+instance Element (Complex Float) where
+    transdata  = transdataAux ctransQ
+    constantD  = constantAux cconstantQ
+
 instance Element (Complex Double) where
     transdata  = transdataAux ctransC
     constantD  = constantAux cconstantC
@@ -314,6 +318,7 @@ transdataAux fun c1 d c2 =
 
 foreign import ccall "transF" ctransF :: TFMFM
 foreign import ccall "transR" ctransR :: TMM
+foreign import ccall "transQ" ctransQ :: TQMQM
 foreign import ccall "transC" ctransC :: TCMCM
 ----------------------------------------------------------------------
 
@@ -342,9 +347,14 @@ constantR :: Double -> Int -> Vector Double
 constantR = constantAux cconstantR
 foreign import ccall "constantR" cconstantR :: Ptr Double -> TV
 
+constantQ :: Complex Float -> Int -> Vector (Complex Float)
+constantQ = constantAux cconstantQ
+foreign import ccall "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
+
 constantC :: Complex Double -> Int -> Vector (Complex Double)
 constantC = constantAux cconstantC
 foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV
+
 ----------------------------------------------------------------------
 
 -- | Extracts a submatrix from a matrix.
diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs
index 78d00fa..8c1c5f6 100644
--- a/lib/Data/Packed/Internal/Signatures.hs
+++ b/lib/Data/Packed/Internal/Signatures.hs
@@ -59,6 +59,8 @@ type TQV = CInt -> PQ -> IO CInt                --
 type TQVQV = CInt -> PQ -> TQV                  --
 type TQVQVQV = CInt -> PQ -> TQVQV              --
 type TQVF = CInt -> PQ -> TF                    --
+type TQM = CInt -> CInt -> PQ -> IO CInt        --
+type TQMQM = CInt -> CInt -> PQ -> TQM          --
 type TCMCV = CInt -> CInt -> PC -> TCV          --
 type TVCV = CInt -> PD -> TCV                   --
 type TCVM = CInt -> PC -> TM                    --
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index c6d8a90..e7ee781 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -452,6 +452,7 @@ class (Element e) => Container c e where
     fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
     comp        :: RealFloat e => c e -> c (Complex e)
     conj        :: RealFloat e => c (Complex e) -> c (Complex e)
+    -- these next two are now weird given we have Floats as well
     real        :: c Double -> c e
     complex     :: c e -> c (Complex Double)
 
@@ -471,6 +472,14 @@ instance Container Vector Double where
     real = id
     complex = comp
 
+instance Container Vector (Complex Float) where
+    toComplex = undefined -- can't match
+    fromComplex = undefined
+    comp = undefined
+    conj = undefined
+    real = comp . mapVector realToFrac
+    complex = mapVector (\(r :+ i) -> realToFrac r :+ realToFrac i)
+
 instance Container Vector (Complex Double) where
     toComplex = undefined -- can't match
     fromComplex = undefined
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 0148c4f..14ba0ff 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -17,9 +17,9 @@ module Numeric.GSL.Vector (
     sumF, sumR, sumQ, sumC,
     dotF, dotR, dotQ, dotC,
     FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
-    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF,
-    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF,
-    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF,
+    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
+    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
+    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
     RandDist(..), randomVector
 ) where
 
@@ -214,6 +214,12 @@ vectorMapF = vectorMapAux c_vectorMapF
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV
+
 -------------------------------------------------------------------
 
 -- | map of real vectors with given function
@@ -234,6 +240,12 @@ vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV
+
 -------------------------------------------------------------------
 
 -- | elementwise operation on real vectors
@@ -254,6 +266,12 @@ vectorZipF = vectorZipAux c_vectorZipF
 
 foreign import ccall safe "gsl-aux.h 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 safe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV
+
 -----------------------------------------------------------------------
 
 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 3f9eeba..689989d 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -393,6 +393,83 @@ int mapC(int code, KCVEC(x), CVEC(r)) {
 }
 
 
+gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
+{
+  gsl_complex c;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c.dat[0] = a.dat[0];
+  c.dat[1] = a.dat[1];
+
+  r = (*cf)(c);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), 
+					gsl_complex_float a,gsl_complex_float b)
+{
+  gsl_complex c1;
+  gsl_complex c2;
+  gsl_complex r;
+
+  gsl_complex_float float_r;
+
+  c1.dat[0] = a.dat[0];
+  c1.dat[1] = a.dat[1];
+
+  c2.dat[0] = b.dat[0];
+  c2.dat[1] = b.dat[1];
+
+  r = (*cf)(c1,c2);
+
+  float_r.dat[0] = r.dat[0];
+  float_r.dat[1] = r.dat[1];
+
+  return float_r;
+}
+
+#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
+#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
+int mapQAux(int code, KGQVEC(x), GQVEC(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)
+        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)
+        OPC(16,gsl_complex_sqrt)
+
+        // gsl_complex_arg
+        // gsl_complex_abs
+        default: ERROR(BAD_CODE);
+    }
+}
+
+int mapQ(int code, KQVEC(x), QVEC(r)) {
+    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
 int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
     int k;
     double val = *pval;
@@ -446,6 +523,27 @@ int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
 }
 
 
+int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
+    int k;
+    gsl_complex_float 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,gsl_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);
+    }
+}
+
+int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
+    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
+}
+
+
 #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 }
 int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
@@ -519,6 +617,34 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
 }
 
 
+#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
+int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
+    REQUIRES(an == bn && an == rn, BAD_SIZE);
+    int k;
+    switch(code) {
+        OPCZE(0,"zipQ Add",gsl_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)
+    }
+    //KCVVIEW(a);
+    //KCVVIEW(b);
+    //CVVIEW(r);
+    //gsl_vector_memcpy(V(r),V(a));
+    //int res;
+    switch(code) {
+        default: ERROR(BAD_CODE);
+    }
+}
+
+
+int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
+    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
+}
+
+
 
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index bba89c8..04a9d88 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -118,6 +118,14 @@ instance Num (Vector (Complex Double)) where
     abs = vectorMapC Abs
     fromInteger = fromList . return . fromInteger
 
+instance Num (Vector (Complex Float)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapQ Sign
+    abs = vectorMapQ Abs
+    fromInteger = fromList . return . fromInteger
+
 instance Linear Matrix a => Eq (Matrix a) where
     (==) = equal
 
@@ -209,6 +217,27 @@ instance Floating (Vector (Complex Double)) where
 
 -----------------------------------------------------------
 
+instance Floating (Vector (Complex Float)) where
+    sin   = vectorMapQ Sin
+    cos   = vectorMapQ Cos
+    tan   = vectorMapQ Tan
+    asin  = vectorMapQ ASin
+    acos  = vectorMapQ ACos
+    atan  = vectorMapQ ATan
+    sinh  = vectorMapQ Sinh
+    cosh  = vectorMapQ Cosh
+    tanh  = vectorMapQ Tanh
+    asinh = vectorMapQ ASinh
+    acosh = vectorMapQ ACosh
+    atanh = vectorMapQ ATanh
+    exp   = vectorMapQ Exp
+    log   = vectorMapQ Log
+    sqrt  = vectorMapQ Sqrt
+    (**)  = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
+    pi    = fromList [pi]
+
+-----------------------------------------------------------
+
 instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
     sin   = liftMatrix sin
     cos   = liftMatrix cos
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index b9c2572..7a40991 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -1063,6 +1063,18 @@ int transR(KDMAT(x),DMAT(t)) {
     OK
 }
 
+int transQ(KQMAT(x),QMAT(t)) {
+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);
+    DEBUGMSG("transQ");
+    int i,j;
+    for (i=0; i<tr; i++) {
+        for (j=0; j<tc; j++) {
+        ((complex*)tp)[i*tc+j] = ((complex*)xp)[j*xc+i];
+        }
+    }
+    OK
+}
+
 int transC(KCMAT(x),CMAT(t)) {
     REQUIRES(xr==tc && xc==tr,BAD_SIZE);
     DEBUGMSG("transC");
@@ -1097,6 +1109,16 @@ int constantR(double * pval, DVEC(r)) {
     OK
 }
 
+int constantQ(complex* pval, QVEC(r)) {
+    DEBUGMSG("constantQ")
+    int k;
+    complex val = *pval;
+    for(k=0;k<rn;k++) {
+        ((complex*)rp)[k]=val;
+    }
+    OK
+}
+
 int constantC(doublecomplex* pval, CVEC(r)) {
     DEBUGMSG("constantC")
     int k;
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 415a6ab..d01d9e5 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -42,16 +42,20 @@ typedef short ftnlen;
 
 #define FVEC(A) int A##n, float*A##p
 #define DVEC(A) int A##n, double*A##p
+#define QVEC(A) int A##n, float*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 QMAT(A) int A##r, int A##c, float* 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 KQVEC(A) int A##n, const float*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 KQMAT(A) int A##r, int A##c, const float* A##p
 #define KCMAT(A) int A##r, int A##c, const double* A##p
 
 /********************************************************/
@@ -61,10 +65,12 @@ 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 transQ(KQMAT(x),QMAT(t));
 int transC(KCMAT(x),CMAT(t));
 
 int constantF(float * pval, FVEC(r));
 int constantR(double * pval, DVEC(r));
+int constantQ(complex* pval, QVEC(r));
 int constantC(doublecomplex* pval, CVEC(r));
 
 int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index aed6a2b..2351ff1 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -132,6 +132,17 @@ instance Linear Vector (Complex Double) where
     equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0
     scalar x = fromList [x]
 
+instance Linear Vector (Complex Float) where
+    scale = vectorMapValQ Scale
+    scaleRecip = vectorMapValQ Recip
+    addConstant = vectorMapValQ AddConstant
+    add = vectorZipQ Add
+    sub = vectorZipQ Sub
+    mul = vectorZipQ Mul
+    divide = vectorZipQ Div
+    equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0
+    scalar x = fromList [x]
+
 instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where
     scale x = liftMatrix (scale x)
     scaleRecip x = liftMatrix (scaleRecip x)
-- 
cgit v1.2.3