re-export changes

12 files changed, 403 insertions, 910 deletions

diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
index 8128667..bfc2d8b 100644
--- a/lib/Data/Packed.hs
+++ b/lib/Data/Packed.hs
@@ -16,11 +16,13 @@ The Vector and Matrix types and some utilities.
 module Data.Packed (
     module Data.Packed.Vector,
     module Data.Packed.Matrix,
-    module Data.Packed.Random,
-    module Data.Complex
+--    module Numeric.Conversion,
+--    module Data.Packed.Random,
+--    module Data.Complex
 ) where
 
 import Data.Packed.Vector
 import Data.Packed.Matrix
-import Data.Packed.Random
-import Data.Complex
+--import Data.Packed.Random
+--import Data.Complex
+--import Numeric.Conversion
\ No newline at end of file
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
index 0bc5350..b30b299 100644
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -20,11 +20,12 @@ module Data.Packed.Random (
 ) where
 
 import Numeric.GSL.Vector
-import Data.Packed.Matrix
-import Numeric.Vector
-import Numeric.LinearAlgebra.Linear
+import Data.Packed
+import Numeric.Container
+import Data.Packed.Internal(constantD)
 import Numeric.LinearAlgebra.Algorithms
 
+constant k v = constantD k v
 
 -- | 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 ae88a92..c5b5a4c 100644
--- a/lib/Graphics/Plot.hs
+++ b/lib/Graphics/Plot.hs
@@ -28,9 +28,7 @@ module Graphics.Plot(
 
 ) where
 
-import Data.Packed
-import Numeric.Vector
-import Numeric.LinearAlgebra.Linear
+import Numeric.Matrix
 import Data.List(intersperse)
 import System.Process (system)
 
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index d6b1342..fc4fb0d 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE UndecidableInstances #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -19,103 +20,35 @@
 -----------------------------------------------------------------------------
 
 module Numeric.Container (
-    Linear(..),
-    Container(..), RealElement, Precision(..), NumericContainer(..), comp, 
-    Convert(..), --AutoReal(..), 
-    RealOf, ComplexOf, SingleOf, DoubleOf, 
-
---    ElementOf, 
+    --Linear(..),
+    Container(..),
+    Vectors(..),
+    Product(..),
+    mXm,mXv,vXm,
+    outer, kronecker,
+
+    RealElement, --Precision,
+    ComplexContainer(toComplex,fromComplex,comp,conj),
+    Convert(..), --AutoReal(..),
+    RealOf, ComplexOf, SingleOf, DoubleOf,
 
     IndexOf,
-
     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 Data.Packed
+import Numeric.Conversion
+import Data.Packed.Internal
+import Numeric.GSL.Vector
 
---import Control.Arrow((***))
 import Data.Complex
+import Control.Monad(ap)
+--import Control.Arrow((***))
 
--------------------------------------------------------------------
-
--- | 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 NumericContainer 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
-
--- | a synonym for "complex'"
-comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e)
-comp x = complex' x
+import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
 
 -------------------------------------------------------------------
 
-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
@@ -123,75 +56,16 @@ type instance IndexOf Matrix = (Int,Int)
 
 -------------------------------------------------------------------
 
-class (Element t, Element (RealOf t)) => Convert t where
-    real    :: NumericContainer c => c (RealOf t) -> c t
-    complex :: NumericContainer c => c t -> c (ComplexOf t)
-    single  :: NumericContainer c => c t -> c (SingleOf t)
-    double  :: NumericContainer 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'' :: NumericContainer c => c Double -> c t
-    complex'' :: NumericContainer 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 for numeric containers
 class (Element e) => Container c e where
-{-
+
     -- | create a structure with a single element
     scalar      :: e -> c e
-    -- | multiply every element by a scalar
     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
-    -- | add a constant to each element
     addConstant :: e -> c e -> c e
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
@@ -200,7 +74,7 @@ class (Element e) => Container c e where
     -- | element by element division
     divide      :: c e -> c e -> c e
     equal       :: c e -> c e -> Bool
--}
+
     -- | cannot implement instance Functor because of Element class constraint
     cmap        :: (Element a, Element b) => (a -> b) -> c a -> c b
     --
