Float matrix product

3 files changed, 69 insertions, 21 deletions

diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs
index 8c1c5f6..b81efa4 100644
--- a/lib/Data/Packed/Internal/Signatures.hs
+++ b/lib/Data/Packed/Internal/Signatures.hs
@@ -24,12 +24,15 @@ type PQ = Ptr (Complex Float)                   --
 type PC = Ptr (Complex Double)                  --
 type TF = CInt -> PF -> IO CInt                 --
 type TFF = CInt -> PF -> TF                     --
+type TFV = CInt -> PF -> TV                     --
+type TVF = CInt -> PD -> 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 TFMFMFM =  CInt -> CInt -> PF -> TFMFM     --
 type TM = CInt -> CInt -> PD -> IO CInt         --
 type TMM =  CInt -> CInt -> PD -> TM            --
 type TVMM = CInt -> PD -> TMM                   --
@@ -61,6 +64,7 @@ type TQVQVQV = CInt -> PQ -> TQVQV              --
 type TQVF = CInt -> PQ -> TF                    --
 type TQM = CInt -> CInt -> PQ -> IO CInt        --
 type TQMQM = CInt -> CInt -> PQ -> TQM          --
+type TQMQMQM = CInt -> CInt -> PQ -> TQMQM      --
 type TCMCV = CInt -> CInt -> PC -> TCV          --
 type TVCV = CInt -> PD -> TCV                   --
 type TCVM = CInt -> PC -> TM                    --
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index ac2d0d7..c8cc2c2 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -21,7 +21,7 @@ module Data.Packed.Internal.Vector (
     mapVectorM, mapVectorM_,
     foldVector, foldVectorG, foldLoop,
     createVector, vec,
-    asComplex, asReal,
+    asComplex, asReal, float2DoubleV, double2FloatV,
     fwriteVector, freadVector, fprintfVector, fscanfVector,
     cloneVector,
     unsafeToForeignPtr,
@@ -274,6 +274,24 @@ asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
 asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
     where (fp,i,n) = unsafeToForeignPtr v
 
+---------------------------------------------------------------
+
+float2DoubleV :: Vector Float -> Vector Double
+float2DoubleV v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_float2double vec v vec r "float2double"
+    return r
+
+double2FloatV :: Vector Double -> Vector Float
+double2FloatV v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 c_double2float vec v vec r "double2float2"
+    return r
+
+
+foreign import ccall "float2double" c_float2double:: TFV
+foreign import ccall "double2float" c_double2float:: TVF
+
 ----------------------------------------------------------------
 
 cloneVector :: Storable t => Vector t -> IO (Vector t)
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 8aa1693..8694249 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -1,6 +1,10 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Matrix
@@ -16,8 +20,9 @@
 -----------------------------------------------------------------------------
 
 module Data.Packed.Matrix (
-    Element, Scalar, Container(..), Convert(..),
-    RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf, AutoReal(..),
+    Element, RealElement, Container(..),
+    Convert(..), RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf,
+    AutoReal(..),
     Matrix,rows,cols,
     (><),
     trans,
@@ -51,7 +56,7 @@ import Data.Binary
 import Foreign.Storable
 import Control.Monad(replicateM)
 import Control.Arrow((***))
-import GHC.Float(double2Float,float2Double)
+--import GHC.Float(double2Float,float2Double)
 
 -------------------------------------------------------------------
 
@@ -468,17 +473,32 @@ toBlocksEvery r c m = toBlocks rs cs m where
 
 -- | conversion utilities
 
-class (Element t, Element (Complex t), Fractional t, RealFloat t) => Scalar t
+class (Element t, Element (Complex t), RealFloat t) => RealElement t
+
+instance RealElement Double
+instance RealElement Float
+
+class (Element s, Element d) => Prec s d | s -> d, d -> s where
+    double2FloatG :: Vector d -> Vector s
+    float2DoubleG :: Vector s -> Vector d
+
+instance Prec Float Double where
+    double2FloatG = double2FloatV
+    float2DoubleG = float2DoubleV
+
+instance Prec (Complex Float) (Complex Double) where
+    double2FloatG = asComplex . double2FloatV . asReal
+    float2DoubleG = asComplex . float2DoubleV . asReal
 
-instance Scalar Double
-instance Scalar Float
 
 class Container c where
-    toComplex   :: (Scalar e) => (c e, c e) -> c (Complex e)
-    fromComplex :: (Scalar e) => c (Complex e) -> (c e, c e)
-    comp        :: (Scalar e) => c e -> c (Complex e)
-    conj        :: (Scalar e) => c (Complex e) -> c (Complex e)
+    toComplex   :: (RealElement e) => (c e, c e) -> c (Complex e)
+    fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e)
+    comp        :: (RealElement e) => c e -> c (Complex e)
+    conj        :: (RealElement e) => c (Complex e) -> c (Complex e)
     cmap        :: (Element a, Element b) => (a -> b) -> c a -> c b
+    single      :: Prec a b => c b -> c a
+    double      :: Prec a b => c a -> c b
 
 instance Container Vector where
     toComplex = toComplexV
@@ -486,6 +506,8 @@ instance Container Vector where
     comp v = toComplex (v,constantD 0 (dim v))
     conj = conjV
     cmap = mapVector
+    single = double2FloatG
+    double = float2DoubleG
 
 instance Container Matrix where
     toComplex = uncurry $ liftMatrix2 $ curry toComplex
@@ -494,6 +516,8 @@ instance Container Matrix where
     comp = liftMatrix comp
     conj = liftMatrix conj
     cmap f = liftMatrix (cmap f)
+    single = liftMatrix single
+    double = liftMatrix double
 
 -------------------------------------------------------------------
 
@@ -534,38 +558,40 @@ type instance ElementOf (Matrix a) = a
 
 -------------------------------------------------------------------
 
+-- | generic conversion functions
 class Convert t where
     real'    :: Container c => c (RealOf t) -> c t
     complex' :: Container c => c t -> c (ComplexOf t)
-    single   :: Container c => c t -> c (SingleOf t)
-    double   :: Container c => c t -> c (DoubleOf t)
+    single'  :: Container c => c t -> c (SingleOf t)
+    double'  :: Container c => c t -> c (DoubleOf t)
 
 instance Convert Double where
     real' = id
     complex' = comp
-    single = cmap double2Float
-    double = id
+    single' = single
+    double' = id
 
 instance Convert Float where
     real' = id
     complex' = comp
-    single = id
-    double = cmap float2Double
+    single' = id
+    double' = double
 
 instance Convert (Complex Double) where
     real' = comp
     complex' = id
-    single = toComplex . (single *** single) . fromComplex
-    double = id
+    single' = single
+    double' = id
 
 instance Convert (Complex Float) where
     real' = comp
     complex' = id
-    single = id
-    double = toComplex . (double *** double) . fromComplex
+    single' = id
+    double' = double
 
 -------------------------------------------------------------------
 
+-- | to be replaced by Convert
 class AutoReal t where
     real :: Container c => c Double -> c t
     complex :: Container c => c t -> c (Complex Double)