diff --git a/hTensor.cabal b/hTensor.cabal
--- a/hTensor.cabal
+++ b/hTensor.cabal
@@ -1,5 +1,5 @@
 Name:               hTensor
-Version:            0.8.2
+Version:            0.9.0
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -35,7 +35,7 @@
 
 library
 
-    Build-Depends:      base<5, hmatrix>=0.16, containers, random
+    Build-Depends:      base<5, hmatrix>=0.17, containers, random
 
     hs-source-dirs:     lib
     Exposed-modules:    Numeric.LinearAlgebra.Array.Util
diff --git a/lib/Numeric/LinearAlgebra/Array.hs b/lib/Numeric/LinearAlgebra/Array.hs
--- a/lib/Numeric/LinearAlgebra/Array.hs
+++ b/lib/Numeric/LinearAlgebra/Array.hs
@@ -28,7 +28,7 @@
 import Numeric.LinearAlgebra.Array.Util
 import Numeric.LinearAlgebra.Array.Internal(namesR)
 import Numeric.LinearAlgebra.Array.Display(printA)
-import Data.Packed(Vector)
+import Numeric.LinearAlgebra.HMatrix(Vector)
 
 -- | Create an 'Array' from a list of parts (@index = 'newIndex' 'None'@).
 index :: Coord t => Name -> [Array t] -> Array t
@@ -51,7 +51,7 @@
     abs _ = error "abs for arrays not defined"
     signum _ = error "signum for arrays not defined"
 
-instance (Coord t, Compat i, Num (NArray i t)) => Fractional (NArray i t) where
+instance (Fractional t, Coord t, Compat i, Num (NArray i t)) => Fractional (NArray i t) where
     fromRational = scalar . fromRational
     (/) = zipArray (/)
     recip = mapArray recip
diff --git a/lib/Numeric/LinearAlgebra/Array/Decomposition.hs b/lib/Numeric/LinearAlgebra/Array/Decomposition.hs
--- a/lib/Numeric/LinearAlgebra/Array/Decomposition.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Decomposition.hs
@@ -24,7 +24,7 @@
 import Numeric.LinearAlgebra.Array.Internal(seqIdx,namesR,sizesR,renameRaw)
 import Numeric.LinearAlgebra.Array.Util
 import Numeric.LinearAlgebra.Array.Solve
-import Numeric.LinearAlgebra hiding (scalar)
+import Numeric.LinearAlgebra.HMatrix hiding (scalar)
 import Data.List
 import System.Random
 --import Control.Parallel.Strategies
@@ -43,7 +43,7 @@
     n = length rs
     dummies = take n $ seqIdx (2*n) "" \\ (namesR t)
     axs = zipWith (\a b->[a,b]) dummies (namesR t)
-    factors = renameRaw core dummies : zipWith renameRaw (map (fromMatrix None None . trans) rs) axs
+    factors = renameRaw core dummies : zipWith renameRaw (map (fromMatrix None None . tr) rs) axs
     core = product $ renameRaw t dummies : zipWith renameRaw (map (fromMatrix None None) rs) axs
 
 {- | Multilinear Singular Value Decomposition (or Tucker's method, see Lathauwer et al.).
@@ -68,14 +68,14 @@
 
 --check trans/ctrans
 usOfSVD m = if rows m < cols m
-        then let (s2,u) = eigSH' $ m <> ctrans m
+        then let (s2,u) = eigSH' $ m <> tr m
                  s = sqrt (abs s2)
               in (u,r s)
-        else let (s2,v) = eigSH' $ ctrans m <> m
+        else let (s2,v) = eigSH' $ tr m <> m
                  s = sqrt (abs s2)
                  u = m <> v <> pinv (diag s)
               in (u,r s)
-    where r s = (ranksv (sqrt eps) (max (rows m) (cols m)) (toList s), s)
+    where r s = (ranksv (sqrt peps) (max (rows m) (cols m)) (toList s), s)
                 -- (rank m, sv m) where sv m = s where (_,s,_) = svd m
 
 
@@ -89,7 +89,7 @@
 
 ------------------------------------------------------------------------
 
-frobT = pnorm PNorm2 . coords
+frobT = norm_2 . coords
 
 ------------------------------------------------------------------------
 
diff --git a/lib/Numeric/LinearAlgebra/Array/Display.hs b/lib/Numeric/LinearAlgebra/Array/Display.hs
--- a/lib/Numeric/LinearAlgebra/Array/Display.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Display.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE FlexibleContexts #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Array.Display
@@ -15,8 +17,7 @@
 ) where
 
 import Numeric.LinearAlgebra.Array.Internal
-import Data.Packed
-import Numeric.Container(format)
+import Numeric.LinearAlgebra.HMatrix
 import Data.List
 import Text.Printf
 
@@ -70,7 +71,7 @@
 formatAux f x = unlines . addds . els . fmt ms $ x where
     fmt [] _ = undefined -- cannot happen
     fmt (g:gs) t
-        | order t == 0 = rect (f (coords t @> 0))
+        | order t == 0 = rect (f (coords t ! 0))
         | order t == 1 =  rect $ unwords $ map f (toList $ coords t)
         | order t == 2 =  decor t $ rect $ w1 $ format " " f (reshape (iDim $ last $ dims t) (coords t))
         | otherwise    = decor t (g ps)
diff --git a/lib/Numeric/LinearAlgebra/Array/Internal.hs b/lib/Numeric/LinearAlgebra/Array/Internal.hs
--- a/lib/Numeric/LinearAlgebra/Array/Internal.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Internal.hs
@@ -52,32 +52,31 @@
     selDims, mapDims,
     takeDiagT, atT,
     firstIdx, fibers, matrixator, matrixatorFree,
-    Coord,
+    Coord,I,
     asMatrix, asVector, asScalar,
     resetCoords,
     debug
 ) where
 
