more refactoring

6 files changed, 69 insertions, 43 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 89a162a..58b325c 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -348,7 +348,7 @@ foreign import ccall safe "aux.h constantC"
 @> constant 2 7
 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
 -}
-constant :: Field a => a -> Int -> Vector a
+constant :: Double -> Int -> Vector Double
 constant = constantD
 
 --------------------------------------------------------------------------
@@ -361,6 +361,11 @@ conj v = asComplex $ cdat $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1
 toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
 toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i]
 
+-- | the inverse of 'toComplex'
+fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
+fromComplex z = (r,i) where
+    [r,i] = toColumns $ reshape 2 $ asReal z
+
 -- | converts a real vector into a complex representation (with zero imaginary parts)
 comp :: Vector Double -> Vector (Complex Double)
 comp v = toComplex (v,constant 0 (dim v))
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 01e8133..fc08ce4 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -75,12 +75,12 @@ diagRect s r c
     | r == c    = diag s
     | r < c     = trans $ diagRect s c r
     | r > c     = joinVert  [diag s , zeros (r-c,c)]
-    where zeros (r,c) = reshape c $ constant 0 (r*c)
+    where zeros (r,c) = reshape c $ constantD 0 (r*c)
 
 takeDiag :: (Field t) => Matrix t -> Vector t
 takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
 
-ident :: (Num t, Field t) => Int -> Matrix t
+ident :: Int -> Matrix Double
 ident n = diag (constant 1 n)
 
 ------------------------------------------------------------
diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
index 15710df..ec8ceea 100644
--- a/lib/GSL/Matrix.hs
+++ b/lib/GSL/Matrix.hs
@@ -289,7 +289,7 @@ luC m = (l,u,p, fromIntegral s') where
     lu = reshape r $ subVector 0 (r*r) v
     s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v
     u = triang r r 0 1 `mul` lu
-    l = (triang r r 0 0 `mul` lu) `add` ident r
+    l = (triang r r 0 0 `mul` lu) `add` liftMatrix comp (ident r)
     add = liftMatrix2 $ vectorZipC Add
     mul = liftMatrix2 $ vectorZipC Mul
 
diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs
index e65ff28..3158458 100644
--- a/lib/GSLHaskell.hs
+++ b/lib/GSLHaskell.hs
@@ -18,6 +18,7 @@ module GSLHaskell(
     module Data.Packed.Vector,
     module Data.Packed.Matrix,
     module LinearAlgebra.Algorithms,
+    module LinearAlgebra.Linear,
     module LAPACK,
     module GSL.Integration,
     module GSL.Differentiation,
@@ -28,18 +29,15 @@ module GSLHaskell(
     module GSL.Special,
     module Graphics.Plot,
     module Complex,
-    Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->), GSLHaskell.constant,
-    fromArray2D, fromComplex, toComplex, GSLHaskell.pnorm, scale, outer
+    Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->),
+    fromArray2D, GSLHaskell.pnorm,
 ) where
 
-
-import LAPACK
 import GSL.Integration
 import GSL.Differentiation
 import GSL.Fourier
 import GSL.Polynomials
 import GSL.Minimization
-import GSL.Matrix
 import Graphics.Plot
 import Complex
 import GSL.Special(setErrorHandlerOff,
@@ -48,14 +46,16 @@ import GSL.Special(setErrorHandlerOff,
     bessel_J0_e,
     exp_e10_e,
     gamma)
-import Data.Packed.Internal hiding (dsp)
-import Data.Packed.Vector hiding (constant)
+--import Data.Packed.Internal hiding (dsp,comp)
+import Data.Packed.Vector
 import Data.Packed.Matrix
 import Data.Packed.Matrix hiding ((><))
 import GSL.Vector
-import LinearAlgebra.Linear
 import qualified LinearAlgebra.Algorithms
+import LAPACK
+import GSL.Matrix
 import LinearAlgebra.Algorithms hiding (pnorm)
+import LinearAlgebra.Linear
 import Complex
 import Numeric(showGFloat)
 import Data.List(transpose,intersperse)
@@ -69,7 +69,7 @@ adaptScalar f1 f2 f3 x y
     | 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 (cdat m1) (cdat m2))
+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
@@ -89,7 +89,7 @@ instance (Linear Vector a) => Num (Vector a) where
     fromInteger = fromList . return . fromInteger
 
 instance (Eq a, Field a) => Eq (Matrix a) where
-    a == b = rows a == rows b && cols a == cols b && cdat a == cdat b && fdat a == fdat b
+    a == b = cols a == cols b && flatten a == flatten b
 
 instance (Field a, Linear Vector a) => Num (Matrix a) where
     (+) = liftMatrix2' (+)
@@ -377,9 +377,6 @@ size = dim
 gmap :: (Storable a, Storable b) => (a->b) -> Vector a -> Vector b
 gmap f v = liftVector f v
 
