refactor Numeric Vector/Matrix and classes

13 files changed, 682 insertions, 292 deletions

diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
index 50a5eb6..8128667 100644
--- a/lib/Data/Packed.hs
+++ b/lib/Data/Packed.hs
@@ -24,4 +24,3 @@ import Data.Packed.Vector
 import Data.Packed.Matrix
 import Data.Packed.Random
 import Data.Complex
-import Numeric.LinearAlgebra.Instances()
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 2855a90..07258f8 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -19,10 +19,7 @@
 -----------------------------------------------------------------------------
 
 module Data.Packed.Matrix (
-    Element, RealElement, Container(..),
-    Convert(..), RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf,
-    Precision(..), comp,
-    AutoReal(..),
+    Element(..),
     Matrix,rows,cols,
     (><),
     trans,
@@ -55,7 +52,7 @@ import Data.Complex
 import Data.Binary
 import Foreign.Storable
 import Control.Monad(replicateM)
-import Control.Arrow((***))
+--import Control.Arrow((***))
 --import GHC.Float(double2Float,float2Double)
 
 -------------------------------------------------------------------
@@ -75,6 +72,33 @@ instance (Binary a, Element a, Storable a) => Binary (Matrix a) where
 
 -------------------------------------------------------------------
 
+instance (Show a, Element a) => (Show (Matrix a)) where
+    show m = (sizes++) . dsp . map (map show) . toLists $ m
+        where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
+
+dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
+    where
+        mt = transpose as
+        longs = map (maximum . map length) mt
+        mtp = zipWith (\a b -> map (pad a) b) longs mt
+        pad n str = replicate (n - length str) ' ' ++ str
+        unwords' = concat . intersperse ", "
+
+------------------------------------------------------------------
+
+instance (Element a, Read a) => Read (Matrix a) where
+    readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
+        where (thing,rest) = breakAt ']' s
+              (dims,listnums) = breakAt ')' thing
+              cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
+              rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
+
+
+breakAt c l = (a++[c],tail b) where
+    (a,b) = break (==c) l
+
+------------------------------------------------------------------
+
 -- | creates a matrix from a vertical list of matrices
 joinVert :: Element t => [Matrix t] -> Matrix t
 joinVert ms = case common cols ms of
@@ -470,150 +494,3 @@ toBlocksEvery r c m = toBlocks rs cs m where
     cs = replicate qc c ++ if rc > 0 then [rc] else []
 
 -------------------------------------------------------------------
-
--- | Supported single-double precision type pairs
-class (Element s, Element d) => Precision s d | s -> d, d -> s where
-    double2FloatG :: Vector d -> Vector s
-    float2DoubleG :: Vector s -> Vector d
-
-instance Precision Float Double where
-    double2FloatG = double2FloatV
-    float2DoubleG = float2DoubleV
-
-instance Precision (Complex Float) (Complex Double) where
-    double2FloatG = asComplex . double2FloatV . asReal
-    float2DoubleG = asComplex . float2DoubleV . asReal
-
--- | Supported real types
-class (Element t, Element (Complex t), RealFloat t
---       , RealOf t ~ t, RealOf (Complex t) ~ t
-       )
-       => RealElement t
-
-instance RealElement Double
-
-instance RealElement Float
-
--- | Conversion utilities
-class Container c where
-    toComplex   :: (RealElement e) => (c e, c e) -> c (Complex e)
-    fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e)
-    complex'    :: (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'      :: Precision a b => c b -> c a
-    double'      :: Precision a b => c a -> c b
-
-comp x = complex' x
-
-instance Container Vector where
-    toComplex = toComplexV
-    fromComplex = fromComplexV
-    complex' v = toComplex (v,constantD 0 (dim v))
-    conj = conjV
-    cmap = mapVector
-    single' = double2FloatG
-    double' = float2DoubleG
-
-instance Container Matrix where
-    toComplex = uncurry $ liftMatrix2 $ curry toComplex
-    fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z
-        where c = cols z
-    complex' = liftMatrix complex'
-    conj = liftMatrix conj
-    cmap f = liftMatrix (cmap f)
-    single' = liftMatrix single'
-    double' = liftMatrix double'
-
--------------------------------------------------------------------
-
-type family RealOf x
-
-type instance RealOf Double = Double
-type instance RealOf (Complex Double) = Double
-
-type instance RealOf Float = Float
-type instance RealOf (Complex Float) = Float
-
-type family ComplexOf x
-
-type instance ComplexOf Double = Complex Double
-type instance ComplexOf (Complex Double) = Complex Double
-
-type instance ComplexOf Float = Complex Float
-type instance ComplexOf (Complex Float) = Complex Float
-
-type family SingleOf x
-
-type instance SingleOf Double = Float
-type instance SingleOf Float  = Float
-
-type instance SingleOf (Complex a) = Complex (SingleOf a)
-
-type family DoubleOf x
-
-type instance DoubleOf Double = Double
-type instance DoubleOf Float  = Double
-
-type instance DoubleOf (Complex a) = Complex (DoubleOf a)
-
-type family ElementOf c
-
-type instance ElementOf (Vector a) = a
-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)
-
-instance Convert Double where
-    real = id
-    complex = complex'
-    single = single'
-    double = id
-
-instance Convert Float where
-    real = id
-    complex = complex'
-    single = id
-    double = double'
-
-instance Convert (Complex Double) where
-    real = complex'
-    complex = id
-    single = single'
-    double = id
-
-instance Convert (Complex Float) where
-    real = complex'
-    complex = id
-    single = id
-    double = double'
-
--------------------------------------------------------------------
-
--- | to be replaced by Convert
-class Convert t => AutoReal t where
-    real'' :: Container c => c Double -> c t
-    complex'' :: Container c => c t -> c (Complex Double)
-
-instance AutoReal Double where
-    real'' = real
-    complex'' = complex
-
-instance AutoReal (Complex Double) where
-    real'' = real
-    complex'' = complex
-
-instance AutoReal Float where
-    real'' = real . single
-    complex'' = double . complex
-
-instance AutoReal (Complex Float) where
-    real'' = real . single
-    complex'' = double . complex
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
index 33a11d7..579d13c 100644
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -22,8 +22,8 @@ module Data.Packed.Random (
 import Numeric.GSL.Vector
 import Data.Packed.Matrix
 import Data.Packed.Vector
-import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Linear
+import Numeric.LinearAlgebra.Algorithms
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- Gaussian distribution.
diff --git a/lib/Graphics/Plot.hs b/lib/Graphics/Plot.hs
index 2dc0553..478b4ad 100644
--- a/lib/Graphics/Plot.hs
+++ b/lib/Graphics/Plot.hs
@@ -150,8 +150,8 @@ matrixToPGM m = header ++ unlines (map unwords ll) where
     r = rows m
     header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"
     maxgray = 255.0
-    maxval = vectorMax $ flatten $ m
-    minval = vectorMin $ flatten $ m
+    maxval = maxElement m
+    minval = minElement m
     scale' = if (maxval == minval) 
         then 0.0
         else maxgray / (maxval - minval)
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
new file mode 100644
index 0000000..0bec2e8
--- /dev/null
+++ b/lib/Numeric/Container.hs
@@ -0,0 +1,219 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Container
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Numeric classes for containers of numbers, including conversion routines
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Container (
+    RealElement, AutoReal(..),
+    Container(..), Linear(..),
+    Convert(..), RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf, IndexOf,
+    Precision(..), comp,
+    module Data.Complex
+) where
+
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
+--import qualified Data.Packed.ST as ST
+
+import Control.Arrow((***))
+
+import Data.Complex
+
+-------------------------------------------------------------------
+
+-- | Supported single-double precision type pairs
+class (Element s, Element d) => Precision s d | s -> d, d -> s where
+    double2FloatG :: Vector d -> Vector s
+    float2DoubleG :: Vector s -> Vector d
+
+instance Precision Float Double where
+    double2FloatG = double2FloatV
+    float2DoubleG = float2DoubleV
+
+instance Precision (Complex Float) (Complex Double) where
+    double2FloatG = asComplex . double2FloatV . asReal
+    float2DoubleG = asComplex . float2DoubleV . asReal
+
+-- | Supported real types
+class (Element t, Element (Complex t), RealFloat t
+--       , RealOf t ~ t, RealOf (Complex t) ~ t
+       )
+    => RealElement t
+
+instance RealElement Double
+
+instance RealElement Float
+
+-- | Conversion utilities
+class Container c where
+    toComplex   :: (RealElement e) => (c e, c e) -> c (Complex e)
+    fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e)
+    complex'    :: (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'      :: Precision a b => c b -> c a
+    double'      :: Precision a b => c a -> c b
+
+comp x = complex' x
+
+instance Container Vector where
+    toComplex = toComplexV
+    fromComplex = fromComplexV
+    complex' v = toComplex (v,constantD 0 (dim v))
+    conj = conjV
+    cmap = mapVector
+    single' = double2FloatG
+    double' = float2DoubleG
+
+instance Container Matrix where
+    toComplex = uncurry $ liftMatrix2 $ curry toComplex
+    fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z
+        where c = cols z
+    complex' = liftMatrix complex'
+    conj = liftMatrix conj
+    cmap f = liftMatrix (cmap f)
+    single' = liftMatrix single'
+    double' = liftMatrix double'
+
+-------------------------------------------------------------------
+
+type family RealOf x
+
+type instance RealOf Double = Double
+type instance RealOf (Complex Double) = Double
+
+type instance RealOf Float = Float
+type instance RealOf (Complex Float) = Float
+
+type family ComplexOf x
+
+type instance ComplexOf Double = Complex Double
+type instance ComplexOf (Complex Double) = Complex Double
+
+type instance ComplexOf Float = Complex Float
+type instance ComplexOf (Complex Float) = Complex Float
+
+type family SingleOf x
+
+type instance SingleOf Double = Float
+type instance SingleOf Float  = Float
+
+type instance SingleOf (Complex a) = Complex (SingleOf a)
+
+type family DoubleOf x
+
+type instance DoubleOf Double = Double
+type instance DoubleOf Float  = Double
+
+type instance DoubleOf (Complex a) = Complex (DoubleOf a)
+
+type family ElementOf c
+
+type instance ElementOf (Vector a) = a
+type instance ElementOf (Matrix a) = a
+
+type family IndexOf c
+
+type instance IndexOf Vector = Int
+type instance IndexOf Matrix = (Int,Int)
+
+-------------------------------------------------------------------
+
+-- | 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)
+
+instance Convert Double where
+    real = id
+    complex = complex'
+    single = single'
+    double = id
+
+instance Convert Float where
+    real = id
+    complex = complex'
+    single = id
+    double = double'
+
+instance Convert (Complex Double) where
+    real = complex'
+    complex = id
+    single = single'
+    double = id
+
+instance Convert (Complex Float) where
+    real = complex'
+    complex = id
+    single = id
+    double = double'
+
+-------------------------------------------------------------------
+
+-- | to be replaced by Convert
+class Convert t => AutoReal t where
+    real'' :: Container c => c Double -> c t
+    complex'' :: Container c => c t -> c (Complex Double)
+
+instance AutoReal Double where
+    real'' = real
+    complex'' = complex
+
+instance AutoReal (Complex Double) where
+    real'' = real
+    complex'' = complex
+
+instance AutoReal Float where
+    real'' = real . single
+    complex'' = double . complex
+
+instance AutoReal (Complex Float) where
+    real'' = real . single
+    complex'' = double . complex
+
+-------------------------------------------------------------------
+
+-- | Basic element-by-element functions.
+class (Element e, Container c) => Linear c e where
+    -- | create a structure with a single element
+    scalar      :: e -> c e
+    scale       :: e -> c e -> c e
+    -- | scale the element by element reciprocal of the object:
+    --
+    -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
+    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
+    -- | element by element multiplication
+    mul         :: c e -> c e -> c e
+    -- | element by element division
+    divide      :: c e -> c e -> c e
+    equal       :: c e -> c e -> Bool
+    --
+    minIndex    :: c e -> IndexOf c
+    maxIndex    :: c e -> IndexOf c
+    minElement  :: c e -> e
+    maxElement  :: c e -> e
+
+
+
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index e8a14d6..3df9bd7 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -13,13 +13,11 @@ This module reexports all normally required functions for Linear Algebra applica
 -}
 -----------------------------------------------------------------------------
 module Numeric.LinearAlgebra (
-    module Data.Packed,
     module Numeric.LinearAlgebra.Algorithms,
     module Numeric.LinearAlgebra.Interface,
     module Numeric.LinearAlgebra.Linear
 ) where
 
