6 files changed, 235 insertions, 145 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
new file mode 100644
index 0000000..dddd269
--- /dev/null
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -0,0 +1,84 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Internal.Common
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable (uses FFI)
+--
+-- Common tools
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Internal.Common where
+
+import Foreign
+import Complex
+import Control.Monad(when)
+import Debug.Trace
+import Data.List(transpose,intersperse)
+import Data.Typeable
+import Data.Maybe(fromJust)
+
+debug x = trace (show x) x
+
+data Vector t = V { dim  :: Int
+                  , fptr :: ForeignPtr t
+                  , ptr  :: Ptr t
+                  } deriving Typeable
+
+----------------------------------------------------------------------
+instance (Storable a, RealFloat a) => Storable (Complex a) where    --
+    alignment x = alignment (realPart x)                            --
+    sizeOf x    = 2 * sizeOf (realPart x)                           --
+    peek p = do                                                     --
+        [re,im] <- peekArray 2 (castPtr p)                          --
+        return (re :+ im)                                           --
+    poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
+----------------------------------------------------------------------
+
+on f g = \x y -> f (g x) (g y)
+
+partit :: Int -> [a] -> [[a]]
+partit _ [] = []
+partit n l  = take n l : partit n (drop n l)
+
+-- | obtains the common value of a property of a list
+common :: (Eq a) => (b->a) -> [b] -> Maybe a
+common f = commonval . map f where
+    commonval :: (Eq a) => [a] -> Maybe a
+    commonval [] = Nothing
+    commonval [a] = Just a
+    commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
+
+xor a b = a && not b || b && not a
+
+(//) :: x -> (x -> y) -> y
+infixl 0 //
+(//) = flip ($)
+
+errorCode 1000 = "bad size"
+errorCode 1001 = "bad function code"
+errorCode 1002 = "memory problem"
+errorCode 1003 = "bad file"
+errorCode 1004 = "singular"
+errorCode 1005 = "didn't converge"
+errorCode n    = "code "++show n
+
+check msg ls f = do
+    err <- f
+    when (err/=0) (error (msg++": "++errorCode err))
+    mapM_ (touchForeignPtr . fptr) ls
+    return ()
+
+class (Storable a, Typeable a) => Field a where
+instance (Storable a, Typeable a) => Field a where
+
+isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
+isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
+
+scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
+scast = fromJust . cast
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index bae56f1..2c0acdf 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Matrix
@@ -15,18 +15,16 @@
 
 module Data.Packed.Internal.Matrix where
 
+import Data.Packed.Internal.Common
 import Data.Packed.Internal.Vector
 
 import Foreign hiding (xor)
 import Complex
 import Control.Monad(when)
-import Debug.Trace
 import Data.List(transpose,intersperse)
 import Data.Typeable
 import Data.Maybe(fromJust)
 
-debug x = trace (show x) x
-
 
 data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
 
@@ -39,7 +37,7 @@ data Matrix t = M { rows    :: Int
                   , order   :: MatrixOrder
                   } deriving Typeable
 
-xor a b = a && not b || b && not a
+
 
 fortran m = order m == ColumnMajor
 
@@ -57,25 +55,12 @@ type t ::> s = Mt t s
 
 mat d m f = f (rows m) (cols m) (ptr (d m))
 
+toLists m = partit (cols m) . toList . cdat $ m
+
 instance (Show a, Storable a) => (Show (Matrix a)) where
     show m = (sizes++) . dsp . map (map show) . toLists $ m
         where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
 
-partit :: Int -> [a] -> [[a]]
-partit _ [] = []
-partit n l  = take n l : partit n (drop n l)
-
--- | obtains the common value of a property of a list
-common :: (Eq a) => (b->a) -> [b] -> Maybe a
-common f = commonval . map f where
-    commonval :: (Eq a) => [a] -> Maybe a
-    commonval [] = Nothing
-    commonval [a] = Just a
-    commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
-
-
-toLists m = partit (cols m) . toList . cdat $ m
-
 dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
     where
         mt = transpose as
@@ -146,62 +131,6 @@ transdata c1 d c2 | isReal baseOf d = scast $ transdataR c1 (scast d) c2
 --{-# RULES "transdataR" transdata=transdataR #-}
 --{-# RULES "transdataC" transdata=transdataC #-}
 
------------------------------------------------------------------------------
-
--- | creates a Matrix from a list of vectors
-fromRows :: Field t => [Vector t] -> Matrix t
-fromRows vs = case common dim vs of
-    Nothing -> error "fromRows applied to [] or to vectors with different sizes"
-    Just c  -> reshape c (join vs)
-
--- | extracts the rows of a matrix as a list of vectors
-toRows :: Storable t => Matrix t -> [Vector t]
-toRows m = toRows' 0 where
-    v = cdat m
-    r = rows m
-    c = cols m
-    toRows' k | k == r*c  = []
-              | otherwise = subVector k c v : toRows' (k+c)
-
--- | Creates a matrix from a list of vectors, as columns
-fromColumns :: Field t => [Vector t] -> Matrix t
-fromColumns m = trans . fromRows $ m
-
--- | Creates a list of vectors from the columns of a matrix
-toColumns :: Field t => Matrix t -> [Vector t]
-toColumns m = toRows . trans $ m
-
--- | creates a matrix from a vertical list of matrices
-joinVert :: Field t => [Matrix t] -> Matrix t
-joinVert ms = case common cols ms of
-    Nothing -> error "joinVert on matrices with different number of columns"
-    Just c  -> reshape c $ join (map cdat ms)
-
--- | creates a matrix from a horizontal list of matrices
-joinHoriz :: Field t => [Matrix t] -> Matrix t
-joinHoriz ms = trans. joinVert . map trans $ ms
-
--- | creates a complex vector from vectors with real and imaginary parts
-toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
-toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i]
-
--- | obtains the complex conjugate of a complex vector
-conj :: Vector (Complex Double) -> Vector (Complex Double)
-conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1])
-    where mulC = multiply RowMajor
-
-comp v = toComplex (v,constant (dim v) 0)
-
-------------------------------------------------------------------------------
-
--- | Reverse rows 
-flipud :: Field t => Matrix t -> Matrix t
-flipud m = fromRows . reverse . toRows $ m
-
--- | Reverse columns
-fliprl :: Field t => Matrix t -> Matrix t
-fliprl m = fromColumns . reverse . toColumns $ m
-
 -----------------------------------------------------------------
 
 liftMatrix f m = m { dat = f (dat m), tdat = f (tdat m) } -- check sizes
@@ -330,13 +259,25 @@ diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..
           delta i j | i==j      = 1
                     | otherwise = 0
 
-diagRect s r c
-    | dim s < min r c = error "diagRect"
-    | 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 (r*c) 0
+-- | creates a Matrix from a list of vectors
+fromRows :: Field t => [Vector t] -> Matrix t
+fromRows vs = case common dim vs of
+    Nothing -> error "fromRows applied to [] or to vectors with different sizes"
+    Just c  -> reshape c (join vs)
 
-takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
+-- | extracts the rows of a matrix as a list of vectors
+toRows :: Storable t => Matrix t -> [Vector t]
+toRows m = toRows' 0 where
+    v = cdat m
+    r = rows m
+    c = cols m
+    toRows' k | k == r*c  = []
+              | otherwise = subVector k c v : toRows' (k+c)
 
-ident n = diag (constant n 1)
+-- | Creates a matrix from a list of vectors, as columns
+fromColumns :: Field t => [Vector t] -> Matrix t
+fromColumns m = trans . fromRows $ m
+
+-- | Creates a list of vectors from the columns of a matrix
+toColumns :: Field t => Matrix t -> [Vector t]
+toColumns m = toRows . trans $ m
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
index 67dcb09..123270d 100644
--- a/lib/Data/Packed/Internal/Tensor.hs
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -1,4 +1,3 @@
---{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Tensor
@@ -9,15 +8,17 @@
 -- Stability   :  provisional
 -- Portability :  portable (uses FFI)
 --
--- Fundamental types
+-- basic tensor operations
 --
 -----------------------------------------------------------------------------
 
 module Data.Packed.Internal.Tensor where
 
+import Data.Packed.Internal
 import Data.Packed.Internal.Vector
 import Data.Packed.Internal.Matrix
 import Foreign.Storable
+import Data.List(sort)
 
 data IdxTp = Covariant | Contravariant deriving (Show,Eq)
 
@@ -99,3 +100,10 @@ compatIdxAux (n1,(t1,_)) (n2, (t2,_)) = t1 /= t2 && n1 == n2
 compatIdx t1 n1 t2 n2 = compatIdxAux d1 d2 where
     d1 = head $ snd $ fst $ findIdx n1 t1
     d2 = head $ snd $ fst $ findIdx n2 t2
+
+names t = sort $ map (snd.snd) (dims t)
+
+normal t = tridx (names t) t
+
+contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
+
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 4836bdb..125df1e 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Internal.Vector
@@ -9,70 +9,16 @@
 -- Stability   :  provisional
 -- Portability :  portable (uses FFI)
 --
--- Fundamental types
+-- Vector implementation
 --
 -----------------------------------------------------------------------------
 
 module Data.Packed.Internal.Vector where
 
+import Data.Packed.Internal.Common
 import Foreign
 import Complex
 import Control.Monad(when)
-import Debug.Trace
-import Data.List(transpose,intersperse)
-import Data.Typeable
-import Data.Maybe(fromJust)
-
-debug x = trace (show x) x
-
-----------------------------------------------------------------------
-instance (Storable a, RealFloat a) => Storable (Complex a) where    --
-    alignment x = alignment (realPart x)                            --
-    sizeOf x    = 2 * sizeOf (realPart x)                           --
-    peek p = do                                                     --
-        [re,im] <- peekArray 2 (castPtr p)                          --
-        return (re :+ im)                                           --
-    poke p (a :+ b) = pokeArray (castPtr p) [a,b]                   --
-----------------------------------------------------------------------
-
-on f g = \x y -> f (g x) (g y)
-
-(//) :: x -> (x -> y) -> y
-infixl 0 //
-(//) = flip ($)
-
-errorCode 1000 = "bad size"
-errorCode 1001 = "bad function code"
-errorCode 1002 = "memory problem"
-errorCode 1003 = "bad file"
-errorCode 1004 = "singular"
-errorCode 1005 = "didn't converge"
-errorCode n    = "code "++show n
-
-check msg ls f = do
-    err <- f
-    when (err/=0) (error (msg++": "++errorCode err))
-    mapM_ (touchForeignPtr . fptr) ls
-    return ()
-
-class (Storable a, Typeable a) => Field a where
-instance (Storable a, Typeable a) => Field a where
-
-isReal w x   = typeOf (undefined :: Double) == typeOf (w x)
-isComp w x = typeOf (undefined :: Complex Double) == typeOf (w x)
-baseOf v = (v `at` 0)
-
-scast :: forall a . forall b . (Typeable a, Typeable b) => a -> b
-scast = fromJust . cast
-
-
-
-----------------------------------------------------------------------
-
-data Vector t = V { dim  :: Int
-                  , fptr :: ForeignPtr t
-                  , ptr  :: Ptr t
-                  } deriving Typeable
 
 type Vc t s = Int -> Ptr t -> s
 infixr 5 :>
@@ -81,6 +27,8 @@ type t :> s = Vc t s
 vec :: Vector t -> (Vc t s) -> s
 vec v f = f (dim v) (ptr v)
 
+baseOf v = (v `at` 0)
+
 createVector :: Storable a => Int -> IO (Vector a)
 createVector n = do
     when (n <= 0) $ error ("trying to createVector of dim "++show n)
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 8d1c8b6..6ea76c0 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -1 +1,71 @@
- 
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Matrices
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Matrix (
+    Matrix(rows,cols), Field,
+    toLists, (><), (>|<),
+    trans,
+    reshape,
+    fromRows, toRows, fromColumns, toColumns,
+    joinVert, joinHoriz,
+    flipud, fliprl,
+    liftMatrix, liftMatrix2,
+    multiply,
+    subMatrix, diag, takeDiag, diagRect, ident
+) where
+
+import Data.Packed.Internal
+
+-- | creates a matrix from a vertical list of matrices
+joinVert :: Field t => [Matrix t] -> Matrix t
+joinVert ms = case common cols ms of
+    Nothing -> error "joinVert on matrices with different number of columns"
+    Just c  -> reshape c $ join (map cdat ms)
+
+-- | creates a matrix from a horizontal list of matrices
+joinHoriz :: Field t => [Matrix t] -> Matrix t
+joinHoriz ms = trans. joinVert . map trans $ ms
+
+-- | Reverse rows 
+flipud :: Field t => Matrix t -> Matrix t
+flipud m = fromRows . reverse . toRows $ m
+
+-- | Reverse columns
+fliprl :: Field t => Matrix t -> Matrix t
+fliprl m = fromColumns . reverse . toColumns $ m
+
+------------------------------------------------------------
+
+diagRect s r c
+    | dim s < min r c = error "diagRect"
+    | 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 (r*c) 0
+
+takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
+
+ident n = diag (constant n 1)
+
+r >< c = f where
+    f l | dim v == r*c = matrixFromVector RowMajor c v
+        | otherwise    = error $ "inconsistent list size = "
+                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+        where v = fromList l
+
+r >|< c = f where
+    f l | dim v == r*c = matrixFromVector ColumnMajor c v
+        | otherwise    = error $ "inconsistent list size = "
+                                 ++show (dim v) ++"in ("++show r++"><"++show c++")"
+        where v = fromList l
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 8d1c8b6..992301a 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -1 +1,40 @@
- 
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Vectors
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Vector (
+    Vector(dim), Field,
+    fromList, toList,
+    at,
+    subVector, join,
+    constant,
+    toComplex, comp,
+    conj,
+    dot
+) where
+
+import Data.Packed.Internal
+import Complex
+
+-- | creates a complex vector from vectors with real and imaginary parts
+toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)
+toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i]
+
+-- | obtains the complex conjugate of a complex vector
+conj :: Vector (Complex Double) -> Vector (Complex Double)
+conj v = asComplex $ cdat $ reshape 2 (asReal v) `mulC` diag (fromList [1,-1])
+    where mulC = multiply RowMajor
+
+comp v = toComplex (v,constant (dim v) 0)
+
+