@@ -217,23 +91,239 @@ class (Element e) => Container c e where
     sumElements :: c e -> e
     prodElements :: c e -> e
 
--- | Basic element-by-element functions.
-class (Element e, Container c e) => 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
+-- -- | Basic element-by-element functions.
+-- class (Element e, Container c e) => Linear c e where
 
 
+--------------------------------------------------------------------------
 
+instance Container 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]
+    --
+--instance Container Vector Float where
+    cmap = mapVector
+    atIndex = (@>)
+    minIndex     = round . toScalarF MinIdx
+    maxIndex     = round . toScalarF MaxIdx
+    minElement  = toScalarF Min
+    maxElement  = toScalarF Max
+    sumElements  = sumF
+    prodElements = prodF
+
+instance Container 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]
+    --
+--instance Container Vector Double where
+    cmap = mapVector
+    atIndex = (@>)
+    minIndex     = round . toScalarR MinIdx
+    maxIndex     = round . toScalarR MaxIdx
+    minElement  = toScalarR Min
+    maxElement  = toScalarR Max
+    sumElements  = sumR
+    prodElements = prodR
+
+instance Container 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]
+    --
+--instance Container Vector (Complex Double) where
+    cmap = mapVector
+    atIndex = (@>)
+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    minElement  = ap (@>) minIndex
+    maxElement  = ap (@>) maxIndex
+    sumElements  = sumC
+    prodElements = prodC
+
+instance Container 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]
+    --
+--instance Container Vector (Complex Float) where
+    cmap = mapVector
+    atIndex = (@>)
+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
+    minElement  = ap (@>) minIndex
+    maxElement  = ap (@>) maxIndex
+    sumElements  = sumQ
+    prodElements = prodQ
+
+---------------------------------------------------------------
+
+instance (Container Vector a) => Container 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]
+    --
+--instance (Container Vector a) => Container Matrix a where
+    cmap f = liftMatrix (mapVector f)
+    atIndex = (@@>)
+    minIndex m = let (r,c) = (rows m,cols m)
+                     i = (minIndex $ flatten m)
+                 in (i `div` c,(i `mod` c) + 1)
+    maxIndex m = let (r,c) = (rows m,cols m)
+                     i = (maxIndex $ flatten m)
+                 in (i `div` c,(i `mod` c) + 1)
+    minElement = ap (@@>) minIndex
+    maxElement = ap (@@>) maxIndex
+    sumElements = sumElements . flatten
+    prodElements = prodElements . flatten
+
+----------------------------------------------------
+
+
+-- | Linear algebraic properties of objects
+class Num e => Vectors a e where
+    -- | dot (inner) product
+    dot        :: a e -> a e -> e
+    -- | sum of absolute value of elements (differs in complex case from @norm1@
+    absSum     :: a e -> e
+    -- | sum of absolute value of elements
+    norm1      :: a e -> e
+    -- | euclidean norm
+    norm2      :: a e -> e
+    -- | element of maximum magnitude
+    normInf    :: a e -> e
+
+instance Vectors Vector Float where
+    norm2      = toScalarF Norm2
+    absSum     = toScalarF AbsSum
+    dot        = dotF
+    norm1      = toScalarF AbsSum
+    normInf    = maxElement . vectorMapF Abs
+
+instance Vectors Vector Double where
+    norm2      = toScalarR Norm2
+    absSum     = toScalarR AbsSum
+    dot        = dotR
+    norm1      = toScalarR AbsSum
+    normInf    = maxElement . vectorMapR Abs
+
+instance Vectors Vector (Complex Float) where
+    norm2      = (:+ 0) . toScalarQ Norm2
+    absSum     = (:+ 0) . toScalarQ AbsSum
+    dot        = dotQ
+    norm1      = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs
+    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs
+
+instance Vectors Vector (Complex Double) where
+    norm2      = (:+ 0) . toScalarC Norm2
+    absSum     = (:+ 0) . toScalarC AbsSum
+    dot        = dotC
+    norm1      = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs
+    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs
+
+----------------------------------------------------
+
+class Element t => Product t where
+    multiply :: Matrix t -> Matrix t -> Matrix t
+    ctrans :: Matrix t -> Matrix t
+
+instance Product Double where
+    multiply = multiplyR
+    ctrans = trans
+
+instance Product (Complex Double) where
+    multiply = multiplyC
+    ctrans = conj . trans
+
+instance Product Float where
+    multiply = multiplyF
+    ctrans = trans
+
+instance Product (Complex Float) where
+    multiply = multiplyQ
+    ctrans = conj . trans
+
+----------------------------------------------------------
+
+-- synonym for matrix product
+mXm :: Product t => Matrix t -> Matrix t -> Matrix t
+mXm = multiply
+
+-- matrix - vector product
+mXv :: Product t => Matrix t -> Vector t -> Vector t
+mXv m v = flatten $ m `mXm` (asColumn v)
+
+-- vector - matrix product
+vXm :: Product t => Vector t -> Matrix t -> Vector t
+vXm v m = flatten $ (asRow v) `mXm` m
+
+{- | Outer product of two vectors.
+
+@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
+(3><3)
+ [  5.0, 2.0, 3.0
+ , 10.0, 4.0, 6.0
+ , 15.0, 6.0, 9.0 ]@
+-}
+outer :: (Product t) => Vector t -> Vector t -> Matrix t
+outer u v = asColumn u `multiply` asRow v
+
+{- | Kronecker product of two matrices.
+
+@m1=(2><3)
+ [ 1.0,  2.0, 0.0
+ , 0.0, -1.0, 3.0 ]
+m2=(4><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0 ]@
+
+@\> kronecker m1 m2
+(8><9)
+ [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
+ ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
+ ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
+ , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
+ ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
+ ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
+ ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
+-}
+kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
+kronecker a b = fromBlocks
+              . splitEvery (cols a)
+              . map (reshape (cols b))
+              . toRows
+              $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index da4a4b5..edfd87e 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -14,14 +14,12 @@ This module reexports all normally required functions for Linear Algebra applica
 -----------------------------------------------------------------------------
 module Numeric.LinearAlgebra (
     module Numeric.LinearAlgebra.Algorithms,
-    module Numeric.LinearAlgebra.Interface,
-    module Numeric.LinearAlgebra.Linear,
+--    module Numeric.LinearAlgebra.Interface,
     module Numeric.Matrix,
-    module Numeric.Vector
+    module Data.Packed.Random
 ) where
 
 import Numeric.LinearAlgebra.Algorithms
-import Numeric.LinearAlgebra.Interface
-import Numeric.LinearAlgebra.Linear
+--import Numeric.LinearAlgebra.Interface
 import Numeric.Matrix
-import Numeric.Vector
+import Data.Packed.Random
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 8306961..fa6c475 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -77,15 +77,16 @@ module Numeric.LinearAlgebra.Algorithms (
 import Data.Packed.Internal hiding ((//))
 import Data.Packed.Matrix
 import Data.Complex
-import Numeric.LinearAlgebra.Linear
+--import Numeric.LinearAlgebra.Linear
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Data.List(foldl1')
 import Data.Array
-import Numeric.Vector
-import Numeric.Matrix()
+import Numeric.Container
+
+constant x = constantD x
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
-class (Product t, Linear Vector t, Container Vector t, Container Matrix t) => Field t where
+class (Product t, Container Vector t, Container Matrix t) => Field t where
     svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     sv'          :: Matrix t -> Vector Double
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
deleted file mode 100644
index 04a9d88..0000000
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ /dev/null
@@ -1,276 +0,0 @@
-{-# LANGUAGE UndecidableInstances, FlexibleInstances #-}
------------------------------------------------------------------------------
-{- |
-Module      :  Numeric.LinearAlgebra.Instances
-Copyright   :  (c) Alberto Ruiz 2006
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
-Stability   :  provisional
-Portability :  portable
-
-This module exports Show, Read, Eq, Num, Fractional, and Floating 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 Numeric.LinearAlgebra.Instances(
-) where
-
-import Numeric.LinearAlgebra.Linear
-import Numeric.GSL.Vector
-import Data.Packed.Matrix
-import Data.Complex
-import Data.List(transpose,intersperse)
-import Data.Packed.Internal.Vector
-
-#ifndef VECTOR
-import Foreign(Storable)
-#endif
-
-------------------------------------------------------------------
-
-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 ", "
-
-#ifndef VECTOR
-
-instance (Show a, Storable a) => (Show (Vector a)) where
-    show v = (show (dim v))++" |> " ++ show (toList v)
-
-#endif
-
-------------------------------------------------------------------
-
-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
-
-#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
-
-------------------------------------------------------------------
-
-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 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, 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 (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
-    fromRational n = (1><1) [fromRational n]
-    (/) = liftMatrix2Auto (/)
-
----------------------------------------------------------
-
-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 (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 (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
-
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
deleted file mode 100644
index fa3e209..0000000
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ /dev/null
@@ -1,118 +0,0 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
-{-# LANGUAGE UndecidableInstances #-}
------------------------------------------------------------------------------
-{- |
-Module      :  Numeric.LinearAlgebra.Interface
-Copyright   :  (c) Alberto Ruiz 2007
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
-Stability   :  provisional
-Portability :  portable
-
-Some useful operators, and Show, Read, Eq, Num, Fractional, and Floating 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 Numeric.LinearAlgebra.Interface(
-    (<>),(<.>),
-    (<\>),
-    (.*),(*/),
-    (<|>),(<->),
-) where
-
-import Numeric.Vector
-import Numeric.Matrix
-import Numeric.LinearAlgebra.Algorithms
-import Numeric.LinearAlgebra.Linear
-
-class Mul a b c | a b -> c where
- infixl 7 <>
- -- | Matrix-matrix, matrix-vector, and vector-matrix products.
- (<>)  :: Product t => a t -> b t -> c t
-
-instance Mul Matrix Matrix Matrix where
-    (<>) = mXm
-
-instance Mul Matrix Vector Vector where
-    (<>) m v = flatten $ m <> (asColumn v)
-
-instance Mul Vector Matrix Vector where
-    (<>) v m = flatten $ (asRow v) <> m
-
----------------------------------------------------
-
--- | Dot product: @u \<.\> v = dot u v@
---(<.>) :: (Field t) => Vector t -> Vector t -> t
-(<.>) :: Vectors Vector t => Vector t -> Vector t -> t
-infixl 7 <.>
-(<.>) = dot
-
-----------------------------------------------------
-
-{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
-
--- -- | @x .* a = scale x a@
--- (.*) :: (Linear c a) => a -> c a -> c a
-infixl 7 .*
-a .* x = scale a x
-
-----------------------------------------------------
-
-{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
-
--- -- | @a *\/ x = scale (recip x) a@
--- (*/) :: (Linear c a) => c a -> a -> c a
-infixl 7 */
-v */ x = scale (recip x) v
-
--- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
-(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
-infixl 7 <\>
-m <\> v = flatten (linearSolveSVD m (reshape 1 v))
-
-------------------------------------------------
-
-{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
-{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
-
-class Joinable a b where
-    joinH :: Element t => a t -> b t -> Matrix t
-    joinV :: Element t => a t -> b t -> Matrix t
-
-instance Joinable Matrix Matrix where
-    joinH m1 m2 = fromBlocks [[m1,m2]]
-    joinV m1 m2 = fromBlocks [[m1],[m2]]
-
-instance Joinable Matrix Vector where
-    joinH m v = joinH m (asColumn v)
-    joinV m v = joinV m (asRow v)
-
-instance Joinable Vector Matrix where
-    joinH v m = joinH (asColumn v) m
-    joinV v m = joinV (asRow v) m
-
-infixl 4 <|>
-infixl 3 <->
-
-{-- - | Horizontal concatenation of matrices and vectors:
-
-@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
-(6><4)
- [ 1.0, 0.0, 0.0, 1.0
- , 0.0, 1.0, 0.0, 2.0
- , 0.0, 0.0, 1.0, 3.0
- , 3.0, 0.0, 0.0, 4.0
- , 0.0, 3.0, 0.0, 5.0
- , 0.0, 0.0, 3.0, 6.0 ]@
--}
--- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
-a <|> b = joinH a b
-
--- -- | Vertical concatenation of matrices and vectors.
--- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
-a <-> b = joinV a b
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 5d0154d..288439f 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -44,8 +44,7 @@ module Numeric.LinearAlgebra.LAPACK (
 import Data.Packed.Internal
 import Data.Packed.Matrix
 import Data.Complex
-import Numeric.Vector()
-import Numeric.Container
+import Numeric.Conversion
 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
deleted file mode 100644
index 775060e..0000000
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ /dev/null
@@ -1,160 +0,0 @@
-{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE TypeFamilies #-}
------------------------------------------------------------------------------
-{- |
-Module      :  Numeric.LinearAlgebra.Linear
-Copyright   :  (c) Alberto Ruiz 2006-7
-License     :  GPL-style
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
-Stability   :  provisional
-Portability :  uses ffi
-
-Basic optimized operations on vectors and matrices.
-
--}
------------------------------------------------------------------------------
-
-module Numeric.LinearAlgebra.Linear (
-    -- * Linear Algebra Typeclasses
-    Vectors(..),
-    -- * Products
-    Product(..),
-    mXm,mXv,vXm,
-    outer, kronecker,
-    -- * Modules
-    --module Numeric.Vector,
-    --module Numeric.Matrix,
-    module Numeric.Container
-) where
-
-import Data.Packed.Internal.Common
-import Data.Packed.Matrix
-import Data.Packed.Vector
-import Data.Complex
-import Numeric.Container
-import Numeric.Vector()
-import Numeric.Matrix()
-import Numeric.GSL.Vector
-import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
-
--- | Linear algebraic properties of objects
-class Num e => Vectors a e where
-    -- | dot (inner) product
-    dot        :: a e -> a e -> e
-    -- | sum of absolute value of elements (differs in complex case from @norm1@
-    absSum     :: a e -> e
-    -- | sum of absolute value of elements
-    norm1      :: a e -> e
-    -- | euclidean norm
-    norm2      :: a e -> e
-    -- | element of maximum magnitude
-    normInf    :: a e -> e
-
-instance Vectors Vector Float where
-    norm2      = toScalarF Norm2
-    absSum     = toScalarF AbsSum
-    dot        = dotF
-    norm1      = toScalarF AbsSum
-    normInf    = maxElement . vectorMapF Abs
-
-instance Vectors Vector Double where
-    norm2      = toScalarR Norm2
-    absSum     = toScalarR AbsSum
-    dot        = dotR
-    norm1      = toScalarR AbsSum
-    normInf    = maxElement . vectorMapR Abs
-
-instance Vectors Vector (Complex Float) where
-    norm2      = (:+ 0) . toScalarQ Norm2
-    absSum     = (:+ 0) . toScalarQ AbsSum
-    dot        = dotQ
-    norm1      = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs
-    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs
-
-instance Vectors Vector (Complex Double) where
-    norm2      = (:+ 0) . toScalarC Norm2
-    absSum     = (:+ 0) . toScalarC AbsSum
-    dot        = dotC
-    norm1      = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs
-    normInf    = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs
-
-----------------------------------------------------
-
-class Element t => Product t where
-    multiply :: Matrix t -> Matrix t -> Matrix t
-    ctrans :: Matrix t -> Matrix t
-
-instance Product Double where
-    multiply = multiplyR
-    ctrans = trans
-
-instance Product (Complex Double) where
-    multiply = multiplyC
-    ctrans = conj . trans
-
-instance Product Float where
-    multiply = multiplyF
-    ctrans = trans
-
-instance Product (Complex Float) where
-    multiply = multiplyQ
-    ctrans = conj . trans
-
-----------------------------------------------------------
-
--- synonym for matrix product
-mXm :: Product t => Matrix t -> Matrix t -> Matrix t
-mXm = multiply
-
--- matrix - vector product
-mXv :: Product t => Matrix t -> Vector t -> Vector t
-mXv m v = flatten $ m `mXm` (asColumn v)
-
--- vector - matrix product
-vXm :: Product t => Vector t -> Matrix t -> Vector t
-vXm v m = flatten $ (asRow v) `mXm` m
-
-{- | Outer product of two vectors.
-
-@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
-(3><3)
- [  5.0, 2.0, 3.0
- , 10.0, 4.0, 6.0
- , 15.0, 6.0, 9.0 ]@
--}
-outer :: (Product t) => Vector t -> Vector t -> Matrix t
-outer u v = asColumn u `multiply` asRow v
-
-{- | Kronecker product of two matrices.
-
-@m1=(2><3)
- [ 1.0,  2.0, 0.0
- , 0.0, -1.0, 3.0 ]
-m2=(4><3)
- [  1.0,  2.0,  3.0
- ,  4.0,  5.0,  6.0
- ,  7.0,  8.0,  9.0
- , 10.0, 11.0, 12.0 ]@
-
-@\> kronecker m1 m2
-(8><9)
- [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
- ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
- ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
- , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
- ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
- ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
- ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
- ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
--}
-kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
-kronecker a b = fromBlocks
-              . splitEvery (cols a)
-              . map (reshape (cols b))
-              . toRows
-              $ flatten a `outer` flatten b
-
-
--------------------------------------------------------------------
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
index fa3f94a..6ba870f 100644
--- a/lib/Numeric/Matrix.hs
+++ b/lib/Numeric/Matrix.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -14,12 +15,22 @@
 -- Stability   :  provisional
 -- Portability :  portable
 --
--- Numeric instances and functions for 'Data.Packed.Matrix's
+-- Numeric instances and functions for '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 Numeric.Matrix (
-                       module Data.Packed.Matrix,
+    -- * Basic functions
+    module Data.Packed.Matrix,
+    module Numeric.Vector,
+    --module Numeric.Container,
+    -- * Operators
+    (<>), (<\>),
+    -- * Deprecated
+    (.*),(*/),(<|>),(<->),
+    vectorMax,vectorMin
                       ) where
 
 -------------------------------------------------------------------
@@ -28,18 +39,18 @@ import Data.Packed.Vector
 import Data.Packed.Matrix
 import Numeric.Container
 --import Numeric.LinearAlgebra.Linear
-import Numeric.Vector()
+import Numeric.Vector
+import Numeric.LinearAlgebra.Algorithms
+--import Control.Monad(ap)
 
-import Control.Monad(ap)
-
-import Control.Arrow((***))
+--import Control.Arrow((***))
 
 -------------------------------------------------------------------
 
-instance Linear Matrix a => Eq (Matrix a) where
+instance Container Matrix a => Eq (Matrix a) where
     (==) = equal
 
-instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
+instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where
     (+) = liftMatrix2Auto (+)
     (-) = liftMatrix2Auto (-)
     negate = liftMatrix negate
@@ -50,13 +61,13 @@ instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
 
 ---------------------------------------------------
 
-instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
+instance (Container 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
+instance (Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
     sin   = liftMatrix sin
     cos   = liftMatrix cos
     tan   = liftMatrix tan
@@ -75,43 +86,93 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
     sqrt  = liftMatrix sqrt
     pi    = (1><1) [pi]
 
----------------------------------------------------------------
-
-instance NumericContainer 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'
-
----------------------------------------------------------------
-
-instance (Linear Vector a, Container Matrix a) => 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]
-    --
-instance (Container Vector a) => Container Matrix a where
-    cmap f = liftMatrix (mapVector f)
-    atIndex = (@@>)
-    minIndex m = let (r,c) = (rows m,cols m)
-                     i = (minIndex $ flatten m)
-                 in (i `div` c,(i `mod` c) + 1)
-    maxIndex m = let (r,c) = (rows m,cols m)
-                     i = (maxIndex $ flatten m)
-                 in (i `div` c,(i `mod` c) + 1)
-    minElement = ap (@@>) minIndex
-    maxElement = ap (@@>) maxIndex
-    sumElements = sumElements . flatten
-    prodElements = prodElements . flatten
+--------------------------------------------------------
+
+class Mul a b c | a b -> c where
+ infixl 7 <>
+ -- | Matrix-matrix, matrix-vector, and vector-matrix products.
+ (<>)  :: Product t => a t -> b t -> c t
+
+instance Mul Matrix Matrix Matrix where
+    (<>) = mXm
+
+instance Mul Matrix Vector Vector where
+    (<>) m v = flatten $ m <> (asColumn v)
+
+instance Mul Vector Matrix Vector where
+    (<>) v m = flatten $ (asRow v) <> m
 
 ----------------------------------------------------
+
+{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
+
+-- -- | @x .* a = scale x a@
+-- (.*) :: (Linear c a) => a -> c a -> c a
+infixl 7 .*
+a .* x = scale a x
+
+----------------------------------------------------
+
+{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
+
+-- -- | @a *\/ x = scale (recip x) a@
+-- (*/) :: (Linear c a) => c a -> a -> c a
+infixl 7 */
+v */ x = scale (recip x) v
+
+-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
+(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
+infixl 7 <\>
+m <\> v = flatten (linearSolveSVD m (reshape 1 v))
+
+------------------------------------------------
+
+{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
+{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
+
+class Joinable a b where
+    joinH :: Element t => a t -> b t -> Matrix t
+    joinV :: Element t => a t -> b t -> Matrix t
+
+instance Joinable Matrix Matrix where
+    joinH m1 m2 = fromBlocks [[m1,m2]]
+    joinV m1 m2 = fromBlocks [[m1],[m2]]
+
+instance Joinable Matrix Vector where
+    joinH m v = joinH m (asColumn v)
+    joinV m v = joinV m (asRow v)
+
+instance Joinable Vector Matrix where
+    joinH v m = joinH (asColumn v) m
+    joinV v m = joinV (asRow v) m
+
+infixl 4 <|>
+infixl 3 <->
+
+{-- - | Horizontal concatenation of matrices and vectors:
+
+@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
+(6><4)
+ [ 1.0, 0.0, 0.0, 1.0
+ , 0.0, 1.0, 0.0, 2.0
+ , 0.0, 0.0, 1.0, 3.0
+ , 3.0, 0.0, 0.0, 4.0
+ , 0.0, 3.0, 0.0, 5.0
+ , 0.0, 0.0, 3.0, 6.0 ]@
+-}
+-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
+a <|> b = joinH a b
+
+-- -- | Vertical concatenation of matrices and vectors.
+-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
+a <-> b = joinV a b
+
+-------------------------------------------------------------------
+
+{-# DEPRECATED vectorMin "use minElement" #-}
+vectorMin :: (Container Vector t, Element t) => Vector t -> t
+vectorMin = minElement
+
+{-# DEPRECATED vectorMax "use maxElement" #-}
+vectorMax :: (Container Vector t, Element t) => Vector t -> t
+vectorMax = maxElement
diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs
index 55d645a..9427243 100644
--- a/lib/Numeric/Vector.hs
+++ b/lib/Numeric/Vector.hs
@@ -14,24 +14,26 @@
 -- Stability   :  provisional
 -- Portability :  portable
 --
--- Numeric instances and functions for 'Data.Packed.Vector's
+-- Numeric instances and functions for 'Vector'.
 --
 -----------------------------------------------------------------------------
 
 module Numeric.Vector (
-                       -- * Vector creation
-                       constant, linspace,
-                       module Data.Packed.Vector
+    -- * Basic vector functions
+    module Data.Packed.Vector,
+    module Numeric.Container,
+    -- * Vector creation
+    constant, linspace,
+    -- * Operators
+    (<.>)
                       ) where
 
 import Data.Complex
 
-import Control.Monad(ap)
-
 import Data.Packed.Vector
-import Data.Packed.Internal.Matrix(Element(..))
-import Data.Packed.Internal.Vector(asComplex,asReal)
-import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape)
+import Data.Packed.Internal.Matrix
+--import Data.Packed.Internal.Vector(asComplex,asReal)
+--import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape)
 import Numeric.GSL.Vector
 
 import Numeric.Container
@@ -92,10 +94,15 @@ 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 :: (Enum e, Container 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)
 
+-- | Dot product: @u \<.\> v = dot u v@
+(<.>) :: Vectors Vector t => Vector t -> Vector t -> t
+infixl 7 <.>
+(<.>) = dot
+
 ------------------------------------------------------------------
 
 adaptScalar f1 f2 f3 x y
@@ -107,7 +114,7 @@ adaptScalar f1 f2 f3 x y
 
 #ifndef VECTOR
 
-instance Linear Vector a => Eq (Vector a) where
+instance Container Vector a => Eq (Vector a) where
     (==) = equal
 
 #endif
@@ -146,7 +153,7 @@ instance Num (Vector (Complex Float)) where
 
 ---------------------------------------------------
 
-instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where
+instance (Container 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
@@ -254,116 +261,6 @@ instance Floating (Vector (Complex Float)) where
 -- instance (NFData a, Element a) => NFData (Matrix a) where
 --     rnf = rnf . flatten
 
----------------------------------------------------------------
-
--- | obtains the complex conjugate of a complex vector
-conjV :: (RealElement a) => Vector (Complex a) -> Vector (Complex a)
-conjV = mapVector conjugate
-
--- | creates a complex vector from vectors with real and imaginary parts
-toComplexV :: (RealElement a) => (Vector a, Vector a) ->  Vector (Complex a)
-toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
-
--- | the inverse of 'toComplex'
-fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)
-fromComplexV z = (r,i) where
-    [r,i] = toColumns $ reshape 2 $ asReal z
 
 --------------------------------------------------------------------------
 
-instance NumericContainer Vector where
-    toComplex = toComplexV
-    fromComplex = fromComplexV
-    complex' v = toComplex (v,constant 0 (dim v))
-    conj = conjV
---    cmap = mapVector
-    single' = double2FloatG
-    double' = float2DoubleG
-
---------------------------------------------------------------------------
-
-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]
-    --
-instance Container Vector Float where
-    cmap = mapVector
-    atIndex = (@>)
-    minIndex     = round . toScalarF MinIdx
-    maxIndex     = round . toScalarF MaxIdx
-    minElement  = toScalarF Min
-    maxElement  = toScalarF Max
-    sumElements  = sumF
-    prodElements = prodF
-
-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]
-    --
-instance Container Vector Double where
-    cmap = mapVector
-    atIndex = (@>)
-    minIndex     = round . toScalarR MinIdx
-    maxIndex     = round . toScalarR MaxIdx
-    minElement  = toScalarR Min
-    maxElement  = toScalarR Max
-    sumElements  = sumR
-    prodElements = prodR
-
-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]
-    --
-instance Container Vector (Complex Double) where
-    cmap = mapVector
-    atIndex = (@>)
-    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    minElement  = ap (@>) minIndex
-    maxElement  = ap (@>) maxIndex
-    sumElements  = sumC
-    prodElements = prodC
-
-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]
-    --
-instance Container Vector (Complex Float) where
-    cmap = mapVector
-    atIndex = (@>)
-    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    minElement  = ap (@>) minIndex
-    maxElement  = ap (@>) maxIndex
-    sumElements  = sumQ
-    prodElements = prodQ
-
----------------------------------------------------------------