more refactoring

9 files changed, 89 insertions, 62 deletions

diff --git a/HSSL.cabal b/HSSL.cabal
index 9fd148a..8090d1c 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -14,7 +14,7 @@ tested-with:        GHC ==6.6.1
 Build-Depends:      base, haskell98
 Extensions:         ForeignFunctionInterface
 --ghc-options:        -Wall
-ghc-options:        -O
+ghc-options:        -O0
 hs-source-dirs:     lib
 Exposed-modules:    Data.Packed.Internal,
                     Data.Packed.Internal.Common,
diff --git a/examples/oldtests.hs b/examples/oldtests.hs
index b01675f..f60f4e2 100644
--- a/examples/oldtests.hs
+++ b/examples/oldtests.hs
@@ -5,14 +5,14 @@ import System.Random(randomRs,mkStdGen)
 realMatrix = fromLists :: [[Double]] -> Matrix Double
 realVector = fromList ::  [Double] -> Vector Double
 
-toComplexM = uncurry $ liftMatrix2 (curry toComplex)
+
 
 infixl 2 =~=
 a =~= b = pnorm 1 (flatten (a - b)) < 1E-6
 
 randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed 
 
-randomMatrixC seed (n,m) = toComplexM (randomMatrix seed (n,m), randomMatrix (seed+1) (n,m))
+randomMatrixC seed (n,m) = toComplex (randomMatrix seed (n,m), randomMatrix (seed+1) (n,m))
 
 besselTest = do
     let (r,e) = bessel_J0_e 5.0
@@ -31,7 +31,7 @@ ms = realMatrix [[1,2,3]
 
 ms' = randomMatrix 27 (50,100)
 
-ms'' = toComplexM (randomMatrix 100 (50,100),randomMatrix 101 (50,100))
+ms'' = toComplex (randomMatrix 100 (50,100),randomMatrix 101 (50,100))
 
 fullsvdTest method mat msg = do
     let (u,s,vt) = method mat
@@ -43,7 +43,7 @@ full_svd_Rd = svdRdd
 
 --------------------------------------------------------------------
 
-mcu = toComplexM (randomMatrix 33 (20,20),randomMatrix 34 (20,20))
+mcu = toComplex (randomMatrix 33 (20,20),randomMatrix 34 (20,20))
 
 mcur = randomMatrix 35 (40,40)
 
@@ -53,7 +53,7 @@ eigTest method m msg = do
     assertBool msg $ m <> v =~= v <> diag s
 
 bigmat = m + trans m where m = randomMatrix 18 (1000,1000)
-bigmatc = mc + conjTrans mc where mc = toComplexM(m,m)
+bigmatc = mc + conjTrans mc where mc = toComplex(m,m)
                                   m = randomMatrix 19 (1000,1000)
 
 --------------------------------------------------------------------
@@ -62,22 +62,22 @@ invTest msg m = do
     assertBool msg $ m <> inv m =~= ident (rows m)
 
 invComplexTest msg m = do
-    assertBool msg $ m <> invC m =~= ident (rows m)
+    assertBool msg $ m <> invC m =~= identC (rows m)
 
-invC m = linearSolveC m (ident (rows m))
+invC m = linearSolveC m (identC (rows m))
 
---identC n = toComplexM(ident n, (0::Double) <>ident n)
+identC = comp . ident
 
 --------------------------------------------------------------------
 
 pinvTest f msg m = do
     assertBool msg $ m <> f m <> m =~= m
 
-pinvC m = linearSolveLSC m (ident (rows m))
+pinvC m = linearSolveLSC m (identC (rows m))
 
 pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
 
-pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m))
+pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
 
 --------------------------------------------------------------------
 
diff --git a/examples/tests.hs b/examples/tests.hs
index 0a09ced..2f3b3d8 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -6,10 +6,9 @@
 
 -----------------------------------------------------------------------------
 
-import Data.Packed.Internal
+import Data.Packed.Internal((>|<), fdat, cdat, multiply', multiplyG, MatrixOrder(..))
 import Data.Packed.Vector
 import Data.Packed.Matrix
-import Data.Packed.Internal.Matrix
 import GSL.Vector
 import GSL.Integration
 import GSL.Differentiation
@@ -95,6 +94,8 @@ asC m = (rows m >< cols m) $ toList (cdat m)
 mulC a b = multiply' RowMajor a b
 mulF a b = multiply' ColumnMajor a b
 
+identC = comp . ident
+
 infixl 7 <>
 a <> b = mulF a b
 
@@ -198,15 +199,15 @@ svdTestR fun m = u <> s <> trans v |~| m
 
 
 svdTestC m = u <> s' <> (trans v) |~| m
-                  && u <> conjTrans u |~| ident (rows m)
-                  && v <> conjTrans v |~| ident (cols m)
+                  && u <> conjTrans u |~| identC (rows m)
+                  && v <> conjTrans v |~| identC (cols m)
     where (u,s,v) = svdC m
           s' = liftMatrix comp s
 
 --svdg' m = (u,s',v) where 
 
 eigTestC (SqM m) = (m <> v) |~| (v <> diag s)
-                        && takeDiag (conjTrans v <> v) |~| constant 1 (rows m) --normalized
+                        && takeDiag (conjTrans v <> v) |~| comp (constant 1 (rows m)) --normalized
     where (s,v) = eigC m
 
 eigTestR (SqM m) = (liftMatrix comp m <> v) |~| (v <> diag s)
@@ -218,7 +219,7 @@ eigTestS (Sym m) = (m <> v) |~| (v <> diag s)
     where (s,v) = eigS m
 
 eigTestH (Her m) = (m <> v) |~| (v <> diag (comp s))
-                        && v <> conjTrans v |~| ident (cols m)
+                        && v <> conjTrans v |~| identC (cols m)
     where (s,v) = eigH m
 
 linearSolveSQTest fun singu (PairSM a b) = singu a || (a <> fun a b) |~| b
@@ -248,11 +249,11 @@ identC n = toComplex(ident n, (0::Double) <>ident n)
 pinvTest f m = (m <> f m <> m) |~| m
 
 pinvR m = linearSolveLSR m (ident (rows m))
-pinvC m = linearSolveLSC m (ident (rows m))
+pinvC m = linearSolveLSC m (identC (rows m))
 
 pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
 
-pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m))
+pinvSVDC m = linearSolveSVDC Nothing m (identC (rows m))
 
 --------------------------------------------------------------------
 
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]