code refactoring

6 files changed, 177 insertions, 158 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index 8090d1c..02f5129 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -48,8 +48,10 @@ Exposed-modules:    Data.Packed.Internal,
                     LAPACK,
                     LinearAlgebra,
                     LinearAlgebra.Linear,
+                    LinearAlgebra.Instances,
                     LinearAlgebra.Algorithms
-                    , GSLHaskell, Graphics.Plot
+                    , Graphics.Plot
+--                    , GSLHaskell
 Other-modules:
 C-sources:          lib/Data/Packed/Internal/aux.c,
                     lib/LAPACK/lapack-aux.c,
diff --git a/examples/tests.hs b/examples/tests.hs
index 2f3b3d8..f352b71 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -20,7 +20,7 @@ import Test.QuickCheck
 import Test.HUnit hiding ((~:))
 import Complex
 import LinearAlgebra.Algorithms
-import LinearAlgebra.Linear
+import LinearAlgebra.Linear hiding ((<>))
 import GSL.Matrix
 import GSLHaskell hiding ((<>),constant)
 
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 8bca510..951cec6 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -172,7 +172,7 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
 ----------------------------------------------------------------
 
 -- | Optimized matrix computations are provided for elements in the Field class.
-class (Storable a, Num a) => Field a where
+class (Storable a, Floating a) => Field a where
     constantD :: a -> Int -> Vector a
     transdata :: Int -> Vector a -> Int -> Vector a
     multiplyD :: Matrix a -> Matrix a -> Matrix a
diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs
index 3158458..254a957 100644
--- a/lib/GSLHaskell.hs
+++ b/lib/GSLHaskell.hs
@@ -9,7 +9,7 @@ Maintainer  :  Alberto Ruiz (aruiz at um dot es)
 Stability   :  provisional
 Portability :  uses -fffi and -fglasgow-exts
 
-GSLHaskell interface, with reasonable numeric instances for Vectors and Matrices. In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
+Old GSLHaskell interface.
 
 -}
 -----------------------------------------------------------------------------
@@ -46,7 +46,6 @@ import GSL.Special(setErrorHandlerOff,
     bessel_J0_e,
     exp_e10_e,
     gamma)
---import Data.Packed.Internal hiding (dsp,comp)
 import Data.Packed.Vector
 import Data.Packed.Matrix
 import Data.Packed.Matrix hiding ((><))
@@ -55,163 +54,14 @@ import qualified LinearAlgebra.Algorithms
 import LAPACK
 import GSL.Matrix
 import LinearAlgebra.Algorithms hiding (pnorm)
-import LinearAlgebra.Linear
+import LinearAlgebra.Linear hiding (Mul,(<>))
+import Data.Packed.Internal.Matrix(multiply)
 import Complex
 import Numeric(showGFloat)
 import Data.List(transpose,intersperse)
 import Foreign(Storable)
 import Data.Array
