43 files changed, 594 insertions, 717 deletions

diff --git a/packages/Makefile b/packages/Makefile
index e9d8586..b00d71f 100644
--- a/packages/Makefile
+++ b/packages/Makefile
@@ -1,22 +1,26 @@
-pkgs=base gsl special glpk tests ../../hTensor ../../easyVision/packages/base
+pkgs=base gsl special glpk tests ../../hTensor ../../easyVision/packages/tools ../../easyVision/packages/base
+
+mkl=--extra-include-dirs=$(MKL) --extra-lib-dirs=$(MKL)
+
+cabalcmd =                         \
+    for p in $(1); do              \
+          if [ -e $$p ]; then          \
+            cd $$p; cabal $(2) ; cd -; \
+          fi;                          \
+        done;                          \
+        cd sparse;                     \
+        cabal $(3) $(2); cd -;
+
 
 all:
-        for p in $(pkgs); do \
-          if [ -e $$p ]; then \
-            cd $$p; cabal install --force-reinstall --enable-documentation ; cd -; \
-          fi; \
-        done
-        cd sparse; \
-        cabal install --extra-include-dirs=$(MKL) --extra-lib-dirs=$(MKL) \
-                  --force-reinstall --enable-documentation ; cd -;
+        $(call cabalcmd, $(pkgs), install --force-reinstall --enable-documentation, $(mkl))
 
 fast:
-        for p in $(pkgs); do \
-          if [ -e $$p ]; then \
-            cd $$p; cabal install --force-reinstall ; cd -; \
-          fi; \
-        done
-        cd sparse; \
-        cabal install --extra-include-dirs=$(MKL) --extra-lib-dirs=$(MKL) \
-                  --force-reinstall; cd -;
+        $(call cabalcmd, $(pkgs), install --force-reinstall, $(mkl))
+
+clean:
+        $(call cabalcmd, $(pkgs), clean)
+
+prof:
+        $(call cabalcmd, $(pkgs), install --force-reinstall --enable-library-profiling, $(mkl))
 
diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index c137285..ce19c97 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -1,3 +1,10 @@
+0.17.0.0
+--------
+
+    * old compatibility modules removed
+
+    * added "unitary" and "pairwiseD2"
+
 0.16.1.0
 --------
 
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 3895dc1..7bfc762 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.16.1.5
+Version:            0.17.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -9,16 +9,6 @@ Homepage:           https://github.com/albertoruiz/hmatrix
 Synopsis:           Numeric Linear Algebra
 Description:        Linear algebra based on BLAS and LAPACK.
                     .
-                    The package is organized as follows:
-                    .
-                    ["Numeric.LinearAlgebra.HMatrix"] Starting point and recommended import module for most applications.
-                    .
-                    ["Numeric.LinearAlgebra.Static"] Experimental alternative interface.
-                    .
-                    ["Numeric.LinearAlgebra.Devel"] Tools for extending the library.
-                    .
-                    (Other modules are exposed with hidden documentation for backwards compatibility.)
-                    .
                     Code examples: <http://dis.um.es/~alberto/hmatrix/hmatrix.html>
 
 Category:           Math
@@ -39,7 +29,7 @@ flag openblas
 
 library
 
-    Build-Depends:      base >= 4 && < 5,
+    Build-Depends:      base >= 4.8 && < 5,
                         binary,
                         array,
                         deepseq,
@@ -51,27 +41,19 @@ library
 
     hs-source-dirs:     src
 
-    exposed-modules:    Data.Packed,
-                        Data.Packed.Vector,
-                        Data.Packed.Matrix,
-                        Data.Packed.Foreign,
-                        Data.Packed.ST,
-                        Data.Packed.Development,
-
-                        Numeric.LinearAlgebra
-                        Numeric.LinearAlgebra.LAPACK
-                        Numeric.LinearAlgebra.Algorithms
-                        Numeric.Container
-                        Numeric.LinearAlgebra.Util
-
+    exposed-modules:    Numeric.LinearAlgebra
                         Numeric.LinearAlgebra.Devel
                         Numeric.LinearAlgebra.Data
                         Numeric.LinearAlgebra.HMatrix
                         Numeric.LinearAlgebra.Static
-                        
-
 
-    other-modules:      Data.Packed.Internal,
+    other-modules:      Data.Packed,
+                        Data.Packed.Vector
+                        Data.Packed.Matrix
+                        Data.Packed.Foreign
+                        Data.Packed.ST
+                        Data.Packed.Development
+                        Data.Packed.Internal
                         Data.Packed.Internal.Common
                         Data.Packed.Internal.Signatures
                         Data.Packed.Internal.Vector
@@ -83,12 +65,16 @@ library
                         Numeric.Matrix
                         Data.Packed.Internal.Numeric
                         Data.Packed.Numeric
+                        Numeric.LinearAlgebra.LAPACK
+                        Numeric.LinearAlgebra.Algorithms
+                        Numeric.Container
                         Numeric.LinearAlgebra.Util.Convolution
                         Numeric.LinearAlgebra.Util.CG
                         Numeric.LinearAlgebra.Random
                         Numeric.Conversion
                         Numeric.Sparse
                         Numeric.LinearAlgebra.Static.Internal
+                        Numeric.LinearAlgebra.Util
 
     C-sources:          src/C/lapack-aux.c
                         src/C/vector-aux.c
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
index 85f1b37..b0999a8 100644
--- a/packages/base/src/Data/Packed/IO.hs
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -22,7 +22,6 @@ import Text.Printf(printf)
 import Data.List(intersperse)
 import Data.Complex
 import Numeric.Vectorized(vectorScan,saveMatrix)