-import Data.Packed
-import Numeric.Container(konst)
-import Data.Complex
+
+import qualified Numeric.LinearAlgebra.HMatrix as LA
+import Numeric.LinearAlgebra.HMatrix hiding (size,scalar,ident)
 import Data.List
-import Numeric.LinearAlgebra((<>),Field,Normed)
-import Control.Applicative
---import Control.Arrow ((***))
 import Data.Function(on)
--- import Foreign(Storable)
--- import Control.Parallel.Strategies
 import Debug.Trace
 
+dim x = LA.size x
+trans x = LA.tr x
+
 ident n = diagRect 0 (konst 1 n) n n
 
 debug m f x = trace (m ++ show (f x)) x
 
 -- | Types that can be elements of the multidimensional arrays.
-class (Num (Vector t), Field t, Normed Vector t, Show t) => Coord t
+class (Num (Vector t), Normed (Vector t), Show t, Numeric t, Indexable (Vector t) t) => Coord t
 instance Coord Double
 instance Coord (Complex Double)
+instance Coord I
 
 -- | indices are denoted by strings, (frequently single-letter)
 type Name = String
@@ -272,9 +271,9 @@
 
 -- | Apply a function (defined on hmatrix 'Vector's) to all elements of a structure.
 -- Use @mapArray (mapVector f)@ for general functions.
-mapArray :: Coord b => (Vector a -> Vector b) -> NArray i a -> NArray i b
+mapArray :: (Coord b) => (Vector a -> Vector b) -> NArray i a -> NArray i b
 mapArray f t
-    | null (dims t) = scalar (f (coords t)@>0)
+    | null (dims t) = scalar (f (coords t)!0)
     | otherwise = mkNArray (dims t) (f (coords t))
 
 liftNA2 f (A d1 v1) (A _d2 v2) = A d1 (f v1 v2)
@@ -393,7 +392,7 @@
 
 -- | Extract the scalar element corresponding to a 0-dimensional array.
 asScalar :: (Coord t) => NArray i t -> t
-asScalar a | order a == 0 = coords a @>0
+asScalar a | order a == 0 = coords a ! 0
            | otherwise = error $ "asScalar requires a 0th order array."
 
 ------------------------------------------------------------------------
@@ -504,3 +503,4 @@
 -- | sequence of n indices with given prefix
 seqIdx :: Int -> String -> [Name]
 seqIdx n prefix = [prefix ++ show k | k <- [1 .. n] ]
+
diff --git a/lib/Numeric/LinearAlgebra/Array/Simple.hs b/lib/Numeric/LinearAlgebra/Array/Simple.hs
--- a/lib/Numeric/LinearAlgebra/Array/Simple.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Simple.hs
@@ -20,7 +20,7 @@
 ) where
 
 import Numeric.LinearAlgebra.Array.Internal
-import Data.Packed
+import Numeric.LinearAlgebra.HMatrix
 import Data.List(intersperse)
 
 
@@ -41,7 +41,7 @@
 type Array t = NArray None t
 
 instance (Coord t) => Show (Array t) where
-    show t | null (dims t) = "scalar "++ show (coords t @>0)
+    show t | null (dims t) = "scalar "++ show (coords t !0)
            | order t == 1 = "index " ++ show n ++" " ++ (show . toList . coords $ t)
            | otherwise = "index "++ show n ++ " [" ++ ps ++ "]"
       where n = head (namesR t)