-import Data.Packed
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Interface
 import Numeric.LinearAlgebra.Linear
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 7d2f84d..14bf5d8 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -77,8 +77,8 @@ module Numeric.LinearAlgebra.Algorithms (
 import Data.Packed.Internal hiding ((//))
 import Data.Packed.Matrix
 import Data.Complex
-import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Numeric.LinearAlgebra.Linear
+import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Data.List(foldl1')
 import Data.Array
 
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
index 13d175b..fa3e209 100644
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -25,8 +25,8 @@ module Numeric.LinearAlgebra.Interface(
     (<|>),(<->),
 ) where
 
-import Data.Packed.Vector
-import Data.Packed.Matrix
+import Numeric.Vector
+import Numeric.Matrix
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Linear
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index eec3035..8888712 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -44,6 +44,7 @@ module Numeric.LinearAlgebra.LAPACK (
 import Data.Packed.Internal
 import Data.Packed.Matrix
 import Data.Complex
+import Numeric.Container
 import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
 import Foreign
 import Foreign.C.Types (CInt)
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 778b976..9048204 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -19,33 +19,31 @@ Basic optimized operations on vectors and matrices.
 module Numeric.LinearAlgebra.Linear (
     -- * Linear Algebra Typeclasses
     Vectors(..),
-    Linear(..),
     -- * Products
     Product(..),
     mXm,mXv,vXm,
     outer, kronecker,
-    -- * Creation of numeric vectors
-    constant, linspace
+    -- * Modules
+    module Numeric.Vector,
+    module Numeric.Matrix,
+    module Numeric.Container
 ) where
 
-import Data.Packed.Internal
+import Data.Packed.Internal.Common
 import Data.Packed.Matrix
 import Data.Complex
+import Numeric.Container
+import Numeric.Vector
+import Numeric.Matrix
 import Numeric.GSL.Vector
 import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
 
-import Control.Monad(ap)
-
 -- | basic Vector functions
 class Num e => Vectors a e where
     -- the C functions sumX are twice as fast as using foldVector
     vectorSum  :: a e -> e
     vectorProd :: a e -> e
     absSum     :: a e -> e
-    vectorMin  :: a e -> e
-    vectorMax  :: a e -> e
-    minIdx     :: a e -> Int
-    maxIdx     :: a e -> Int
     dot        :: a e -> a e -> e
     norm1      :: a e -> e
     norm2      :: a e -> e
@@ -57,153 +55,36 @@ instance Vectors Vector Float where
     vectorProd = prodF
     norm2      = toScalarF Norm2
     absSum     = toScalarF AbsSum
-    vectorMin  = toScalarF Min
-    vectorMax  = toScalarF Max
-    minIdx     = round . toScalarF MinIdx
-    maxIdx     = round . toScalarF MaxIdx
     dot        = dotF
     norm1      = toScalarF AbsSum
-    normInf    = vectorMax . vectorMapF Abs
+    normInf    = maxElement . vectorMapF Abs
 
 instance Vectors Vector Double where
     vectorSum  = sumR
     vectorProd = prodR
     norm2      = toScalarR Norm2
     absSum     = toScalarR AbsSum
-    vectorMin  = toScalarR Min
-    vectorMax  = toScalarR Max
-    minIdx     = round . toScalarR MinIdx
-    maxIdx     = round . toScalarR MaxIdx
     dot        = dotR
     norm1      = toScalarR AbsSum
-    normInf    = vectorMax . vectorMapR Abs
+    normInf    = maxElement . vectorMapR Abs
 
 instance Vectors Vector (Complex Float) where
     vectorSum  = sumQ
     vectorProd = prodQ
     norm2      = (:+ 0) . toScalarQ Norm2
     absSum     = (:+ 0) . toScalarQ AbsSum
-    vectorMin  = ap (@>) minIdx
-    vectorMax  = ap (@>) maxIdx
-    minIdx     = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIdx     = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     dot        = dotQ
     norm1      = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
-    normInf    = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
+    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs
 
 instance Vectors Vector (Complex Double) where
     vectorSum  = sumC
     vectorProd = prodC
     norm2      = (:+ 0) . toScalarC Norm2
     absSum     = (:+ 0) . toScalarC AbsSum
-    vectorMin  = ap (@>) minIdx
-    vectorMax  = ap (@>) maxIdx
-    minIdx     = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIdx     = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
     dot        = dotC
     norm1      = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
-    normInf    = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
-
-----------------------------------------------------
-
--- | Basic element-by-element functions.
-class (Element e, Container c) => Linear c e where
-    -- | create a structure with a single element
-    scalar      :: e -> c e
-    scale       :: e -> c e -> c e
-    -- | scale the element by element reciprocal of the object:
-    --
-    -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
-    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
-    -- | element by element multiplication
-    mul         :: c e -> c e -> c e
-    -- | element by element division
-    divide      :: c e -> c e -> c e
-    equal       :: c e -> c e -> Bool
-
-
-instance Linear Vector Float where
-    scale = vectorMapValF Scale
-    scaleRecip = vectorMapValF Recip
-    addConstant = vectorMapValF AddConstant
-    add = vectorZipF Add
-    sub = vectorZipF Sub
-    mul = vectorZipF Mul
-    divide = vectorZipF Div
-    equal u v = dim u == dim v && vectorMax (vectorMapF Abs (sub u v)) == 0.0
-    scalar x = fromList [x]
-
-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
-    equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0
-    scalar x = fromList [x]
-
-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
-    equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0
-    scalar x = fromList [x]
-
-instance Linear Vector (Complex Float) where
-    scale = vectorMapValQ Scale
-    scaleRecip = vectorMapValQ Recip
-    addConstant = vectorMapValQ AddConstant
-    add = vectorZipQ Add
-    sub = vectorZipQ Sub
-    mul = vectorZipQ Mul
-    divide = vectorZipQ Div
-    equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0
-    scalar x = fromList [x]
-
-instance (Linear Vector a, Container Matrix) => (Linear Matrix a) 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
-    equal a b = cols a == cols b && flatten a `equal` flatten b
-    scalar x = (1><1) [x]
-
-
-----------------------------------------------------
-
-{- | creates a vector with a given number of equal components:
-
-@> constant 2 7
-7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
--}
-constant :: Element a => a -> Int -> Vector a
--- constant x n = runSTVector (newVector x n)
-constant = constantD -- about 2x faster
-
-{- | Creates a real vector containing a range of values:
-
-@\> linspace 5 (-3,7)
-5 |> [-3.0,-0.5,2.0,4.5,7.0]@
-
-Logarithmic spacing can be defined as follows:
-
-@logspace n (a,b) = 10 ** linspace n (a,b)@
--}
-linspace :: (Enum e, Linear Vector e) => Int -> (e, e) -> Vector e
-linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
-    where s = (b-a)/fromIntegral (n-1)
+    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs
 
 ----------------------------------------------------
 
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 43a62e5..5b42226 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -21,6 +21,7 @@ module Numeric.LinearAlgebra.Tests(
 --, runBigTests
 ) where
 
+import Data.Packed.Random
 import Numeric.LinearAlgebra
 import Numeric.LinearAlgebra.LAPACK
 import Numeric.LinearAlgebra.Tests.Instances
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
new file mode 100644
index 0000000..8d3764a
--- /dev/null
+++ b/lib/Numeric/Matrix.hs
@@ -0,0 +1,97 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Numeric instances and functions for 'Data.Packed.Matrix's
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Matrix (
+                       module Data.Packed.Matrix,
+                      ) where
+
+-------------------------------------------------------------------
+
+import Data.Packed.Vector
+import Data.Packed.Matrix
+import Numeric.Container
+import Numeric.Vector()
+
+import Control.Monad(ap)
+
+-------------------------------------------------------------------
+
+instance Linear Matrix a => Eq (Matrix a) where
+    (==) = equal
+
+instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
+    (+) = liftMatrix2Auto (+)
+    (-) = liftMatrix2Auto (-)
+    negate = liftMatrix negate
+    (*) = liftMatrix2Auto (*)
+    signum = liftMatrix signum
+    abs = liftMatrix abs
+    fromInteger = (1><1) . return . fromInteger
+
+---------------------------------------------------
+
+instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
+    fromRational n = (1><1) [fromRational n]
+    (/) = liftMatrix2Auto (/)
+
+---------------------------------------------------------
+
+instance (Linear Vector a, Floating (Vector a), Fractional (Matrix 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
+    (**)  = liftMatrix2Auto (**)
+    sqrt  = liftMatrix sqrt
+    pi    = (1><1) [pi]
+
+---------------------------------------------------------------
+
+instance (Linear Vector a, Container Matrix) => (Linear Matrix a) 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
+    equal a b = cols a == cols b && flatten a `equal` flatten b
+    scalar x = (1><1) [x]
+    minIndex m = let (r,c) = (rows m,cols m)
+                     i = 1 + (minIndex $ flatten m)
+                 in (i `div` r,i `mod` r)
+    maxIndex m = let (r,c) = (rows m,cols m)
+                     i = 1 + (maxIndex $ flatten m)
+                 in (i `div` r,i `mod` r)
+    minElement = ap (@@>) minIndex
+    maxElement = ap (@@>) maxIndex
+
+----------------------------------------------------
+
diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs
new file mode 100644
index 0000000..ced202f
--- /dev/null
+++ b/lib/Numeric/Vector.hs
@@ -0,0 +1,317 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+--{-# LANGUAGE FunctionalDependencies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Numeric instances and functions for 'Data.Packed.Vector's
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Vector (
+                       -- * Vector creation
+                       constant, linspace,
+                       module Data.Packed.Vector
+                      ) where
+
+import Data.Complex
+
+import Control.Monad(ap)
+
+import Data.Packed.Vector
+import Data.Packed.Matrix(Element(..))
+import Numeric.GSL.Vector
+
+import Numeric.Container
+
+-------------------------------------------------------------------
+
+#ifndef VECTOR
+import Foreign(Storable)
+#endif
+
+------------------------------------------------------------------
+
+#ifndef VECTOR
+
+instance (Show a, Storable a) => (Show (Vector a)) where
+    show v = (show (dim v))++" |> " ++ show (toList v)
+
+#endif
+
+#ifdef VECTOR
+
+instance (Element a, Read a) => Read (Vector a) where
+    readsPrec _ s = [(fromList . read $ listnums, rest)]
+        where (thing,trest) = breakAt ']' s
+              (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)
+              rest = drop 31 trest
+#else
+
+instance (Element a, Read a) => Read (Vector a) where
+    readsPrec _ s = [((d |>) . read $ listnums, rest)]
+        where (thing,rest) = breakAt ']' s
+              (dims,listnums) = breakAt '>' thing
+              d = read . init . fst . breakAt '|' $ dims
+
+#endif
+
+breakAt c l = (a++[c],tail b) where
+    (a,b) = break (==c) l
+
+------------------------------------------------------------------
+
+{- | creates a vector with a given number of equal components:
+
+@> constant 2 7
+7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
+-}
+constant :: Element a => a -> Int -> Vector a
+-- constant x n = runSTVector (newVector x n)
+constant = constantD -- about 2x faster
+
+{- | Creates a real vector containing a range of values:
+
+@\> linspace 5 (-3,7)
+5 |> [-3.0,-0.5,2.0,4.5,7.0]@
+
+Logarithmic spacing can be defined as follows:
+
+@logspace n (a,b) = 10 ** linspace n (a,b)@
+-}
+linspace :: (Enum e, Linear Vector e) => Int -> (e, e) -> Vector e
+linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
+    where s = (b-a)/fromIntegral (n-1)
+
+------------------------------------------------------------------
+
+adaptScalar f1 f2 f3 x y
+    | dim x == 1 = f1   (x@>0) y
+    | dim y == 1 = f3 x (y@>0)
+    | otherwise = f2 x y
+
+------------------------------------------------------------------
+
+#ifndef VECTOR
+
+instance Linear Vector a => Eq (Vector a) where
+    (==) = equal
+
+#endif
+
+instance Num (Vector Float) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapF Sign
+    abs = vectorMapF Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector Double) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapR Sign
+    abs = vectorMapR Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector (Complex Double)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapC Sign
+    abs = vectorMapC Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector (Complex Float)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapQ Sign
+    abs = vectorMapQ Abs
+    fromInteger = fromList . return . fromInteger
+
+---------------------------------------------------
+
+instance (Linear Vector a, Num (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 Floating (Vector Float) where
+    sin   = vectorMapF Sin
+    cos   = vectorMapF Cos
+    tan   = vectorMapF Tan
+    asin  = vectorMapF ASin
+    acos  = vectorMapF ACos
+    atan  = vectorMapF ATan
+    sinh  = vectorMapF Sinh
+    cosh  = vectorMapF Cosh
+    tanh  = vectorMapF Tanh
+    asinh = vectorMapF ASinh
+    acosh = vectorMapF ACosh
+    atanh = vectorMapF ATanh
+    exp   = vectorMapF Exp
+    log   = vectorMapF Log
+    sqrt  = vectorMapF Sqrt
+    (**)  = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))
+    pi    = fromList [pi]
+
+-------------------------------------------------------------
+
+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 Floating (Vector (Complex Float)) where
+    sin   = vectorMapQ Sin
+    cos   = vectorMapQ Cos
+    tan   = vectorMapQ Tan
+    asin  = vectorMapQ ASin
+    acos  = vectorMapQ ACos
+    atan  = vectorMapQ ATan
+    sinh  = vectorMapQ Sinh
+    cosh  = vectorMapQ Cosh
+    tanh  = vectorMapQ Tanh
+    asinh = vectorMapQ ASinh
+    acosh = vectorMapQ ACosh
+    atanh = vectorMapQ ATanh
+    exp   = vectorMapQ Exp
+    log   = vectorMapQ Log
+    sqrt  = vectorMapQ Sqrt
+    (**)  = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
+    pi    = fromList [pi]
+
+-----------------------------------------------------------
+
+
+-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
+--     mempty = 0 { idim = 0 }
+--     mappend a b = mconcat [a,b]
+--     mconcat = j . filter ((>0).dim)
+--         where j [] = mempty
+--               j l  = join l
+
+---------------------------------------------------------------
+
+-- instance (NFData a, Storable a) => NFData (Vector a) where
+--     rnf = rnf . (@>0)
+--
+-- instance (NFData a, Element a) => NFData (Matrix a) where
+--     rnf = rnf . flatten
+
+---------------------------------------------------------------
+
+instance Linear Vector Float where
+    scale = vectorMapValF Scale
+    scaleRecip = vectorMapValF Recip
+    addConstant = vectorMapValF AddConstant
+    add = vectorZipF Add
+    sub = vectorZipF Sub
+    mul = vectorZipF Mul
+    divide = vectorZipF Div
+    equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
+    scalar x = fromList [x]
+    minIndex     = round . toScalarF MinIdx
+    maxIndex     = round . toScalarF MaxIdx
+    minElement  = toScalarF Min
+    maxElement  = toScalarF Max
+
+
+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
+    equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
+    scalar x = fromList [x]
+    minIndex     = round . toScalarR MinIdx
+    maxIndex     = round . toScalarR MaxIdx
+    minElement  = toScalarR Min
+    maxElement  = toScalarR Max
+
+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
+    equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
+    scalar x = fromList [x]
+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    minElement  = ap (@>) minIndex
+    maxElement  = ap (@>) maxIndex
+
+instance Linear Vector (Complex Float) where
+    scale = vectorMapValQ Scale
+    scaleRecip = vectorMapValQ Recip
+    addConstant = vectorMapValQ AddConstant
+    add = vectorZipQ Add
+    sub = vectorZipQ Sub
+    mul = vectorZipQ Mul
+    divide = vectorZipQ Div
+    equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
+    scalar x = fromList [x]
+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    minElement  = ap (@>) minIndex
+    maxElement  = ap (@>) maxIndex
+
+---------------------------------------------------------------