-
-
-adaptScalar f1 f2 f3 x y
-    | dim x == 1 = f1   (x@>0) y
-    | dim y == 1 = f3 x (y@>0)
-    | otherwise = f2 x y
-
-liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
-liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
-                     | otherwise    = error "nonconformant matrices in liftMatrix2'"
-
-compat' :: Matrix a -> Matrix b -> Bool
-compat' m1 m2 = rows m1 == 1 && cols m1 == 1
-             || rows m2 == 1 && cols m2 == 1
-             || rows m1 == rows m2 && cols m1 == cols m2
-
-instance (Eq a, Field a) => Eq (Vector a) where
-    a == b = dim a == dim b && toList a == toList b
-
-instance (Linear Vector a) => Num (Vector a) where
-    (+) = adaptScalar addConstant add (flip addConstant)
-    negate = scale (-1)
-    (*) = adaptScalar scale mul (flip scale)
-    signum = liftVector signum
-    abs = liftVector abs
-    fromInteger = fromList . return . fromInteger
-
-instance (Eq a, Field a) => Eq (Matrix a) where
-    a == b = cols a == cols b && flatten a == flatten b
-
-instance (Field a, Linear Vector a) => Num (Matrix a) where
-    (+) = liftMatrix2' (+)
-    (-) = liftMatrix2' (-)
-    negate = liftMatrix negate
-    (*) = liftMatrix2' (*)
-    signum = liftMatrix signum
-    abs = liftMatrix abs
-    fromInteger = (1><1) . return . fromInteger
-
----------------------------------------------------
-
-instance Fractional (Vector Double) where
-    fromRational n = fromList [fromRational n]
-    (/) = adaptScalar f (vectorZipR Div) g where
-        r `f` v = vectorMapValR Recip r v
-        v `g` r = scale (recip r) v
-
--------------------------------------------------------
-
-instance Fractional (Vector (Complex Double)) where
-    fromRational n = fromList [fromRational n]
-    (/) = adaptScalar f (vectorZipC Div) g where
-        r `f` v = vectorMapValC Recip r v
-        v `g` r = scale (recip r) v
-
-------------------------------------------------------
-
-instance Fractional (Matrix Double) where
-    fromRational n = (1><1) [fromRational n]
-    (/) = liftMatrix2' (/)
-
--------------------------------------------------------
-
-instance Fractional (Matrix (Complex Double)) where
-    fromRational n = (1><1) [fromRational n]
-    (/) = liftMatrix2' (/)
-
----------------------------------------------------------
-
-instance Floating (Vector Double) where
-    sin   = vectorMapR Sin
-    cos   = vectorMapR Cos
-    tan   = vectorMapR Tan
-    asin  = vectorMapR ASin
-    acos  = vectorMapR ACos
-    atan  = vectorMapR ATan
-    sinh  = vectorMapR Sinh
-    cosh  = vectorMapR Cosh
-    tanh  = vectorMapR Tanh
-    asinh = vectorMapR ASinh
-    acosh = vectorMapR ACosh
-    atanh = vectorMapR ATanh
-    exp   = vectorMapR Exp
-    log   = vectorMapR Log
-    sqrt  = vectorMapR Sqrt
-    (**)  = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
-    pi    = fromList [pi]
-
------------------------------------------------------------
-
-instance Floating (Matrix Double) where
-    sin   = liftMatrix sin
-    cos   = liftMatrix cos
-    tan   = liftMatrix tan
-    asin  = liftMatrix asin
-    acos  = liftMatrix acos
-    atan  = liftMatrix atan
-    sinh  = liftMatrix sinh
-    cosh  = liftMatrix cosh
-    tanh  = liftMatrix tanh
-    asinh = liftMatrix asinh
-    acosh = liftMatrix acosh
-    atanh = liftMatrix atanh
-    exp   = liftMatrix exp
-    log   = liftMatrix log
-    (**)  = liftMatrix2' (**)
-    sqrt  = liftMatrix sqrt
-    pi    = (1><1) [pi]
--------------------------------------------------------------
-
-instance Floating (Vector (Complex Double)) where
-    sin   = vectorMapC Sin
-    cos   = vectorMapC Cos
-    tan   = vectorMapC Tan
-    asin  = vectorMapC ASin
-    acos  = vectorMapC ACos
-    atan  = vectorMapC ATan
-    sinh  = vectorMapC Sinh
-    cosh  = vectorMapC Cosh
-    tanh  = vectorMapC Tanh
-    asinh = vectorMapC ASinh
-    acosh = vectorMapC ACosh
-    atanh = vectorMapC ATanh
-    exp   = vectorMapC Exp
-    log   = vectorMapC Log
-    sqrt  = vectorMapC Sqrt
-    (**)  = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
-    pi    = fromList [pi]
-
----------------------------------------------------------------
-
-instance Floating (Matrix (Complex Double)) where
-    sin   = liftMatrix sin
-    cos   = liftMatrix cos
-    tan   = liftMatrix tan
-    asin  = liftMatrix asin
-    acos  = liftMatrix acos
-    atan  = liftMatrix atan
-    sinh  = liftMatrix sinh
-    cosh  = liftMatrix cosh
-    tanh  = liftMatrix tanh
-    asinh = liftMatrix asinh
-    acosh = liftMatrix acosh
-    atanh = liftMatrix atanh
-    exp   = liftMatrix exp
-    log   = liftMatrix log
-    (**)  = liftMatrix2' (**)
-    sqrt  = liftMatrix sqrt
-    pi    = (1><1) [pi]
-
----------------------------------------------------------------
+import LinearAlgebra.Instances
 
 
 class Mul a b c | a b -> c where