diff --git a/lib/Numeric/LinearAlgebra/Array/Solve.hs b/lib/Numeric/LinearAlgebra/Array/Solve.hs
--- a/lib/Numeric/LinearAlgebra/Array/Solve.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Solve.hs
@@ -23,14 +23,15 @@
 -- ** Factorized
     solveFactors, solveFactorsH,
 -- * Utilities
-    eps, eqnorm, infoRank,
+    eqnorm, infoRank,
     solve', solveHomog', solveHomog1', solveP'
 ) where
 
 import Numeric.LinearAlgebra.Array.Util
 import Numeric.LinearAlgebra.Exterior
 import Numeric.LinearAlgebra.Array.Internal(mkNArray, selDims, debug, namesR)
-import Numeric.LinearAlgebra hiding (scalar)
+import Numeric.LinearAlgebra.HMatrix hiding (scalar,size)
+--import qualified Numeric.LinearAlgebra.HMatrix as LA
 import Data.List
 import System.Random
 
@@ -38,7 +39,7 @@
 -- | Solution of the linear system a x = b, where a and b are
 -- general multidimensional arrays. The structure and dimension names
 -- of the result are inferred from the arguments.
-solve :: (Compat i, Coord t)
+solve :: (Compat i, Coord t, Field t)
         => NArray i t -- ^ coefficients (a)
         -> NArray i t -- ^ target       (b)
         -> NArray i t -- ^ result       (x)
@@ -59,7 +60,7 @@
 -- general multidimensional array.
 --
 -- If the system is overconstrained we may provide the theoretical rank to get a MSE solution.
-solveHomog :: (Compat i, Coord t)
+solveHomog :: (Compat i, Coord t, Field t)
            =>  NArray i t    -- ^ coefficients (a)
            -> [Name]         -- ^ desired dimensions for the result
                              --   (a subset selected from the target).
@@ -77,7 +78,7 @@
 
 -- | A simpler way to use 'solveHomog', which returns just one solution.
 -- If the system is overconstrained it returns the MSE solution.
-solveHomog1 :: (Compat i, Coord t)
+solveHomog1 :: (Compat i, Coord t, Field t)
             => NArray i t
             -> [Name]
             -> NArray i t
@@ -87,7 +88,7 @@
     where k = product $ map iDim $ selDims (dims m) ns
 
 -- | 'solveHomog1' for single letter index names.
-solveH :: (Compat i, Coord t) => NArray i t -> [Char] -> NArray i t
+solveH :: (Compat i, Coord t, Field t) => NArray i t -> [Char] -> NArray i t
 solveH m ns = solveHomog1 m (map return ns)
 
 
@@ -153,7 +154,7 @@
     g v = v / atT v [n]
     n = size h t - 1
 
-frobT t = realToFrac . pnorm PNorm2 . coords $ t
+frobT t = realToFrac . norm_2 . coords $ t
 --unitT t = t / scalar (frobT t)
 
 dropElemPos k xs = take k xs ++ drop (k+1) xs
@@ -171,7 +172,7 @@
 
 -- | Solution of a multilinear system a x y z ... = b based on alternating least squares.
 mlSolve
-  :: (Compat i, Coord t, Num (NArray i t), Show (NArray i t))
+  :: (Compat i, Coord t, Field t, Num (NArray i t), Show (NArray i t))
      => ALSParam i t     -- ^ optimization parameters
      -> [NArray i t]  -- ^ coefficients (a), given as a list of factors.
      -> [NArray i t]  -- ^ initial solution [x,y,z...]
@@ -191,7 +192,7 @@
 
 -- | Solution of the homogeneous multilinear system a x y z ... = 0 based on alternating least squares.
 mlSolveH
-  :: (Compat i, Coord t, Num (NArray i t), Show (NArray i t))
+  :: (Compat i, Coord t, Field t, Num (NArray i t), Show (NArray i t))
      => ALSParam  i t    -- ^ optimization parameters
      -> [NArray i t]  -- ^ coefficients (a), given as a list of factors.
      -> [NArray i t]  -- ^ initial solution [x,y,z...]
@@ -231,7 +232,7 @@
 {- | Given two arrays a (source) and  b (target), we try to compute linear transformations x,y,z,... for each dimension, such that product [a,x,y,z,...] == b.
 (We can use 'eqnorm' for 'post' processing, or 'id'.)
 -}
-solveFactors :: (Coord t, Random t, Compat i, Num (NArray i t), Show (NArray i t))
+solveFactors :: (Coord t, Field t, Random t, Compat i, Num (NArray i t), Show (NArray i t))
              => Int          -- ^ seed for random initialization
              -> ALSParam i t     -- ^ optimization parameters
              -> [NArray i t] -- ^ source (also factorized)