-constant :: Double -> Int -> Vector Double
-constant = constantR
-
 -- shows a Double with n digits after the decimal point    
 shf :: (RealFloat a) => Int -> a -> String     
 shf dec n | abs n < 1e-10 = "0."
@@ -473,27 +470,6 @@ fromArray2D m = (r><c) (elems m)
           r = r1-r0+1
           c = c1-c0+1
 
--- | creates a complex vector from vectors with real and imaginary parts
-toComplexV :: (Vector Double, Vector Double) ->  Vector (Complex Double)
-toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
-
--- | extracts the real and imaginary parts of a complex vector
-fromComplexV :: Vector (Complex Double) -> (Vector Double, Vector Double)
-fromComplexV m = (a,b) where [a,b] = toColumns $ reshape 2 $ asReal m
-
--- | creates a complex matrix from matrices with real and imaginary parts
-toComplexM :: (Matrix Double, Matrix Double) ->  Matrix (Complex Double)
-toComplexM (r,i) = reshape (cols r) $ asComplex $ flatten $ fromColumns [flatten r, flatten i]
-
--- | extracts the real and imaginary parts of a complex matrix
-fromComplexM :: Matrix (Complex Double) -> (Matrix Double, Matrix Double)
-fromComplexM m = (reshape c a, reshape c b)
-    where c = cols m
-          [a,b] = toColumns $ reshape 2 $ asReal $ flatten m
-
-fromComplex :: Matrix (Complex Double) -> (Matrix Double, Matrix Double)
-fromComplex = fromComplexM
-
 
 pnorm :: (Normed t1, Num t) => t -> t1 -> Double
 pnorm 0 = LinearAlgebra.Algorithms.pnorm Infinity
diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs
index b0c8b9d..3b56fc4 100644
--- a/lib/LinearAlgebra.hs
+++ b/lib/LinearAlgebra.hs
@@ -13,6 +13,19 @@ Some linear algebra algorithms, implemented by means of BLAS, LAPACK or GSL.
 -}
 -----------------------------------------------------------------------------
 module LinearAlgebra (
-
+    module Data.Packed.Vector,
+    module Data.Packed.Matrix,
+    module LinearAlgebra.Linear,
+    module LAPACK,
+    module GSL.Matrix,
+    module LinearAlgebra.Algorithms,
+    module Complex
 ) where
 
+import LinearAlgebra.Linear
+import LinearAlgebra.Algorithms
+import LAPACK
+import GSL.Matrix
+import Data.Packed.Matrix
+import Data.Packed.Vector
+import Complex
\ No newline at end of file
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index 53e011f..a2071ed 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -15,8 +15,6 @@ Portability :  uses ffi
 
 module LinearAlgebra.Linear (
     Linear(..),
-    toComplex, comp,
-    conj,
     multiply, dot, outer
 ) where
 
@@ -33,6 +31,10 @@ class (Field e) => Linear c e where
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
     mul         :: 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)
+    conj        :: RealFloat e => c (Complex e) -> c (Complex e)
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
@@ -40,6 +42,10 @@ instance Linear Vector Double where
     add = vectorZipR Add
     sub = vectorZipR Sub
     mul = vectorZipR Mul
+    toComplex = Data.Packed.Internal.toComplex
+    fromComplex = Data.Packed.Internal.fromComplex
+    comp = Data.Packed.Internal.comp
+    conj = Data.Packed.Internal.conj
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
@@ -47,6 +53,34 @@ instance Linear Vector (Complex Double) where
     add = vectorZipC Add
     sub = vectorZipC Sub
     mul = vectorZipC Mul
+    toComplex = undefined -- can't match
+    fromComplex = undefined
+    comp = undefined
+    conj = undefined
+
+instance Linear Matrix Double where
+    scale x = liftMatrix (scale x)
+    addConstant x = liftMatrix (addConstant x)
+    add = liftMatrix2 add
+    sub = liftMatrix2 sub
+    mul = liftMatrix2 mul
+    toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex
+    fromComplex z = (reshape c r, reshape c i)
+        where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z)
+              c = cols z
+    comp = liftMatrix Data.Packed.Internal.comp
+    conj = liftMatrix Data.Packed.Internal.conj
+
+instance Linear Matrix (Complex Double) where
+    scale x = liftMatrix (scale x)
+    addConstant x = liftMatrix (addConstant x)
+    add = liftMatrix2 add
+    sub = liftMatrix2 sub
+    mul = liftMatrix2 mul
+    toComplex = undefined
+    fromComplex = undefined
+    comp = undefined
+    conj = undefined
 
 --------------------------------------------------
 
@@ -58,8 +92,6 @@ dot u v = dat (multiply r c) `at` 0
           c = asColumn v
 
 
-
-
 {- | Outer product of two vectors.
 
 @\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]