-import Control.Applicative((<$>))
 import Data.Packed.Internal
 
 {- | Creates a string from a matrix given a separator and a function to show each entry. Using
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 257ad73..7a4dd29 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -48,7 +48,6 @@ import Numeric.Conversion
 import Data.Packed.Development
 import Numeric.Vectorized
 import Data.Complex
-import Control.Applicative((<*>))
 
 import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
 import Data.Packed.Internal
diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs
index d0bc143..b49f379 100644
--- a/packages/base/src/Data/Packed/Internal/Vector.hs
+++ b/packages/base/src/Data/Packed/Internal/Vector.hs
@@ -35,7 +35,6 @@ import Foreign.Ptr(Ptr)
 import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
 import Foreign.C.Types
 import Data.Complex
-import Control.Monad(when)
 import System.IO.Unsafe(unsafePerformIO)
 
 #if __GLASGOW_HASKELL__ >= 605
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
index 6027f43..6d62f22 100644
--- a/packages/base/src/Data/Packed/Numeric.hs
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -39,7 +39,7 @@ module Data.Packed.Numeric (
     step, cond, find, assoc, accum,
     Transposable(..), Linear(..),
     -- * Matrix product
-    Product(..), udot, dot, (<·>), (#>), app,
+    Product(..), udot, dot, (<·>), (#>), (<#), app,
     Mul(..),
     (<.>),
     optimiseMult,
@@ -71,7 +71,6 @@ import Data.Packed
 import Data.Packed.Internal.Numeric
 import Data.Complex
 import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
-import Data.Monoid(Monoid(mconcat))
 import Data.Packed.IO
 import Numeric.LinearAlgebra.Random
 
@@ -142,6 +141,11 @@ fromList [140.0,320.0]
 app :: Numeric t => Matrix t -> Vector t -> Vector t
 app = (#>)
 
+infixl 8 <#
+-- | dense vector-matrix product
+(<#) :: Numeric t => Vector t -> Matrix t -> Vector t
+(<#) = vXm
+
 --------------------------------------------------------------------------------
 
 class Mul a b c | a b -> c where
diff --git a/packages/base/src/Numeric/Chain.hs b/packages/base/src/Numeric/Chain.hs
index 443bd28..64c09c0 100644
--- a/packages/base/src/Numeric/Chain.hs
+++ b/packages/base/src/Numeric/Chain.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE FlexibleContexts #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Chain
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index ad315e4..4ba0c98 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -1,22 +1,235 @@
---------------------------------------------------------------------------------
+-----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra
-Copyright   :  (c) Alberto Ruiz 2006-14
+Copyright   :  (c) Alberto Ruiz 2006-15
 License     :  BSD3
 Maintainer  :  Alberto Ruiz
 Stability   :  provisional
 
 -}
---------------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
-
+-----------------------------------------------------------------------------
 module Numeric.LinearAlgebra (
-    module Numeric.Container,
-    module Numeric.LinearAlgebra.Algorithms
+
+    -- * Basic types and data processing
+    module Numeric.LinearAlgebra.Data,
+
+    -- * Arithmetic and numeric classes
+    -- |
+    -- The standard numeric classes are defined elementwise:
+    --
+    -- >>>  vector [1,2,3] * vector [3,0,-2]
+    -- fromList [3.0,0.0,-6.0]
+    --
+    -- >>> matrix 3 [1..9] * ident 3
+    -- (3><3)
+    --  [ 1.0, 0.0, 0.0
+    --  , 0.0, 5.0, 0.0
+    --  , 0.0, 0.0, 9.0 ]
+    --
+    -- In arithmetic operations single-element vectors and matrices
+    -- (created from numeric literals or using 'scalar') automatically
+    -- expand to match the dimensions of the other operand:
+    --
+    -- >>> 5 + 2*ident 3 :: Matrix Double
+    -- (3><3)
+    --  [ 7.0, 5.0, 5.0
+    --  , 5.0, 7.0, 5.0
+    --  , 5.0, 5.0, 7.0 ]
+    --
+    -- >>> matrix 3 [1..9] + matrix 1 [10,20,30]
+    -- (3><3)
+    --  [ 11.0, 12.0, 13.0
+    --  , 24.0, 25.0, 26.0
+    --  , 37.0, 38.0, 39.0 ]
+    --
+
+    -- * Products
+    -- ** dot
+    dot, (<·>),
+    -- ** matrix-vector
+    app, (#>), (<#), (!#>),
+    -- ** matrix-matrix
+    mul, (<>),
+    -- | The matrix product is also implemented in the "Data.Monoid" instance, where
+    -- single-element matrices (created from numeric literals or using 'scalar')
+    -- are used for scaling.
+    --
+    -- >>> import Data.Monoid as M
+    -- >>>  let m = matrix 3 [1..6]
+    -- >>> m M.<> 2 M.<> diagl[0.5,1,0]
+    -- (2><3)
+    --  [ 1.0,  4.0, 0.0
+    --  , 4.0, 10.0, 0.0 ]
+    --
+    -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
+
+
+    -- ** other
+    outer, kronecker, cross,
+    scale,
+    sumElements, prodElements,
+
+    -- * Linear Systems
+    (<\>),
+    linearSolve,
+    linearSolveLS,
+    linearSolveSVD,
+    luSolve,
+    cholSolve,
+    cgSolve,
+    cgSolve',
+
+    -- * Inverse and pseudoinverse
+    inv, pinv, pinvTol,
+
+    -- * Determinant and rank
+    rcond, rank,
+    det, invlndet,
+
+    -- * Norms
+    Normed(..),
+    norm_Frob, norm_nuclear,
+
+    -- * Nullspace and range
+    orth,
+    nullspace, null1, null1sym,
+
+    -- * SVD
+    svd,
+    thinSVD,
+    compactSVD,
+    singularValues,
+    leftSV, rightSV,
+
+    -- * Eigensystems
+    eig, eigSH, eigSH',
+    eigenvalues, eigenvaluesSH, eigenvaluesSH',
+    geigSH',
+
+    -- * QR
+    qr, rq, qrRaw, qrgr,
+
+    -- * Cholesky
+    chol, cholSH, mbCholSH,
+
+    -- * Hessenberg
+    hess,
+
+    -- * Schur
+    schur,
+
+    -- * LU
+    lu, luPacked,
+
+    -- * Matrix functions
+    expm,
+    sqrtm,
+    matFunc,
+
+    -- * Correlation and convolution
+    corr, conv, corrMin, corr2, conv2,
+
+    -- * Random arrays
+
+    Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
+
+    -- * Misc
+    meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
+    ℝ,ℂ,iC,
+    -- * Auxiliary classes
+    Element, Container, Product, Numeric, LSDiv,
+    Complexable, RealElement,
+    RealOf, ComplexOf, SingleOf, DoubleOf,
+    IndexOf,
+    Field,
+--    Normed,
+    Transposable,
+    CGState(..),
+    Testable(..)
 ) where
 
-import Numeric.Container
-import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Data
+
 import Numeric.Matrix()
 import Numeric.Vector()
+import Data.Packed.Numeric hiding ((<>), mul)
+import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
+import qualified Numeric.LinearAlgebra.Algorithms as A
+import Numeric.LinearAlgebra.Util
+import Numeric.LinearAlgebra.Random
+import Numeric.Sparse((!#>))
+import Numeric.LinearAlgebra.Util.CG
+
+{- | infix synonym of 'mul'
+
+>>> let a = (3><5) [1..]
+>>> a
+(3><5)
+ [  1.0,  2.0,  3.0,  4.0,  5.0
+ ,  6.0,  7.0,  8.0,  9.0, 10.0
+ , 11.0, 12.0, 13.0, 14.0, 15.0 ]
+
+>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
+>>> b
+(5><2)
+ [  1.0, 3.0
+ ,  0.0, 2.0
+ , -1.0, 5.0
+ ,  7.0, 7.0
+ ,  6.0, 0.0 ]
+
+>>> a <> b
+(3><2)
+ [  56.0,  50.0
+ , 121.0, 135.0
+ , 186.0, 220.0 ]
+
+-}
+(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t
+(<>) = mXm
+infixr 8 <>
+
+-- | dense matrix product
+mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
+mul = mXm
+
+
+{- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
+
+@
+a = (2><2)
+ [ 1.0, 2.0
+ , 3.0, 5.0 ]
+@
+
+@
+b = (2><3)
+ [  6.0, 1.0, 10.0
+ , 15.0, 3.0, 26.0 ]
+@
+
+>>> linearSolve a b
+Just (2><3)
+ [ -1.4802973661668753e-15,     0.9999999999999997, 1.999999999999997
+ ,       3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ]
+
+>>> let Just x = it
+>>> disp 5 x
+2x3
+-0.00000  1.00000  2.00000
+ 3.00000  0.00000  4.00000
+
+>>> a <> x
+(2><3)
+ [  6.0, 1.0, 10.0
+ , 15.0, 3.0, 26.0 ]
+
+-}
+linearSolve m b = A.mbLinearSolve m b
+
+-- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
+nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
+
+-- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
+orth m = orthSVD (Left (1*eps)) m (leftSV m)
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index 02ac6a0..0e085c1 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -228,7 +228,9 @@ fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
 
 -}
 svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
-svd = {-# SCC "svd" #-} svd'
+svd = {-# SCC "svd" #-} g . svd'
+  where
+    g (u,s,v) = (u,s,tr v)
 
 {- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
 
@@ -272,7 +274,10 @@ fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
 
 -}
 thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
-thinSVD = {-# SCC "thinSVD" #-} thinSVD'
+thinSVD = {-# SCC "thinSVD" #-} g . thinSVD'
+  where
+    g (u,s,v) = (u,s,tr v)
+
 
 -- | Singular values only.
 singularValues :: Field t => Matrix t -> Vector Double
@@ -587,7 +592,7 @@ m = (3><3) [ 1, 0,    0
 -}
 
 pinvTol :: Field t => Double -> Matrix t -> Matrix t
-pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where
+pinvTol t m = v' `mXm` diag s' `mXm` ctrans u' where
     (u,s,v) = thinSVD m
     sl@(g:_) = toList s
     s' = real . fromList . map rec $ sl
@@ -649,7 +654,7 @@ nullspaceSVD hint a (s,v) = vs where
     k = case hint of
         Right t -> t
         _       -> rankSVD tol a s
-    vs = conj (dropColumns k v)
+    vs = dropColumns k v
 
 
 -- | The nullspace of a matrix. See also 'nullspaceSVD'.
@@ -935,10 +940,9 @@ relativeError' x y = dig (norm (x `sub` y) / norm x)
           dig r = round $ -logBase 10 (realToFrac r :: Double)
 
 
-relativeError :: (Normed c t, Num (c t)) => NormType -> c t -> c t -> Double
-relativeError t a b = realToFrac r
+relativeError :: Num a => (a -> Double) -> a -> a -> Double
+relativeError norm a b = r
   where
-    norm = pnorm t
     na = norm a
     nb = norm b
     nab = norm (a-b)
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index 677f9ee..8e67eb4 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -1,4 +1,4 @@
------------------------------------------------------------------------------
+--------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.HMatrix
 Copyright   :  (c) Alberto Ruiz 2006-14
@@ -7,229 +7,11 @@ Maintainer  :  Alberto Ruiz
 Stability   :  provisional
 
 -}
------------------------------------------------------------------------------
-module Numeric.LinearAlgebra.HMatrix (
-
-    -- * Basic types and data processing
-    module Numeric.LinearAlgebra.Data,
-
-    -- * Arithmetic and numeric classes
-    -- |
-    -- The standard numeric classes are defined elementwise:
-    --
-    -- >>>  vector [1,2,3] * vector [3,0,-2]
-    -- fromList [3.0,0.0,-6.0]
-    --
-    -- >>> matrix 3 [1..9] * ident 3
-    -- (3><3)
-    --  [ 1.0, 0.0, 0.0
-    --  , 0.0, 5.0, 0.0
-    --  , 0.0, 0.0, 9.0 ]
-    --
-    -- In arithmetic operations single-element vectors and matrices
-    -- (created from numeric literals or using 'scalar') automatically
-    -- expand to match the dimensions of the other operand:
-    --
-    -- >>> 5 + 2*ident 3 :: Matrix Double
-    -- (3><3)
-    --  [ 7.0, 5.0, 5.0
-    --  , 5.0, 7.0, 5.0
-    --  , 5.0, 5.0, 7.0 ]
-    --
-    -- >>> matrix 3 [1..9] + matrix 1 [10,20,30]
-    -- (3><3)
-    --  [ 11.0, 12.0, 13.0
-    --  , 24.0, 25.0, 26.0
-    --  , 37.0, 38.0, 39.0 ]
-    --
-
-    -- * Products
-    -- ** dot
-    dot, (<·>),
-    -- ** matrix-vector
-    app, (#>), (!#>),
-    -- ** matrix-matrix
-    mul, (<>),
-    -- | The matrix product is also implemented in the "Data.Monoid" instance, where
-    -- single-element matrices (created from numeric literals or using 'scalar')
-    -- are used for scaling.
-    --
-    -- >>> import Data.Monoid as M
-    -- >>>  let m = matrix 3 [1..6]
-    -- >>> m M.<> 2 M.<> diagl[0.5,1,0]
-    -- (2><3)
-    --  [ 1.0,  4.0, 0.0
-    --  , 4.0, 10.0, 0.0 ]
-    --
-    -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
-
-
-    -- ** other
-    outer, kronecker, cross,
-    scale,
-    sumElements, prodElements,
-
-    -- * Linear Systems
-    (<\>),
-    linearSolve,
-    linearSolveLS,
-    linearSolveSVD,
-    luSolve,
-    cholSolve,
-    cgSolve,
-    cgSolve',
-
-    -- * Inverse and pseudoinverse
-    inv, pinv, pinvTol,
-
-    -- * Determinant and rank
-    rcond, rank,
-    det, invlndet,
-
-    -- * Norms
-    Normed(..),
-    norm_Frob, norm_nuclear,
-
-    -- * Nullspace and range
-    orth,
-    nullspace, null1, null1sym,
-
-    -- * SVD
-    svd,
-    thinSVD,
-    compactSVD,
-    singularValues,
-    leftSV, rightSV,
-
-    -- * Eigensystems
-    eig, eigSH, eigSH',
-    eigenvalues, eigenvaluesSH, eigenvaluesSH',
-    geigSH',
-
-    -- * QR
-    qr, rq, qrRaw, qrgr,
-
-    -- * Cholesky
-    chol, cholSH, mbCholSH,
-
-    -- * Hessenberg
-    hess,
-
-    -- * Schur
-    schur,
-
-    -- * LU
-    lu, luPacked,
-
-    -- * Matrix functions
-    expm,
-    sqrtm,
-    matFunc,
-
-    -- * Correlation and convolution
-    corr, conv, corrMin, corr2, conv2,
-
-    -- * Random arrays
+--------------------------------------------------------------------------------
 
-    Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
-
-    -- * Misc
-    meanCov, rowOuters, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
-    ℝ,ℂ,iC,
-    -- * Auxiliary classes
-    Element, Container, Product, Numeric, LSDiv,
-    Complexable, RealElement,
-    RealOf, ComplexOf, SingleOf, DoubleOf,
-    IndexOf,
-    Field,
---    Normed,
-    Transposable,
-    CGState(..),
-    Testable(..)
+module Numeric.LinearAlgebra.HMatrix (
+    module Numeric.LinearAlgebra
 ) where
 
-import Numeric.LinearAlgebra.Data
-
-import Numeric.Matrix()
-import Numeric.Vector()
-import Data.Packed.Numeric hiding ((<>), mul)
-import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
-import qualified Numeric.LinearAlgebra.Algorithms as A
-import Numeric.LinearAlgebra.Util
-import Numeric.LinearAlgebra.Random
-import Numeric.Sparse((!#>))
-import Numeric.LinearAlgebra.Util.CG
-
-{- | infix synonym of 'mul'
-
->>> let a = (3><5) [1..]
->>> a
-(3><5)
- [  1.0,  2.0,  3.0,  4.0,  5.0
- ,  6.0,  7.0,  8.0,  9.0, 10.0
- , 11.0, 12.0, 13.0, 14.0, 15.0 ]
-
->>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
->>> b
-(5><2)
- [  1.0, 3.0
- ,  0.0, 2.0
- , -1.0, 5.0
- ,  7.0, 7.0
- ,  6.0, 0.0 ]
-
->>> a <> b
-(3><2)
- [  56.0,  50.0
- , 121.0, 135.0
- , 186.0, 220.0 ]
-
--}
-(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t
-(<>) = mXm
-infixr 8 <>
-
--- | dense matrix product
-mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
-mul = mXm
-
-
-{- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
-
-@
-a = (2><2)
- [ 1.0, 2.0
- , 3.0, 5.0 ]
-@
-
-@
-b = (2><3)
- [  6.0, 1.0, 10.0
- , 15.0, 3.0, 26.0 ]
-@
-
->>> linearSolve a b
-Just (2><3)
- [ -1.4802973661668753e-15,     0.9999999999999997, 1.999999999999997
- ,       3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ]
-
->>> let Just x = it
->>> disp 5 x
-2x3
--0.00000  1.00000  2.00000
- 3.00000  0.00000  4.00000
-
->>> a <> x
-(2><3)
- [  6.0, 1.0, 10.0
- , 15.0, 3.0, 26.0 ]
-
--}
-linearSolve m b = A.mbLinearSolve m b
-
--- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
-nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
-
--- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
-orth m = orthSVD (Left (1*eps)) m (leftSV m)
+import Numeric.LinearAlgebra
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
index e088fdc..d96e14f 100644
--- a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
@@ -115,7 +115,7 @@ svdAux f st x = unsafePerformIO $ do
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
     app4 f mat x mat u vec s mat v st
-    return (u,s,trans v)
+    return (u,s,v)
   where r = rows x
         c = cols x
 
@@ -141,7 +141,7 @@ thinSVDAux f st x = unsafePerformIO $ do
     s <- createVector q
     v <- createMatrix ColumnMajor q c
     app4 f mat x mat u vec s mat v st
-    return (u,s,trans v)
+    return (u,s,v)
   where r = rows x
         c = cols x
         q = min r c
@@ -185,7 +185,7 @@ rightSVAux f st x = unsafePerformIO $ do
     s <- createVector q
     v <- createMatrix ColumnMajor c c
     app3 g mat x vec s mat v st
-    return (s,trans v)
+    return (s,v)
   where r = rows x
         c = cols x
         q = min r c
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 3398e6a..25b10b4 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -1,5 +1,3 @@
-#if __GLASGOW_HASKELL__ >= 708
-
 {-# LANGUAGE DataKinds #-}
 {-# LANGUAGE KindSignatures #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
@@ -13,7 +11,6 @@
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE ViewPatterns #-}
 {-# LANGUAGE GADTs #-}
-{-# LANGUAGE OverlappingInstances #-}
 {-# LANGUAGE TypeFamilies #-}
 
 
@@ -63,13 +60,13 @@ module Numeric.LinearAlgebra.Static(
 
 
 import GHC.TypeLits
-import Numeric.LinearAlgebra.HMatrix hiding (
+import Numeric.LinearAlgebra hiding (
     (<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),
     row,col,vector,matrix,linspace,toRows,toColumns,
     (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
     eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
     qr,size,app,mul,dot,chol)
-import qualified Numeric.LinearAlgebra.HMatrix as LA
+import qualified Numeric.LinearAlgebra as LA
 import Data.Proxy(Proxy)
 import Numeric.LinearAlgebra.Static.Internal
 import Control.Arrow((***))
@@ -184,8 +181,9 @@ a ¦ b = tr (tr a —— tr b)
 type Sq n  = L n n
 --type CSq n = CL n n
 
-type GL = forall n m. (KnownNat n, KnownNat m) => L m n
-type GSq = forall n. KnownNat n => Sq n
+
+type GL = forall n m . (KnownNat n, KnownNat m) => L m n
+type GSq = forall n . KnownNat n => Sq n
 
 isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (ℝ,(Int,Int))
 isKonst s@(unwrap -> x)
@@ -618,23 +616,3 @@ instance (KnownNat n', KnownNat m') => Testable (L n' m')
   where
     checkT _ = test
 
-#else
-
-{- |
-Module      :  Numeric.LinearAlgebra.Static
-Copyright   :  (c) Alberto Ruiz 2014
-License     :  BSD3
-Stability   :  experimental
-
-Experimental interface with statically checked dimensions.
-
-This module requires GHC >= 7.8
-
--}
-
-module Numeric.LinearAlgebra.Static
-{-# WARNING "This module requires GHC >= 7.8" #-}
-where
-
-#endif
-
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
index ec02cf6..a5fc29b 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
@@ -24,8 +24,8 @@ module Numeric.LinearAlgebra.Static.Internal where
 
 
 import GHC.TypeLits
-import qualified Numeric.LinearAlgebra.HMatrix as LA
-import Numeric.LinearAlgebra.HMatrix hiding (konst,size)
+import qualified Numeric.LinearAlgebra as LA
+import Numeric.LinearAlgebra hiding (konst,size)
 import Data.Packed as D
 import Data.Packed.ST
 import Data.Proxy(Proxy)
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 043aa21..370ca27 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -16,7 +16,6 @@ Stability   :  provisional
 
 -}
 -----------------------------------------------------------------------------
-{-# OPTIONS_HADDOCK hide #-}
 
 module Numeric.LinearAlgebra.Util(
 
@@ -53,18 +52,7 @@ module Numeric.LinearAlgebra.Util(
     -- ** 1D
     corr, conv, corrMin,
     -- ** 2D
-    corr2, conv2, separable,
-    -- * Tools for the Kronecker product
-    --
-    -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in
-    --  3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132)
-
-    --
-    -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@
-    vec,
-    vech,
-    dup,
-    vtrans
+    corr2, conv2, separable
 ) where
 
 import Data.Packed.Numeric
@@ -227,10 +215,11 @@ infixl 9 ¿
 (¿)= flip extractColumns
 
 
-cross :: Vector Double -> Vector Double -> Vector Double
--- ^ cross product (for three-element real vectors)
+cross :: Product t => Vector t -> Vector t -> Vector t
+-- ^ cross product (for three-element vectors)
 cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
-          | otherwise = error $ "cross ("++show x++") ("++show y++")"
+          | otherwise = error $ "the cross product requires 3-element vectors (sizes given: "
+                                ++show (dim x)++" and "++show (dim y)++")"
   where
     [x1,x2,x3] = toList x
     [y1,y2,y3] = toList y
@@ -238,6 +227,9 @@ cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
     z2 = x3*y1-x1*y3
     z3 = x1*y2-x2*y1
 
+{-# SPECIALIZE cross :: Vector Double -> Vector Double -> Vector Double #-}
+{-# SPECIALIZE cross :: Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) #-}
+
 norm :: Vector Double -> Double
 -- ^ 2-norm of real vector
 norm = pnorm PNorm2
@@ -403,40 +395,6 @@ null1sym = last . toColumns . snd . eigSH'
 
 --------------------------------------------------------------------------------
 
-vec :: Element t => Matrix t -> Vector t
--- ^ stacking of columns
-vec = flatten . trans
-
-
-vech :: Element t => Matrix t -> Vector t
--- ^ half-vectorization (of the lower triangular part)
-vech m = vjoin . zipWith f [0..] . toColumns $ m
-  where
-    f k v = subVector k (dim v - k) v
-
-
-dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t
--- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k)
-dup k = trans $ fromRows $ map f es
-  where
-    rs = zip [0..] (toRows (ident (k^(2::Int))))
-    es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ]
-    f (i,j) | i == j = g (k*j + i)
-            | otherwise = g (k*j + i) + g (k*i + j)
-    g j = v
-      where
-        Just v = lookup j rs
-
-
-vtrans :: Element t => Int -> Matrix t -> Matrix t
--- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
-vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m
-           | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows"
-  where
-    (q,r) = divMod (rows m) p
-
---------------------------------------------------------------------------------
-
 infixl 0 ~!~
 c ~!~ msg = when c (error msg)
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index b82c74f..899a5bf 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -9,7 +9,7 @@ module Numeric.LinearAlgebra.Util.CG(
 import Data.Packed.Numeric
 import Numeric.Sparse
 import Numeric.Vector()
-import Numeric.LinearAlgebra.Algorithms(linearSolveLS, relativeError, NormType(..))
+import Numeric.LinearAlgebra.Algorithms(linearSolveLS, relativeError, pnorm, NormType(..))
 import Control.Arrow((***))
 
 {-
@@ -142,13 +142,13 @@ instance Testable GMatrix
             print s3; print d3
             print s4; print d4
             print s5; print d5
-            print $ relativeError Infinity s5 d5
+            print $ relativeError (pnorm Infinity) s5 d5
 
         ok = s1==d1
           && s2==d2
           && s3==d3
           && s4==d4
-          && relativeError Infinity s5 d5 < 1E-10
+          && relativeError (pnorm Infinity) s5 d5 < 1E-10
 
         disp = putStr . dispf 2
 
diff --git a/packages/base/src/Numeric/Vectorized.hs b/packages/base/src/Numeric/Vectorized.hs
index 6f0d240..405ae01 100644
--- a/packages/base/src/Numeric/Vectorized.hs
+++ b/packages/base/src/Numeric/Vectorized.hs
@@ -37,7 +37,6 @@ import Foreign.C.String
 import System.IO.Unsafe(unsafePerformIO)
 
 import Control.Monad(when)
-import Control.Applicative((<$>))
 
 
 
diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index 5a1b59c..8593e0a 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-glpk
-Version:            0.4.1.0
+Version:            0.5.0.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -23,7 +23,7 @@ extra-source-files:     examples/simplex1.hs
                         examples/simplex5.hs
 
 library
-    Build-Depends:      base <5, hmatrix >= 0.16, containers >= 0.5.4.0
+    Build-Depends:      base <5, hmatrix >= 0.17, containers
 
     hs-source-dirs:     src
 
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index 6a0c47d..7bf4279 100644
--- a/packages/glpk/src/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -85,8 +85,8 @@ module Numeric.LinearProgramming(
     Solution(..)
 ) where
 
-import Data.Packed
-import Data.Packed.Development
+import Numeric.LinearAlgebra.HMatrix
+import Numeric.LinearAlgebra.Devel hiding (Dense)
 import Foreign(Ptr)
 import System.IO.Unsafe(unsafePerformIO)
 import Foreign.C.Types
@@ -180,16 +180,17 @@ exact opt constr@(General _) bnds = exact opt (sparseOfGeneral constr) bnds
 
 adapt :: Optimization -> (Int, Double, [Double])
 adapt opt = case opt of
-    Maximize x -> (size x, 1 ,x)
-    Minimize x -> (size x, -1, (map negate x))
- where size x | null x = error "simplex: objective function with zero variables"
-              | otherwise = length x
+    Maximize x -> (sz x, 1 ,x)
+    Minimize x -> (sz x, -1, (map negate x))
+  where
+    sz x | null x = error "simplex: objective function with zero variables"
+         | otherwise = length x
 
 extract :: Double -> Vector Double -> Solution
 extract sg sol = r where
-    z = sg * (sol@>1)
-    v = toList $ subVector 2 (dim sol -2) sol
-    r = case round(sol@>0)::Int of
+    z = sg * (sol!1)
+    v = toList $ subVector 2 (size sol -2) sol
+    r = case round(sol!0)::Int of
           1 -> Undefined
           2 -> Feasible (z,v)
           3 -> Infeasible (z,v)
@@ -261,7 +262,7 @@ mkConstrD n f b1 | ok = fromLists (ob ++ co)
        ok = all (==n) ls
        den = fromLists cs
        ob = map (([0,0]++).return) f
-       co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
+       co = [[fromIntegral i, fromIntegral j,den `atIndex` (i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
 
 mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
 mkConstrS n objfun b1 = fromLists (ob ++ co) where
diff --git a/packages/glpk/src/Numeric/LinearProgramming/L1.hs b/packages/glpk/src/Numeric/LinearProgramming/L1.hs
index f55c721..d7f1258 100644
--- a/packages/glpk/src/Numeric/LinearProgramming/L1.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming/L1.hs
@@ -14,7 +14,7 @@ module Numeric.LinearProgramming.L1 (
     l1SolveU,
 ) where
 
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.LinearProgramming
 
 -- | L_inf solution of overconstrained system Ax=b.
diff --git a/packages/gsl/CHANGELOG b/packages/gsl/CHANGELOG
new file mode 100644
index 0000000..b83973a
--- /dev/null
+++ b/packages/gsl/CHANGELOG
@@ -0,0 +1,10 @@
+0.17.0.0
+--------
+
+    * Added interpolation modules
+
+0.16.0.0
+--------
+
+    * The modules Numeric.GSL.* have been moved from hmatrix to the new package hmatrix-gsl.
+
diff --git a/packages/gsl/THANKS.md b/packages/gsl/THANKS.md
new file mode 100644
index 0000000..b066185
--- /dev/null
+++ b/packages/gsl/THANKS.md
@@ -0,0 +1,4 @@
+Matt Peddie wrote the interpolation interface.
+
+(see also the THANKS file of the hmatrix package)
+
diff --git a/packages/gsl/hmatrix-gsl.cabal b/packages/gsl/hmatrix-gsl.cabal
index 2f6f51b..7028bae 100644
--- a/packages/gsl/hmatrix-gsl.cabal
+++ b/packages/gsl/hmatrix-gsl.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-gsl
-Version:            0.16.0.3
+Version:            0.17.0.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -25,7 +25,7 @@ flag onlygsl
 
 library
 
-    Build-Depends:      base<5, hmatrix>=0.16, array, vector,
+    Build-Depends:      base<5, hmatrix>=0.17, array, vector,
                         process, random
 
 
diff --git a/packages/gsl/src/Graphics/Plot.hs b/packages/gsl/src/Graphics/Plot.hs
index 0ea41ac..d2ea192 100644
--- a/packages/gsl/src/Graphics/Plot.hs
+++ b/packages/gsl/src/Graphics/Plot.hs
@@ -27,13 +27,13 @@ module Graphics.Plot(
 
 ) where
 
-import Numeric.Container
+import Numeric.LinearAlgebra.HMatrix
 import Data.List(intersperse)
 import System.Process (system)
 
 -- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions.
 meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)
-meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2)
+meshdom r1 r2 = (outer r1 (konst 1 (size r2)), outer (konst 1 (size r1)) r2)
 
 
 {- | Draws a 3D surface representation of a real matrix.
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
index 0a92373..db9d82f 100644
--- a/packages/gsl/src/Numeric/GSL/Fitting.hs
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE FlexibleContexts #-}
+
 {- |
 Module      :  Numeric.GSL.Fitting
 Copyright   :  (c) Alberto Ruiz 2010
@@ -50,7 +52,7 @@ module Numeric.GSL.Fitting (
     fitModelScaled, fitModel
 ) where
 
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL.Internal
 
 import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
@@ -80,13 +82,13 @@ nlFitting :: FittingMethod
 nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
 
 nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
-    let p   = dim xiv
-        n   = dim (f xiv)
+    let p   = size xiv
+        n   = size (f xiv)
     fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
     jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
     rawpath <- createMatrix RowMajor maxit (2+p)
     app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toRows $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -99,7 +101,7 @@ foreign import ccall safe "nlfit"
 -------------------------------------------------------
 
 checkdim1 n _p v
-    | dim v == n = v
+    | size v == n = v
     | otherwise = error $ "Error: "++ show n
                         ++ " components expected in the result of the function supplied to nlFitting"
 
@@ -114,9 +116,9 @@ err (model,deriv) dat vsol = zip sol errs where
     sol = toList vsol
     c = max 1 (chi/sqrt (fromIntegral dof))
     dof = length dat - (rows cov)
-    chi = norm2 (fromList $ cost (resMs model) dat sol)
+    chi = norm_2 (fromList $ cost (resMs model) dat sol)
     js = fromLists $ jacobian (resDs deriv) dat sol
-    cov = inv $ trans js <> js
+    cov = inv $ tr js <> js
     errs = toList $ scalar c * sqrt (takeDiag cov)
 
 
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
index 734325b..d824b4f 100644
--- a/packages/gsl/src/Numeric/GSL/Fourier.hs
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -16,14 +16,13 @@ module Numeric.GSL.Fourier (
     ifft
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL.Internal
-import Data.Complex
 import Foreign.C.Types
 import System.IO.Unsafe (unsafePerformIO)
 
 genfft code v = unsafePerformIO $ do
-    r <- createVector (dim v)
+    r <- createVector (size v)
     app2 (c_fft code) vec v vec r "fft"
     return r
 
diff --git a/packages/gsl/src/Numeric/GSL/IO.hs b/packages/gsl/src/Numeric/GSL/IO.hs
index 0d6031a..936f6bf 100644
--- a/packages/gsl/src/Numeric/GSL/IO.hs
+++ b/packages/gsl/src/Numeric/GSL/IO.hs
@@ -14,7 +14,7 @@ module Numeric.GSL.IO (
     fileDimensions, loadMatrix, fromFile
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix hiding(saveMatrix, loadMatrix)
 import Numeric.GSL.Vector
 import System.Process(readProcess)
 
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
index a1c4e0c..a269224 100644
--- a/packages/gsl/src/Numeric/GSL/Internal.hs
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -22,14 +22,13 @@ module Numeric.GSL.Internal(
     aux_vTom,
     createV,
     createMIO,
-    module Data.Packed.Development,
+    module Numeric.LinearAlgebra.Devel,
     check,
     Res,TV,TM,TCV,TCM
 ) where
 
-import Data.Packed
-import Data.Packed.Development hiding (check)
-import Data.Complex
+import Numeric.LinearAlgebra.HMatrix
+import Numeric.LinearAlgebra.Devel hiding (check)
 
 import Foreign.Marshal.Array(copyArray)
 import Foreign.Ptr(Ptr, FunPtr)
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
index 17e2258..cb78bf4 100644
--- a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -15,7 +15,7 @@ module Numeric.GSL.LinearAlgebra (
     fileDimensions, loadMatrix, fromFile
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix hiding (RandDist,randomVector,saveMatrix,loadMatrix)
 import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
 
 import Foreign.Marshal.Alloc(free)
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
index 056d463..00e0619 100644
--- a/packages/gsl/src/Numeric/GSL/Minimization.hs
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+
 {- |
 Module      :  Numeric.GSL.Minimization
 Copyright   :  (c) Alberto Ruiz 2006-9
@@ -56,7 +59,7 @@ module Numeric.GSL.Minimization (
 ) where
 
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix hiding(step)
 import Numeric.GSL.Internal
 
 import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
@@ -99,7 +102,7 @@ uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
     rawpath <- createMIO maxit 4
                          (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
                          "uniMinimize"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -137,13 +140,13 @@ minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList s
 ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
 
 minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
-    let n   = dim xiv
+    let n   = size xiv
     fp <- mkVecfun (iv f)
     rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
                    createMIO maxit (n+3)
                          (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
                          "minimize"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         sol = flatten $ dropColumns 3 $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -191,7 +194,7 @@ minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
 
 
 minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
-    let n = dim xiv
+    let n = size xiv
         f' = f
         df' = (checkdim1 n . df)
     fp <- mkVecfun (iv f')
@@ -200,7 +203,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
                     createMIO maxit (n+2)
                          (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
                          "minimizeD"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         sol = flatten $ dropColumns 2 $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -217,6 +220,6 @@ foreign import ccall safe "gsl-aux.h minimizeD"
 ---------------------------------------------------------------------
 
 checkdim1 n v
-    | dim v == n = v
+    | size v == n = v
     | otherwise = error $ "Error: "++ show n
                         ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 7549a65..3258b83 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+
 {- |
 Module      :  Numeric.GSL.ODE
 Copyright   :  (c) Alberto Ruiz 2010
@@ -32,7 +35,7 @@ module Numeric.GSL.ODE (
     odeSolve, odeSolveV, ODEMethod(..), Jacobian
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL.Internal
 
 import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
@@ -68,7 +71,7 @@ odeSolve
     -> Vector Double   -- ^ desired solution times
     -> Matrix Double   -- ^ solution
 odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
-    where hi = (ts@>1 - ts@>0)/100
+    where hi = (ts!1 - ts!0)/100
           epsAbs = 1.49012e-08
           epsRel = 1.49012e-08
           l2v f = \t -> fromList  . f t . toList
@@ -107,14 +110,14 @@ odeSolveV'
     -> Vector Double     -- ^ desired solution times
     -> Matrix Double     -- ^ solution
 odeSolveV' method mbjac h epsAbs epsRel f  xiv ts = unsafePerformIO $ do
-    let n   = dim xiv
+    let n   = size xiv
     fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
     jp <- case mbjac of
         Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
         Nothing  -> return nullFunPtr
     sol <- vec xiv $ \xiv' ->
             vec (checkTimes ts) $ \ts' ->
-             createMIO (dim ts) n
+             createMIO (size ts) n
               (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
               "ode"
     freeHaskellFunPtr fp
@@ -126,7 +129,7 @@ foreign import ccall safe "ode"
 -------------------------------------------------------
 
 checkdim1 n v
-    | dim v == n = v
+    | size v == n = v
     | otherwise = error $ "Error: "++ show n
                         ++ " components expected in the result of the function supplied to odeSolve"
 
@@ -135,6 +138,7 @@ checkdim2 n m
     | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
                         ++ " Jacobian expected in odeSolve"
 
-checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
+checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
               | otherwise = error "odeSolve requires increasing times"
     where ts' = toList ts
+
diff --git a/packages/gsl/src/Numeric/GSL/Polynomials.hs b/packages/gsl/src/Numeric/GSL/Polynomials.hs
index b1be85d..246e301 100644
--- a/packages/gsl/src/Numeric/GSL/Polynomials.hs
+++ b/packages/gsl/src/Numeric/GSL/Polynomials.hs
@@ -16,9 +16,8 @@ module Numeric.GSL.Polynomials (
     polySolve
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL.Internal
-import Data.Complex
 import System.IO.Unsafe (unsafePerformIO)
 
 #if __GLASGOW_HASKELL__ >= 704
@@ -47,8 +46,8 @@ polySolve :: [Double] -> [Complex Double]
 polySolve = toList . polySolve' . fromList
 
 polySolve' :: Vector Double -> Vector (Complex Double)
-polySolve' v | dim v > 1 = unsafePerformIO $ do
-    r <- createVector (dim v-1)
+polySolve' v | size v > 1 = unsafePerformIO $ do
+    r <- createVector (size v-1)
     app2 c_polySolve vec v vec r "polySolve"
     return r
              | otherwise = error "polySolve on a polynomial of degree zero"
diff --git a/packages/gsl/src/Numeric/GSL/Random.hs b/packages/gsl/src/Numeric/GSL/Random.hs
index f1f49e5..139c921 100644
--- a/packages/gsl/src/Numeric/GSL/Random.hs
+++ b/packages/gsl/src/Numeric/GSL/Random.hs
@@ -21,11 +21,13 @@ module Numeric.GSL.Random (
 ) where
 
 import Numeric.GSL.Vector
-import Numeric.LinearAlgebra(cholSH)
-import Numeric.Container hiding (
+import Numeric.LinearAlgebra.HMatrix hiding (
     randomVector,
     gaussianSample,
-    uniformSample
+    uniformSample,
+    Seed,
+    rand,
+    randn
     )
 import System.Random(randomIO)
 
@@ -40,10 +42,10 @@ gaussianSample :: Seed
                -> Matrix Double -- ^ covariance matrix
                -> Matrix Double -- ^ result
 gaussianSample seed n med cov = m where
-    c = dim med
+    c = size med
     meds = konst 1 n `outer` med
     rs = reshape c $ randomVector seed Gaussian (c * n)
-    m = rs `mXm` cholSH cov `add` meds
+    m = rs <> cholSH cov + meds
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- uniform distribution.
@@ -55,10 +57,10 @@ uniformSample seed n rgs = m where
     (as,bs) = unzip rgs
     a = fromList as
     cs = zipWith subtract as bs
-    d = dim a
+    d = size a
     dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
     am = konst 1 n `outer` a
-    m = fromColumns (zipWith scale cs dat) `add` am
+    m = fromColumns (zipWith scale cs dat) + am
 
 -- | pseudorandom matrix with uniform elements between 0 and 1
 randm :: RandDist
diff --git a/packages/gsl/src/Numeric/GSL/Root.hs b/packages/gsl/src/Numeric/GSL/Root.hs
index b9f3b94..724f32f 100644
--- a/packages/gsl/src/Numeric/GSL/Root.hs
+++ b/packages/gsl/src/Numeric/GSL/Root.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE  FlexibleContexts #-}
+
 {- |
 Module      :  Numeric.GSL.Root
 Copyright   :  (c) Alberto Ruiz 2009
@@ -39,7 +41,7 @@ module Numeric.GSL.Root (
     rootJ, RootMethodJ(..),
 ) where
 
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL.Internal
 import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
 import Foreign.C.Types
@@ -69,7 +71,7 @@ uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
     rawpath <- createMIO maxit 4
                          (c_root m fp epsrel (fi maxit) xl xu)
                          "root"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -100,7 +102,7 @@ uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
     rawpath <- createMIO maxit 2
                          (c_rootj m fp dfp epsrel (fi maxit) x)
                          "rootj"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -132,13 +134,13 @@ root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit ep
 
 rootGen m f xi epsabs maxit = unsafePerformIO $ do
     let xiv = fromList xi
-        n   = dim xiv
+        n   = size xiv
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
                          (c_multiroot m fp epsabs (fi maxit) // xiv')
                          "multiroot"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -169,14 +171,14 @@ rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun ja
 
 rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
     let xiv = fromList xi
-        n   = dim xiv
+        n   = size xiv
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
     jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
                          (c_multirootj m fp jp epsabs (fi maxit) // xiv')
                          "multiroot"
-    let it = round (rawpath @@> (maxit-1,0))
+    let it = round (rawpath `atIndex` (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
     freeHaskellFunPtr fp
@@ -189,7 +191,7 @@ foreign import ccall safe "multirootj"
 -------------------------------------------------------
 
 checkdim1 n v
-    | dim v == n = v
+    | size v == n = v
     | otherwise = error $ "Error: "++ show n
                         ++ " components expected in the result of the function supplied to root"
 
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
index af79f32..0cd99eb 100644
--- a/packages/gsl/src/Numeric/GSL/Vector.hs
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -14,8 +14,7 @@ module Numeric.GSL.Vector (
     fwriteVector, freadVector, fprintfVector, fscanfVector
 ) where
 
-import Data.Packed
-import Numeric.LinearAlgebra(RandDist(..))
+import Numeric.LinearAlgebra.HMatrix hiding(randomVector, saveMatrix)
 import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
 
 import Foreign.Marshal.Alloc(free)
diff --git a/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
index 8608394..d75fec2 100644
--- a/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
+++ b/packages/sparse/src/Numeric/LinearAlgebra/Sparse.hs
@@ -13,8 +13,11 @@ import System.IO.Unsafe(unsafePerformIO)
 import Foreign(Ptr)
 import Numeric.LinearAlgebra.HMatrix
 import Text.Printf
-import Numeric.LinearAlgebra.Util((~!~))
+import Control.Monad(when)
 
+(???) :: Bool -> String -> IO ()
+infixl 0 ???
+c ??? msg = when c (error msg)
 
 type IV t = CInt -> Ptr CInt   -> t
 type  V t = CInt -> Ptr Double -> t
@@ -22,7 +25,7 @@ type SMxV = V (IV (IV (V (V (IO CInt)))))
 
 dss :: CSR -> Vector Double -> Vector Double
 dss CSR{..} b = unsafePerformIO $ do
-    size b /= csrNRows ~!~ printf "dss: incorrect sizes: (%d,%d) x %d" csrNRows csrNCols (size b)
+    size b /= csrNRows ??? printf "dss: incorrect sizes: (%d,%d) x %d" csrNRows csrNCols (size b)
     r <- createVector csrNCols
     app5 c_dss vec csrVals vec csrCols vec csrRows vec b vec r "dss"
     return r
diff --git a/packages/special/hmatrix-special.cabal b/packages/special/hmatrix-special.cabal
index 28b294b..3b122c8 100644
--- a/packages/special/hmatrix-special.cabal
+++ b/packages/special/hmatrix-special.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-special
-Version:            0.3.0.1
+Version:            0.4.0.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -27,7 +27,7 @@ flag safe-cheap
     default:        False
 
 library
-    Build-Depends:  base <5, hmatrix, hmatrix-gsl
+    Build-Depends:  base <5, hmatrix>=0.17, hmatrix-gsl
 
     Extensions:     ForeignFunctionInterface,
                     CPP
diff --git a/packages/special/lib/Numeric/GSL/Special/Internal.hsc b/packages/special/lib/Numeric/GSL/Special/Internal.hsc
index e7c38e8..a9aab9b 100644
--- a/packages/special/lib/Numeric/GSL/Special/Internal.hsc
+++ b/packages/special/lib/Numeric/GSL/Special/Internal.hsc
@@ -33,7 +33,7 @@ import Foreign.Storable
 import Foreign.Ptr
 import Foreign.Marshal
 import System.IO.Unsafe(unsafePerformIO)
-import Data.Packed.Development(check,(//))
+import Numeric.LinearAlgebra.Devel(check,(//))
 import Foreign.C.Types
 
 data Precision = PrecDouble | PrecSingle | PrecApprox
diff --git a/packages/tests/hmatrix-tests.cabal b/packages/tests/hmatrix-tests.cabal
index 0514843..de796e8 100644
--- a/packages/tests/hmatrix-tests.cabal
+++ b/packages/tests/hmatrix-tests.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-tests
-Version:            0.4.1.0
+Version:            0.5.0.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -28,9 +28,9 @@ library
 
     Build-Depends:      base >= 4 && < 5,
                         QuickCheck >= 2, HUnit, random,
-                        hmatrix >= 0.16
+                        hmatrix >= 0.17
     if flag(gsl)
-      Build-Depends:    hmatrix-gsl >= 0.16
+      Build-Depends:    hmatrix-gsl >= 0.17
 
     hs-source-dirs:     src
 
diff --git a/packages/tests/src/Numeric/GSL/Tests.hs b/packages/tests/src/Numeric/GSL/Tests.hs
index 9dff6f5..e5d205d 100644
--- a/packages/tests/src/Numeric/GSL/Tests.hs
+++ b/packages/tests/src/Numeric/GSL/Tests.hs
@@ -19,7 +19,7 @@ import System.Exit (exitFailure)
 
 import Test.HUnit (runTestTT, failures, Test(..), errors)
 
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.HMatrix
 import Numeric.GSL
 import Numeric.LinearAlgebra.Tests (qCheck, utest)
 import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~))
@@ -42,7 +42,7 @@ fittingTest = utest "levmar" (ok1 && ok2)
     sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
 
     ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d
-    ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5
+    ok2 = norm_2 (fromList (map fst sols) - fromList sol) < 1E-5
 
 ---------------------------------------------------------------------
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 713af79..71c7c16 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -28,12 +28,9 @@ module Numeric.LinearAlgebra.Tests(
 --, runBigTests
 ) where
 
-import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.HMatrix hiding ((<>),linearSolve)
+import Numeric.LinearAlgebra hiding (unitary)
+import Numeric.LinearAlgebra.Devel hiding (vec)
 import Numeric.LinearAlgebra.Static(L)
-import Numeric.LinearAlgebra.Util(col,row)
-import Data.Packed
-import Numeric.LinearAlgebra.LAPACK
 import Numeric.LinearAlgebra.Tests.Instances
 import Numeric.LinearAlgebra.Tests.Properties
 import Test.HUnit hiding ((~:),test,Testable,State)
@@ -44,16 +41,13 @@ import qualified Prelude
 import System.CPUTime
 import System.Exit
 import Text.Printf
-import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)
+import Numeric.LinearAlgebra.Devel(unsafeFromForeignPtr,unsafeToForeignPtr)
 import Control.Arrow((***))
 import Debug.Trace
 import Control.Monad(when)
-import Numeric.LinearAlgebra.Util hiding (ones,row,col)
 import Control.Applicative
 import Control.Monad(ap)
 
-import Data.Packed.ST
-
 import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
                       ,sized,classify,Testable,Property
                       ,quickCheckWithResult,maxSize,stdArgs,shrink)
@@ -89,7 +83,7 @@ detTest1 = det m == 26
         mc = (3><3)
             [ 1, 2, 3
             , 4, 5, 7
-            , 2, 8, i
+            , 2, 8, iC
             ]
 
 detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]
@@ -140,7 +134,7 @@ randomTestGaussian = c :~1~: snd (meanCov dat) where
                 2,4,0,
                -2,2,1]
     m = 3 |> [1,2,3]
-    c = a <> trans a
+    c = a <> tr a
     dat = gaussianSample 7 (10^6) m c
 
 randomTestUniform = c :~1~: snd (meanCov dat) where
@@ -174,54 +168,54 @@ offsetTest = y == y' where
 
 normsVTest = TestList [
     utest "normv2CD" $ norm2PropC v
-  , utest "normv2CF" $ norm2PropC (single v)
+--  , utest "normv2CF" $ norm2PropC (single v)
 #ifndef NONORMVTEST
   , utest "normv2D"  $ norm2PropR x
-  , utest "normv2F"  $ norm2PropR (single x)
+--  , utest "normv2F"  $ norm2PropR (single x)
 #endif
-  , utest "normv1CD" $ norm1 v          == 8
-  , utest "normv1CF" $ norm1 (single v) == 8
-  , utest "normv1D"  $ norm1 x          == 6
-  , utest "normv1F"  $ norm1 (single x) == 6
+  , utest "normv1CD" $ norm_1 v          == 8
+--  , utest "normv1CF" $ norm_1 (single v) == 8
+  , utest "normv1D"  $ norm_1 x          == 6
+--  , utest "normv1F"  $ norm_1 (single x) == 6
 
-  , utest "normvInfCD" $ normInf v          == 5
-  , utest "normvInfCF" $ normInf (single v) == 5
-  , utest "normvInfD"  $ normInf x          == 3
-  , utest "normvInfF"  $ normInf (single x) == 3
+  , utest "normvInfCD" $ norm_Inf v          == 5
+--  , utest "normvInfCF" $ norm_Inf (single v) == 5
+  , utest "normvInfD"  $ norm_Inf x          == 3
+--  , utest "normvInfF"  $ norm_Inf (single x) == 3
 
  ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)
          x = fromList [1,2,-3] :: Vector Double
 #ifndef NONORMVTEST
-         norm2PropR a = norm2 a =~= sqrt (udot a a)
+         norm2PropR a = norm_2 a =~= sqrt (udot a a)
 #endif
-         norm2PropC a = norm2 a =~= realPart (sqrt (a <.> a))
+         norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a))
          a =~= b = fromList [a] |~| fromList [b]
 
 normsMTest = TestList [
-    utest "norm2mCD" $ pnorm PNorm2 v          =~= 8.86164970498005
-  , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005
-  , utest "norm2mD"  $ pnorm PNorm2 x          =~= 5.96667765076216
-  , utest "norm2mF"  $ pnorm PNorm2 (single x) =~= 5.96667765076216
-
-  , utest "norm1mCD" $ pnorm PNorm1 v          == 9
-  , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9
-  , utest "norm1mD"  $ pnorm PNorm1 x          == 7
-  , utest "norm1mF"  $ pnorm PNorm1 (single x) == 7
-
-  , utest "normmInfCD" $ pnorm Infinity v          == 12
-  , utest "normmInfCF" $ pnorm Infinity (single v) == 12
-  , utest "normmInfD"  $ pnorm Infinity x          == 8
-  , utest "normmInfF"  $ pnorm Infinity (single x) == 8
-
-  , utest "normmFroCD" $ pnorm Frobenius v          =~= 8.88819441731559
-  , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559
-  , utest "normmFroD"  $ pnorm Frobenius x          =~= 6.24499799839840
-  , utest "normmFroF"  $ pnorm Frobenius (single x) =~~= 6.24499799839840
-
- ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double)
+    utest "norm2mCD" $ norm_2 v          =~= 8.86164970498005
+--  , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005
+  , utest "norm2mD"  $ norm_2 x          =~= 5.96667765076216
+--  , utest "norm2mF"  $ norm_2 (single x) =~= 5.96667765076216
+
+  , utest "norm1mCD" $ norm_1 v          == 9
+--  , utest "norm1mCF" $ norm_1 (single v) == 9
+  , utest "norm1mD"  $ norm_1 x          == 7
+--  , utest "norm1mF"  $ norm_1 (single x) == 7
+
+  , utest "normmInfCD" $ norm_Inf v          == 12
+--  , utest "normmInfCF" $ norm_Inf (single v) == 12
+  , utest "normmInfD"  $ norm_Inf x          == 8
+--  , utest "normmInfF"  $ norm_Inf (single x) == 8
+
+  , utest "normmFroCD" $ norm_Frob v          =~= 8.88819441731559
+--  , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559
+  , utest "normmFroD"  $ norm_Frob x          =~= 6.24499799839840
+--  , utest "normmFroF"  $ norm_Frob (single x) =~~= 6.24499799839840
+
+ ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Double)
          x = (2><2) [1,2,-3,5] :: Matrix Double
          a =~= b = fromList [a] :~10~: fromList [b]
-         a =~~= b = fromList [a] :~5~: fromList [b]
+--       a =~~= b = fromList [a] :~5~: fromList [b]
 
 ---------------------------------------------------------------------
 
@@ -236,7 +230,7 @@ sumprodTest = TestList [
   , utest "prodD"  $ prodProp v
   , utest "prodF"  $ prodProp (single v)
  ] where v = fromList [1,2,3] :: Vector Double
-         z = fromList [1,2-i,3+i]
+         z = fromList [1,2-iC,3+iC]
          prodProp x = prodElements x == product (toList x)
 
 ---------------------------------------------------------------------
@@ -250,7 +244,7 @@ chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where
 
 ---------------------------------------------------------------------
 
-conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)
+conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m)
 
 ---------------------------------------------------------------------
 
@@ -306,7 +300,7 @@ lift_maybe m = MaybeT $ do
 
 -- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
 --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
-successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0)
+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0)
    where stp e  = do
                   ep <- lift_maybe $ state_get
                   if t e ep
@@ -315,7 +309,7 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s
 
 -- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
 --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
-successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0)
+successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0)
    where stp  e = do
                   ep <- state_get
                   state_put e
@@ -377,23 +371,6 @@ convolutionTest = utest "convolution" ok
 
 --------------------------------------------------------------------------------
 
-kroneckerTest = utest "kronecker" ok
-  where
-    a,x,b :: Matrix Double
-    a = (3><4) [1..]
-    x = (4><2) [3,5..]
-    b = (2><5) [0,5..]
-    v1 = vec (a <> x <> b)
-    v2 = (trans b `kronecker` a) <> vec x
-    s = trans b <> b
-    v3 = vec s
-    v4 = (dup 5 :: Matrix Double) <> vech s
-    ok = v1 == v2 && v3 == v4
-      && vtrans 1 a == trans a
-      && vtrans (rows a) a == asColumn (vec a)
-
---------------------------------------------------------------------------------
-
 sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))
 
 --------------------------------------------------------------------------------
@@ -435,11 +412,11 @@ runTests n = do
     test (multProp1 10 . cConsist)
     test (multProp2 10 . rConsist)
     test (multProp2 10 . cConsist)
-    putStrLn "------ mult Float"
-    test (multProp1  6 . (single *** single) . rConsist)
-    test (multProp1  6 . (single *** single) . cConsist)
-    test (multProp2  6 . (single *** single) . rConsist)
-    test (multProp2  6 . (single *** single) . cConsist)
+--    putStrLn "------ mult Float"
+--    test (multProp1  6 . (single *** single) . rConsist)
+--    test (multProp1  6 . (single *** single) . cConsist)
+--    test (multProp2  6 . (single *** single) . rConsist)
+--    test (multProp2  6 . (single *** single) . cConsist)
     putStrLn "------ sub-trans"
     test (subProp . rM)
     test (subProp . cM)
@@ -472,16 +449,16 @@ runTests n = do
     putStrLn "------ svd"
     test (svdProp1  . rM)
     test (svdProp1  . cM)
-    test (svdProp1a svdR)
-    test (svdProp1a svdC)
-    test (svdProp1a svdRd)
-    test (svdProp1b svdR)
-    test (svdProp1b svdC)
-    test (svdProp1b svdRd)
-    test (svdProp2 thinSVDR)
-    test (svdProp2 thinSVDC)
-    test (svdProp2 thinSVDRd)
-    test (svdProp2 thinSVDCd)
+    test (svdProp1a svd . rM)
+    test (svdProp1a svd . cM)
+--    test (svdProp1a svdRd)
+    test (svdProp1b svd . rM)
+    test (svdProp1b svd . cM)
+--    test (svdProp1b svdRd)
+    test (svdProp2 thinSVD . rM)
+    test (svdProp2 thinSVD . cM)
+--    test (svdProp2 thinSVDRd)
+--    test (svdProp2 thinSVDCd)
     test (svdProp3  . rM)
     test (svdProp3  . cM)
     test (svdProp4  . rM)
@@ -492,12 +469,12 @@ runTests n = do
     test (svdProp6b)
     test (svdProp7  . rM)
     test (svdProp7  . cM)
-    putStrLn "------ svdCd"
+--    putStrLn "------ svdCd"
 #ifdef NOZGESDD
-    putStrLn "Omitted"
+--    putStrLn "Omitted"
 #else
-    test (svdProp1a svdCd)
-    test (svdProp1b svdCd)
+--    test (svdProp1a svdCd)
+--    test (svdProp1b svdCd)
 #endif
     putStrLn "------ eig"
     test (eigSHProp . rHer)
@@ -515,10 +492,10 @@ runTests n = do
     test (qrProp     . rM)
     test (qrProp     . cM)
     test (rqProp     . rM)
-    test (rqProp     . cM)
+--    test (rqProp     . cM)
     test (rqProp1     . cM)
     test (rqProp2     . cM)
-    test (rqProp3     . cM)
+--    test (rqProp3     . cM)
     putStrLn "------ hess"
     test (hessProp   . rSq)
     test (hessProp   . cSq)
@@ -539,12 +516,12 @@ runTests n = do
     test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
     test (\u -> cos u * tan u |~| sin (u::RM))
     test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary
-    putStrLn "------ vector operations - Float"
-    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))
-    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary
-    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))
-    test (\u -> cos u * tan u |~~| sin (u::FM))
-    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary
+--    putStrLn "------ vector operations - Float"
+--    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))
+--    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary
+--    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))
+--    test (\u -> cos u * tan u |~~| sin (u::FM))
+--    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary
     putStrLn "------ read . show"
     test (\m -> (m::RM) == read (show m))
     test (\m -> (m::CM) == read (show m))
@@ -562,8 +539,8 @@ runTests n = do
         , utest "expm1" (expmTest1)
         , utest "expm2" (expmTest2)
         , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
-        , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM)
-        , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0
+        , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM)
+        , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0
         , utest "<\\>"   $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
 --        , utest "gamma" (gamma 5 == 24.0)
 --        , besselTest
@@ -571,10 +548,10 @@ runTests n = do
         , utest "randomGaussian" randomTestGaussian
         , utest "randomUniform" randomTestUniform
         , utest "buildVector/Matrix" $
-                        complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral
-                     && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)
-        , utest "rank" $  rank ((2><3)[1,0,0,1,5*eps,0]) == 1
-                       && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
+                        complex (10 |> [0::Double ..]) == build 10 id
+                     && ident 5 == build (5,5) (\r c -> if r==c then 1::Double else 0)
+        , utest "rank" $  rank ((2><3)[1,0,0,1,5*peps,0::Double]) == 1
+                       && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
         , mbCholTest
         , utest "offset" offsetTest
@@ -588,7 +565,6 @@ runTests n = do
         , conformTest
         , accumTest
         , convolutionTest
-        , kroneckerTest
         , sparseTest
         , staticTest
         ]
@@ -597,12 +573,12 @@ runTests n = do
 
 
 -- single precision approximate equality
-infixl 4 |~~|
-a |~~| b = a :~6~: b
+-- infixl 4 |~~|
+-- a |~~| b = a :~6~: b
 
 makeUnitary v | realPart n > 1    = v / scalar n
               | otherwise = v
-    where n = sqrt (v <.> v)
+    where n = sqrt (v `dot` v)
 
 -- -- | Some additional tests on big matrices. They take a few minutes.
 -- runBigTests :: IO ()
@@ -668,9 +644,9 @@ manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs
 
 
 manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs
-manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs
+manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs
 
-manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs
+manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs
 
 vec3 :: Double -> Double -> Double -> Vector Double
 vec3 a b c = runSTVector $ do
@@ -695,11 +671,11 @@ mkVecBench = do
 
 subBench = do
     putStrLn ""
-    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))
-    time "0.1M subVector   " (g (konst 1 (1+10^5) :: Vector Double) @> 0)
+    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v))
+    time "0.1M subVector   " (g (konst 1 (1+10^5) :: Vector Double) ! 0)
     let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
-    time "subVector-join  3" (f (ident  3 :: Matrix Double) @@>(0,0))
-    time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))
+    time "subVector-join  3" (f (ident  3 :: Matrix Double) `atIndex` (0,0))
+    time "subVector-join 10" (f (ident 10 :: Matrix Double) `atIndex` (0,0))
 
 --------------------------------
 
@@ -724,7 +700,7 @@ multBench = do
 
 eigBench = do
     let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
-        s = m + trans m
+        s = m + tr m
     m `seq` s `seq` putStrLn ""
     time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH' m)
     time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)
@@ -736,7 +712,7 @@ eigBench = do
 svdBench = do
     let a = reshape 500  (randomVector 777 Uniform (3000*500))
         b = reshape 1000 (randomVector 777 Uniform (1000*1000))
-        fv (_,_,v) = v@@>(0,0)
+        fv (_,_,v) = v `atIndex` (0,0)
     a `seq` b `seq` putStrLn ""
     time "singular values  3000x500" (singularValues a)
     time "thin svd         3000x500" (fv $ thinSVD a)
@@ -748,7 +724,7 @@ svdBench = do
 
 solveBenchN n = do
     let x = uniformSample 777 (2*n) (replicate n (-1,1))
-        a = trans x <> x
+        a = tr x <> x
         b = asColumn $ randomVector 666 Uniform n
     a `seq` b `seq` putStrLn ""
     time ("svd solve " ++ show n) (linearSolveSVD a b)
@@ -765,7 +741,7 @@ solveBench = do
 
 cholBenchN n = do
     let x = uniformSample 777 (2*n) (replicate n (-1,1))
-        a = trans x <> x
+        a = tr x <> x
     a `seq` putStr ""
     time ("chol " ++ show n) (chol a)
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 53fc4d2..904ae05 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,5 +1,4 @@
-{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-}
-{-# OPTIONS_GHC -fno-warn-unused-imports #-}
+{-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Tests.Instances
@@ -26,15 +25,11 @@ module Numeric.LinearAlgebra.Tests.Instances(
 
 import System.Random
 
-import Numeric.LinearAlgebra
-import Numeric.LinearAlgebra.Devel
-import Numeric.Container
+import Numeric.LinearAlgebra.HMatrix hiding (vector)
 import Control.Monad(replicateM)
-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
-                      ,sized,classify,Testable,Property
-                      ,quickCheckWith,maxSize,stdArgs,shrink)
+import Test.QuickCheck(Arbitrary,arbitrary,choose,vector,sized,shrink)
+
 
-#if MIN_VERSION_QuickCheck(2,0,0)
 shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
 shrinkListElementwise []     = []
 shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]
@@ -42,25 +37,6 @@ shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]
 
 shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
 shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]
-#endif
-
-#if MIN_VERSION_QuickCheck(2,1,1)
-#else
-instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
-    arbitrary = do
-        re <- arbitrary
-        im <- arbitrary
-        return (re :+ im)
-
-#if MIN_VERSION_QuickCheck(2,0,0)
-    shrink (re :+ im) = 
-        [ u :+ v | (u,v) <- shrinkPair (re,im) ]
-#else
-    -- this has been moved to the 'Coarbitrary' class in QuickCheck 2
-    coarbitrary = undefined 
-#endif
-
-#endif
 
 chooseDim = sized $ \m -> choose (1,max 1 m)
 
@@ -68,15 +44,9 @@ instance (Field a, Arbitrary a) => Arbitrary (Vector a) where
     arbitrary = do m <- chooseDim
                    l <- vector m
                    return $ fromList l
-
-#if MIN_VERSION_QuickCheck(2,0,0)
     -- shrink any one of the components
     shrink = map fromList . shrinkListElementwise . toList
 
-#else
-    coarbitrary = undefined
-#endif
-
 instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 
     arbitrary = do
         m <- chooseDim
@@ -84,16 +54,11 @@ instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
         l <- vector (m*n)
         return $ (m><n) l
 
-#if MIN_VERSION_QuickCheck(2,0,0)
     -- shrink any one of the components
     shrink a = map (rows a >< cols a)
                . shrinkListElementwise
                . concat . toLists 
                      $ a
-#else
-    coarbitrary = undefined
-#endif
-
 
 -- a square matrix
 newtype (Sq a) = Sq (Matrix a) deriving Show
@@ -103,11 +68,7 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
         l <- vector (n*n)
         return $ Sq $ (n><n) l
 
-#if MIN_VERSION_QuickCheck(2,0,0)
     shrink (Sq a) = [ Sq b | b <- shrink a ]
-#else
-    coarbitrary = undefined
-#endif
 
 
 -- a unitary matrix
@@ -118,11 +79,6 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
         let (q,_) = qr m
         return (Rot q)
 
-#if MIN_VERSION_QuickCheck(2,0,0)
-#else
-    coarbitrary = undefined
-#endif
-
 
 -- a complex hermitian or real symmetric matrix
 newtype (Her a) = Her (Matrix a) deriving Show
@@ -130,12 +86,8 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
     arbitrary = do
         Sq m <- arbitrary
         let m' = m/2
-        return $ Her (m' + ctrans m')
+        return $ Her (m' + tr m')
 
-#if MIN_VERSION_QuickCheck(2,0,0)
-#else
-    coarbitrary = undefined
-#endif
 
 class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
 instance ArbitraryField Double
@@ -144,7 +96,7 @@ instance ArbitraryField (Complex Double)
 
 -- a well-conditioned general matrix (the singular values are between 1 and 100)
 newtype (WC a) = WC (Matrix a) deriving Show
-instance (ArbitraryField a) => Arbitrary (WC a) where
+instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where
     arbitrary = do
         m <- arbitrary
         let (u,_,v) = svd m
@@ -153,34 +105,24 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
             n = min r c
         sv' <- replicateM n (choose (1,100))
         let s = diagRect 0 (fromList sv') r c
-        return $ WC (u `mXm` real s `mXm` trans v)
-
-#if MIN_VERSION_QuickCheck(2,0,0)
-#else
-    coarbitrary = undefined
-#endif
+        return $ WC (u <> real s <> tr v)
 
 
 -- a well-conditioned square matrix (the singular values are between 1 and 100)
 newtype (SqWC a) = SqWC (Matrix a) deriving Show
-instance (ArbitraryField a) => Arbitrary (SqWC a) where
+instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where
     arbitrary = do
         Sq m <- arbitrary
         let (u,_,v) = svd m
             n = rows m
         sv' <- replicateM n (choose (1,100))
         let s = diag (fromList sv')
-        return $ SqWC (u `mXm` real s `mXm` trans v)
-
-#if MIN_VERSION_QuickCheck(2,0,0)
-#else
-    coarbitrary = undefined
-#endif
+        return $ SqWC (u <> real s <> tr v)
 
 
 -- a positive definite square matrix (the eigenvalues are between 0 and 100)
 newtype (PosDef a) = PosDef (Matrix a) deriving Show
-instance (ArbitraryField a, Num (Vector a)) 
+instance (Numeric a, ArbitraryField a, Num (Vector a)) 
     => Arbitrary (PosDef a) where
     arbitrary = do
         Her m <- arbitrary
@@ -188,13 +130,8 @@ instance (ArbitraryField a, Num (Vector a))
             n = rows m
         l <- replicateM n (choose (0,100))
         let s = diag (fromList l)
-            p = v `mXm` real s `mXm` ctrans v
-        return $ PosDef (0.5 * p + 0.5 * ctrans p)
-
-#if MIN_VERSION_QuickCheck(2,0,0)
-#else
-    coarbitrary = undefined
-#endif
+            p = v <> real s <> tr v
+        return $ PosDef (0.5 * p + 0.5 * tr p)
 
 
 -- a pair of matrices that can be multiplied
@@ -208,11 +145,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
         lb <- vector (k*m)
         return $ Consistent ((n><k) la, (k><m) lb)
 
-#if MIN_VERSION_QuickCheck(2,0,0)
     shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]
-#else
-    coarbitrary = undefined
-#endif
 
 
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index a5c37f4..207a303 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -1,5 +1,4 @@
-{-# LANGUAGE CPP, FlexibleContexts #-}
-{-# OPTIONS_GHC -fno-warn-unused-imports #-}
+{-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE TypeFamilies #-}
 
 -----------------------------------------------------------------------------
@@ -29,7 +28,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     pinvProp,
     detProp,
     nullspaceProp,
-    bugProp,
+--    bugProp,
     svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
     eigProp, eigSHProp, eigProp2, eigSHProp2,
@@ -43,20 +42,15 @@ module Numeric.LinearAlgebra.Tests.Properties (
     linearSolveProp, linearSolveProp2
 ) where
 
-import Numeric.Container
-import Numeric.LinearAlgebra --hiding (real,complex)
-import Numeric.LinearAlgebra.LAPACK
-import Debug.Trace
-import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
-                      ,sized,classify,Testable,Property
-                      ,quickCheckWith,maxSize,stdArgs,shrink)
+import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
+import Test.QuickCheck
 
 trivial :: Testable a => Bool -> a -> Property
 trivial = (`classify` "trivial")
 
 -- relative error
-dist :: (Normed c t, Num (c t)) => c t -> c t -> Double
-dist = relativeError Infinity
+dist :: (Num a, Normed a) => a -> a -> Double
+dist = relativeError norm_Inf
 
 infixl 4 |~|
 a |~| b = a :~10~: b
@@ -73,11 +67,11 @@ a :~n~: b = dist a b < 10^^(-n)
 square m = rows m == cols m
 
 -- orthonormal columns
-orthonormal m = ctrans m <> m |~| ident (cols m)
+orthonormal m = tr m <> m |~| ident (cols m)
 
 unitary m = square m && orthonormal m
 
-hermitian m = square m && m |~| ctrans m
+hermitian m = square m && m |~| tr m
 
 wellCond m = rcond m > 1/100
 
@@ -85,12 +79,12 @@ positiveDefinite m = minimum (toList e) > 0
     where (e,_v) = eigSH m
 
 upperTriang m = rows m == 1 || down == z
-    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
-          z = konst 0 (dim down)
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
+          z = konst 0 (size down)
 
 upperHessenberg m = rows m < 3 || down == z
-    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
-          z = konst 0 (dim down)
+    where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
+          z = konst 0 (size down)
 
 zeros (r,c) = reshape c (konst 0 (r*c))
 
@@ -118,81 +112,94 @@ detProp m = s d1 |~| s d2
           s x = fromList [x]
 
 nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
-                                     && orthonormal (fromColumns nl))
-    where nl = nullspacePrec 1 m
-          n = fromColumns nl
+                                     && orthonormal n)
+    where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
+          nl = toColumns n
           r = rows m
           c = cols m - rank m
 
 ------------------------------------------------------------------
-
+{-
 -- testcase for nonempty fpu stack
 -- uncommenting unitary' signature eliminates the problem
-bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v
-    where (u,d,v) = fullSVD m
+bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
+    where (u,d,v) = svd m
           -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
           unitary' a = unitary a
-
+-}
 ------------------------------------------------------------------
 
 -- fullSVD
-svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
-    where (u,d,v) = fullSVD m
+svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
+  where
+    (u,s,v) = svd m
+    d = diagRect 0 s (rows m) (cols m)
 
-svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where
+svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
+  where
     (u,s,v) = svdfun m
     d = diagRect 0 s (rows m) (cols m)
 
-svdProp1b svdfun m = unitary u && unitary v where
+svdProp1b svdfun m = unitary u && unitary v
+  where
     (u,_,v) = svdfun m
 
 -- thinSVD
-svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)
-    where (u,s,v) = thinSVDfun m
+svdProp2 thinSVDfun m
+    =  m |~| u <> diag (real s) <> tr v
+    && orthonormal u && orthonormal v
+    && size s == min (rows m) (cols m)
+  where
+    (u,s,v) = thinSVDfun m
 
 -- compactSVD
-svdProp3 m = (m |~| u <> real (diag s) <> trans v
+svdProp3 m = (m |~| u <> real (diag s) <> tr v
              && orthonormal u && orthonormal v)
-    where (u,s,v) = compactSVD m
+  where
+    (u,s,v) = compactSVD m
 
-svdProp4 m' = m |~| u <> real (diag s) <> trans v
+svdProp4 m' = m |~| u <> real (diag s) <> tr v
            && orthonormal u && orthonormal v
-           && (dim s == r || r == 0 && dim s == 1)
-    where (u,s,v) = compactSVD m
-          m = fromBlocks [[m'],[m']]
-          r = rank m'
-
-svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where
-    s1       = svR  m
-    s2       = svRd m
-    (_,s3,_) = svdR m
-    (_,s4,_) = svdRd m
-    (_,s5,_) = thinSVDR m
-    (_,s6,_) = thinSVDRd m
-
-svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where
-    s1       = svC  m
-    s2       = svCd m
-    (_,s3,_) = svdC m
-    (_,s4,_) = svdCd m
-    (_,s5,_) = thinSVDC m
-    (_,s6,_) = thinSVDCd m
+           && (size s == r || r == 0 && size s == 1)
+  where
+    (u,s,v) = compactSVD m
+    m = fromBlocks [[m'],[m']]
+    r = rank m'
+
+svdProp5a m = all (s1|~|) [s3,s5] where
+    s1       = singularValues (m :: Matrix Double)
+--  s2       = svRd m
+    (_,s3,_) = svd m
+--  (_,s4,_) = svdRd m
+    (_,s5,_) = thinSVD m
+--  (_,s6,_) = thinSVDRd m
+
+svdProp5b m = all (s1|~|) [s3,s5] where
+    s1       = singularValues (m :: Matrix (Complex Double))
+--  s2       = svCd m
+    (_,s3,_) = svd m
+--  (_,s4,_) = svdCd m
+    (_,s5,_) = thinSVD m
+--  (_,s6,_) = thinSVDCd m
 
 svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
-    where (u,s,v) = svdR m
-          (s',v') = rightSVR m
-          (u',s'') = leftSVR m
+  where
+    (u,s,v) = svd (m :: Matrix Double)
+    (s',v') = rightSV m
+    (u',s'') = leftSV m
 
 svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
-    where (u,s,v) = svdC m
-          (s',v') = rightSVC m
-          (u',s'') = leftSVC m
+  where
+    (u,s,v) = svd (m :: Matrix (Complex Double))
+    (s',v') = rightSV m
+    (u',s'') = leftSV m
 
 svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
-    where (u,s,v) = svd m
-          (s',v') = rightSV m
-          (u',_s'') = leftSV m
-          s''' = singularValues m
+  where
+    (u,s,v) = svd m
+    (s',v') = rightSV m
+    (u',_s'') = leftSV m
+    s''' = singularValues m
 
 ------------------------------------------------------------------
 
@@ -201,7 +208,7 @@ eigProp m = complex m <> v |~| v <> diag s
 
 eigSHProp m = m <> v |~| v <> real (diag s)
               && unitary v
-              && m |~| v <> real (diag s) <> ctrans v
+              && m |~| v <> real (diag s) <> tr v
     where (s, v) = eigSH m
 
 eigProp2 m = fst (eig m) |~| eigenvalues m
@@ -226,19 +233,19 @@ rqProp3 m = upperTriang' r
     where (r,_q) = rq m
 
 upperTriang' r = upptr (rows r) (cols r) * r |~| r
-    where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
-              where t = f-c
+    where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
+              where t = fromIntegral (f-c)
 
-hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
+hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
     where (p,h) = hess m
 
-schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
+schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
     where (u,s) = schur m
 
-schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
+schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
     where (u,s) = schur m
 
-cholProp m = m |~| ctrans c <> c && upperTriang c
+cholProp m = m |~| tr c <> c && upperTriang c
     where c = chol m
 
 exactProp m = chol m == chol (m+0)
@@ -252,7 +259,7 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
 
 multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
 
-multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)
+multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
 
 linearSolveProp f m = f m m |~| ident (rows m)
 
@@ -261,5 +268,5 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
           b = a <> x
           wc = rank a == q
 
-subProp m = m == (trans . fromColumns . toRows) m
+subProp m = m == (conj . tr . fromColumns . toRows) m