diff --git a/lib/LinearAlgebra/Instances.hs b/lib/LinearAlgebra/Instances.hs
new file mode 100644
index 0000000..3dbe5a7
--- /dev/null
+++ b/lib/LinearAlgebra/Instances.hs
@@ -0,0 +1,140 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  LinearAlgebra.Instances
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Numeric instances for Vector and Matrix.
+
+In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
+
+-}
+-----------------------------------------------------------------------------
+
+module LinearAlgebra.Instances(
+) where
+
+import LinearAlgebra.Linear
+import GSL.Vector
+import Data.Packed.Matrix
+import Data.Packed.Vector
+import Complex
+
+adaptScalar f1 f2 f3 x y
+    | dim x == 1 = f1   (x@>0) y
+    | dim y == 1 = f3 x (y@>0)
+    | otherwise = f2 x y
+
+liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
+                     | otherwise    = error "nonconformant matrices in liftMatrix2'"
+
+compat' :: Matrix a -> Matrix b -> Bool
+compat' m1 m2 = rows m1 == 1 && cols m1 == 1
+             || rows m2 == 1 && cols m2 == 1
+             || rows m1 == rows m2 && cols m1 == cols m2
+
+instance (Eq a, Field a) => Eq (Vector a) where
+    a == b = dim a == dim b && toList a == toList b
+
+instance (Linear Vector a) => Num (Vector a) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = liftVector signum
+    abs = liftVector abs
+    fromInteger = fromList . return . fromInteger
+
+instance (Eq a, Field a) => Eq (Matrix a) where
+    a == b = cols a == cols b && flatten a == flatten b
+
+instance (Linear Vector a) => Num (Matrix a) where
+    (+) = liftMatrix2' (+)
+    (-) = liftMatrix2' (-)
+    negate = liftMatrix negate
+    (*) = liftMatrix2' (*)
+    signum = liftMatrix signum
+    abs = liftMatrix abs
+    fromInteger = (1><1) . return . fromInteger
+
+---------------------------------------------------
+
+instance (Linear Vector a) => Fractional (Vector a) where
+    fromRational n = fromList [fromRational n]
+    (/) = adaptScalar f divide g where
+        r `f` v = scaleRecip r v
+        v `g` r = scale (recip r) v
+
+-------------------------------------------------------
+
+instance (Linear Vector a, Fractional (Vector a)) => Fractional (Matrix a) where
+    fromRational n = (1><1) [fromRational n]
+    (/) = liftMatrix2' (/)
+
+---------------------------------------------------------
+
+instance Floating (Vector Double) where
+    sin   = vectorMapR Sin
+    cos   = vectorMapR Cos
+    tan   = vectorMapR Tan
+    asin  = vectorMapR ASin
+    acos  = vectorMapR ACos
+    atan  = vectorMapR ATan
+    sinh  = vectorMapR Sinh
+    cosh  = vectorMapR Cosh
+    tanh  = vectorMapR Tanh
+    asinh = vectorMapR ASinh
+    acosh = vectorMapR ACosh
+    atanh = vectorMapR ATanh
+    exp   = vectorMapR Exp
+    log   = vectorMapR Log
+    sqrt  = vectorMapR Sqrt
+    (**)  = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
+    pi    = fromList [pi]
+
+-------------------------------------------------------------
+
+instance Floating (Vector (Complex Double)) where
+    sin   = vectorMapC Sin
+    cos   = vectorMapC Cos
+    tan   = vectorMapC Tan
+    asin  = vectorMapC ASin
+    acos  = vectorMapC ACos
+    atan  = vectorMapC ATan
+    sinh  = vectorMapC Sinh
+    cosh  = vectorMapC Cosh
+    tanh  = vectorMapC Tanh
+    asinh = vectorMapC ASinh
+    acosh = vectorMapC ACosh
+    atanh = vectorMapC ATanh
+    exp   = vectorMapC Exp
+    log   = vectorMapC Log
+    sqrt  = vectorMapC Sqrt
+    (**)  = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
+    pi    = fromList [pi]
+
+-----------------------------------------------------------
+
+instance (Linear Vector a, Floating (Vector a)) => Floating (Matrix a) where
+    sin   = liftMatrix sin
+    cos   = liftMatrix cos
+    tan   = liftMatrix tan
+    asin  = liftMatrix asin
+    acos  = liftMatrix acos
+    atan  = liftMatrix atan
+    sinh  = liftMatrix sinh
+    cosh  = liftMatrix cosh
+    tanh  = liftMatrix tanh
+    asinh = liftMatrix asinh
+    acosh = liftMatrix acosh
+    atanh = liftMatrix atanh
+    exp   = liftMatrix exp
+    log   = liftMatrix log
+    (**)  = liftMatrix2' (**)
+    sqrt  = liftMatrix sqrt
+    pi    = (1><1) [pi]
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index a2071ed..148bbd3 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -15,7 +15,8 @@ Portability :  uses ffi
 
 module LinearAlgebra.Linear (
     Linear(..),
-    multiply, dot, outer
+    dot, outer,
+    Mul(..)
 ) where
 
 