@@ -263,7 +264,7 @@
 -- [\"pi\",\"qj\", \"rk\", etc.], we try to compute linear transformations
 -- x!\"pi\", y!\"pi\", z!\"rk\", etc. such that product [a,x,y,z,...] == 0.
 solveFactorsH
-  :: (Coord t, Random t, Compat i, Num (NArray i t), Show (NArray i t))
+  :: (Coord t, Random t, Field t, Compat i, Num (NArray i t), Show (NArray i t))
      => Int -- ^ seed for random initialization
      -> ALSParam  i t    -- ^ optimization parameters
      -> [NArray i t] -- ^ coefficient array (a), (also factorized)
diff --git a/lib/Numeric/LinearAlgebra/Array/Util.hs b/lib/Numeric/LinearAlgebra/Array/Util.hs
--- a/lib/Numeric/LinearAlgebra/Array/Util.hs
+++ b/lib/Numeric/LinearAlgebra/Array/Util.hs
@@ -45,7 +45,7 @@
 
 import Numeric.LinearAlgebra.Array.Internal
 import Numeric.LinearAlgebra.Array.Display
-import Data.Packed(Matrix)
+import Numeric.LinearAlgebra.HMatrix(Matrix)
 import Numeric.LinearAlgebra.Array.Simple
 import Data.List(intersperse,sort,foldl1')
 
diff --git a/lib/Numeric/LinearAlgebra/Exterior.hs b/lib/Numeric/LinearAlgebra/Exterior.hs
--- a/lib/Numeric/LinearAlgebra/Exterior.hs
+++ b/lib/Numeric/LinearAlgebra/Exterior.hs
@@ -60,7 +60,7 @@
 -- | The exterior (wedge) product of two tensors. Obtains the union of subspaces.
 --
 --   Implemented as the antisymmetrization of the tensor product.
-(/\) :: (Coord t)
+(/\) :: (Coord t, Fractional t)
      => Tensor t
      -> Tensor t
      -> Tensor t
@@ -98,7 +98,7 @@
                 Nothing -> error $ "dual with different dimensions"
 
 -- | Euclidean inner product of multivectors.
-inner :: (Coord t)
+inner :: (Coord t, Fractional t)
       => Tensor t
       -> Tensor t
       -> Tensor t
diff --git a/lib/Numeric/LinearAlgebra/Multivector.hs b/lib/Numeric/LinearAlgebra/Multivector.hs
--- a/lib/Numeric/LinearAlgebra/Multivector.hs
+++ b/lib/Numeric/LinearAlgebra/Multivector.hs
@@ -23,7 +23,7 @@
     fromTensor
 ) where
 
-import Numeric.LinearAlgebra(toList,reshape,(<\>),(@>))
+import Numeric.LinearAlgebra.HMatrix(toList,reshape,(<\>),atIndex)
 import Numeric.LinearAlgebra.Array.Internal hiding (scalar,coords)
 import Numeric.LinearAlgebra.Array.Display (showBases)
 import Numeric.LinearAlgebra.Tensor hiding (scalar,vector)
@@ -246,7 +246,7 @@
 -- (We do not check that the tensor is actually antisymmetric.)
 fromTensor :: Tensor Double -> Multivector
 fromTensor t = MV $ filter ((/=0.0).fst) $ zip vals basis
-    where vals = map ((@> 0). Array.coords .foldl' partF t) (map (map pred) basis)
+    where vals = map ((`atIndex` 0). Array.coords .foldl' partF t) (map (map pred) basis)
           r = length (dims t)
           n = iDim . head . dims $ t
           partF s i = part s (name,i) where name = iName . head . dims $ s
diff --git a/lib/Numeric/LinearAlgebra/Tensor.hs b/lib/Numeric/LinearAlgebra/Tensor.hs
--- a/lib/Numeric/LinearAlgebra/Tensor.hs
+++ b/lib/Numeric/LinearAlgebra/Tensor.hs
@@ -26,7 +26,7 @@
 ) where
 
 import Numeric.LinearAlgebra.Array.Internal
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.HMatrix hiding (vector)
 import Numeric.LinearAlgebra.Array
 import Data.List(intersperse)
 
@@ -43,7 +43,7 @@
     show (Idx Contra n s) = s ++ "^" ++ show n
 
 instance (Coord t) => Show (Tensor t) where
-    show t | null (dims t) = "scalar "++ show (coords t @>0)
+    show t | null (dims t) = "scalar "++ show (coords t `atIndex` 0)
            | order t == 1 = ixn ++ show n ++" " ++ (show . toList . coords $ t)
            | otherwise = ixn ++ show n ++ " [" ++ ps ++ "]"
       where n = head (namesR t)