@@ -27,10 +28,12 @@ import Complex
 
 class (Field e) => Linear c e where
     scale       :: e -> c e -> c e
+    scaleRecip  :: e -> c e -> c e
     addConstant :: e -> c e -> c e
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
     mul         :: c e -> c e -> c e
+    divide      :: c e -> c e -> c e
     toComplex   :: RealFloat e => (c e, c e) -> c (Complex e)
     fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
     comp        :: RealFloat e => c e -> c (Complex e)
@@ -38,10 +41,12 @@ class (Field e) => Linear c e where
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
+    scaleRecip = vectorMapValR Recip
     addConstant = vectorMapValR AddConstant
     add = vectorZipR Add
     sub = vectorZipR Sub
     mul = vectorZipR Mul
+    divide = vectorZipR Div
     toComplex = Data.Packed.Internal.toComplex
     fromComplex = Data.Packed.Internal.fromComplex
     comp = Data.Packed.Internal.comp
@@ -49,10 +54,12 @@ instance Linear Vector Double where
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
+    scaleRecip = vectorMapValC Recip
     addConstant = vectorMapValC AddConstant
     add = vectorZipC Add
     sub = vectorZipC Sub
     mul = vectorZipC Mul
+    divide = vectorZipC Div
     toComplex = undefined -- can't match
     fromComplex = undefined
     comp = undefined
@@ -60,10 +67,12 @@ instance Linear Vector (Complex Double) where
 
 instance Linear Matrix Double where
     scale x = liftMatrix (scale x)
+    scaleRecip x = liftMatrix (scaleRecip x)
     addConstant x = liftMatrix (addConstant x)
     add = liftMatrix2 add
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
+    divide = liftMatrix2 divide
     toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex
     fromComplex z = (reshape c r, reshape c i)
         where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z)
@@ -73,10 +82,12 @@ instance Linear Matrix Double where
 
 instance Linear Matrix (Complex Double) where
     scale x = liftMatrix (scale x)
+    scaleRecip x = liftMatrix (scaleRecip x)
     addConstant x = liftMatrix (addConstant x)
     add = liftMatrix2 add
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
+    divide = liftMatrix2 divide
     toComplex = undefined
     fromComplex = undefined
     comp = undefined
@@ -102,3 +113,19 @@ dot u v = dat (multiply r c) `at` 0
 -}
 outer :: (Field t) => Vector t -> Vector t -> Matrix t
 outer u v = asColumn u `multiply` asRow v
+
+
+class Mul a b c | a b -> c where
+ infixl 7 <>
+ -- | matrix product
+ (<>) :: Field t => a t -> b t -> c t
+
+instance Mul Matrix Matrix Matrix where
+    (<>) = multiply
+
+instance Mul Matrix Vector Vector where
+    (<>) m v = flatten $ m <> (asColumn v)
+
+instance Mul Vector Matrix Vector where
+    (<>) v m = flatten $ (asRow v) <> m
+