diff --git a/.travis.yml b/.travis.yml
--- a/.travis.yml
+++ b/.travis.yml
@@ -7,7 +7,7 @@
   - cabal update
 
   # Try installing some of the build-deps with apt-get for speed.
-  - travis/cabal-apt-install $mode --force-reinstalls --enable-documentation
+  - travis/cabal-apt-install $mode --force-reinstalls
 
 install:
   - cabal configure -flib-Werror $mode
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,7 @@
+0.9
+---
+* `lens 4` support
+
 0.7.0.1
 -------
 * Consolidated the tutorial under on "slug" on the FP Complete site, this broke all of the URLs in the 0.7 documentation.
diff --git a/sparse.cabal b/sparse.cabal
--- a/sparse.cabal
+++ b/sparse.cabal
@@ -1,6 +1,6 @@
 name:          sparse
 category:      Data, Vector
-version:       0.7.0.1
+version:       0.9
 license:       BSD3
 cabal-version: >= 1.8
 license-file:  LICENSE
@@ -62,7 +62,7 @@
     contravariant     >= 0.4.2 && < 1,
     deepseq           >= 1.1   && < 1.4,
     hybrid-vectors    >= 0.1   && < 1,
-    lens              >= 3.9   && < 4,
+    lens              >= 4     && < 5,
     primitive         >= 0.5   && < 0.6,
     transformers      >= 0.3   && < 0.4,
     vector            >= 0.10  && < 0.11,
@@ -72,10 +72,10 @@
 
   exposed-modules:
     Sparse.Matrix
+    Sparse.Matrix.Internal.Array
     Sparse.Matrix.Internal.Fusion
     Sparse.Matrix.Internal.Heap
     Sparse.Matrix.Internal.Key
-    Sparse.Matrix.Internal.Vectored
 
   ghc-options: -Wall
 
diff --git a/src/Sparse/Matrix.hs b/src/Sparse/Matrix.hs
--- a/src/Sparse/Matrix.hs
+++ b/src/Sparse/Matrix.hs
@@ -49,7 +49,7 @@
   , addWith
   , multiplyWith
   -- * Storage
-  , Vectored(..)
+  , Arrayed(..)
   -- * Lenses
   , _Mat, keys, values
   ) where
@@ -73,7 +73,7 @@
 import Prelude hiding (head, last, null)
 import Sparse.Matrix.Internal.Fusion as Fusion
 import Sparse.Matrix.Internal.Key
-import Sparse.Matrix.Internal.Vectored as I
+import Sparse.Matrix.Internal.Array as I
 import Sparse.Matrix.Internal.Heap as Heap hiding (head)
 import Text.Read
 
@@ -82,7 +82,7 @@
 
 -- * Distinguishable Zero
 
-class (Vectored a, Num a) => Eq0 a where
+class (Arrayed a, Num a) => Eq0 a where
   -- | Return whether or not the element is 0.
   --
   -- It may be okay to never return 'True', but you won't be
@@ -153,25 +153,25 @@
 -- * Sparse Matrices
 
 -- invariant: all vectors are the same length
-data Mat a = Mat {-# UNPACK #-} !Int !(U.Vector Word) !(U.Vector Word) !(I.Vector a)
+data Mat a = Mat {-# UNPACK #-} !Int !(U.Vector Word) !(U.Vector Word) !(I.Array a)
  --  deriving (Eq,Ord)
 
-deriving instance (Vectored a, Eq (I.Vector a)) => Eq (Mat a)
+deriving instance (Arrayed a, Eq (I.Array a)) => Eq (Mat a)
 -- Mat n xs ys vs == Mat n' xs' ys' vs' = n == n' && xs == xs' && ys == ys' && vs == vs'
 
-deriving instance (Vectored a, Ord (I.Vector a)) => Ord (Mat a)
+deriving instance (Arrayed a, Ord (I.Array a)) => Ord (Mat a)
 
-instance (Vectored a, Show a) => Show (Mat a) where
+instance (Arrayed a, Show a) => Show (Mat a) where
   showsPrec d m = G.showsPrec d (m^._Mat)
 
-instance (Vectored a, Read a) => Read (Mat a) where
+instance (Arrayed a, Read a) => Read (Mat a) where
   readPrec = (_Mat # ) <$> G.readPrec
 
-instance NFData (I.Vector a) => NFData (Mat a) where
+instance NFData (I.Array a) => NFData (Mat a) where
   rnf (Mat _ xs ys vs) = rnf xs `seq` rnf ys `seq` rnf vs `seq` ()
 
 -- | bundle up the matrix in a form suitable for vector-algorithms
-_Mat :: Vectored a => Iso' (Mat a) (H.Vector U.Vector (Vec a) (Key, a))
+_Mat :: Arrayed a => Iso' (Mat a) (H.Vector U.Vector (Arr a) (Key, a))
 _Mat = iso (\(Mat n xs ys vs) -> H.V (V_Key n xs ys) vs)
            (\(H.V (V_Key n xs ys) vs) -> Mat n xs ys vs)
 {-# INLINE _Mat #-}
@@ -182,7 +182,7 @@
 {-# INLINE keys #-}
 
 -- | Access the keys of a matrix
-values :: Lens (Mat a) (Mat b) (I.Vector a) (I.Vector b)
+values :: Lens (Mat a) (Mat b) (I.Array a) (I.Array b)
 values f (Mat n xs ys vs) = Mat n xs ys <$> f vs
 {-# INLINE values #-}
 
@@ -193,42 +193,39 @@
 eachV :: (Applicative f, G.Vector v a, G.Vector v b) => (a -> f b) -> v a -> f (v b)
 eachV f v = G.fromListN (G.length v) <$> traverse f (G.toList v)
 
-instance (Applicative f, Vectored a, a ~ b) => Each f (Mat a) (Mat b) a b where
-  each f = _Mat $ eachV $ \(k,v) -> (,) k <$> indexed f k v
+instance (Arrayed a, a ~ b) => Each (Mat a) (Mat b) a b where
+  each f = _Mat $ eachV $ \(k,v) -> (,) k <$> f v
   {-# INLINE each #-}
 
-instance (Functor f, Contravariant f, Vectored a) => Contains f (Mat a) where
-  contains = containsIx
-
-instance (Applicative f, Vectored a) => Ixed f (Mat a) where
+instance Arrayed a => Ixed (Mat a) where
   ix ij@(Key i j) f m@(Mat n xs ys vs)
     | Just i' <- xs U.!? l, i == i'
-    , Just j' <- ys U.!? l, j == j' = indexed f ij (vs G.! l) <&> \v -> Mat n xs ys (vs G.// [(l,v)])
+    , Just j' <- ys U.!? l, j == j' = f (vs G.! l) <&> \v -> Mat n xs ys (vs G.// [(l,v)])
     | otherwise = pure m
     where l = search (\k -> Key (xs U.! k) (ys U.! k) >= ij) 0 n
   {-# INLINE ix #-}
 
-instance Vectored a => Vectored (Mat a) where
-  type Vec (Mat a) = V.Vector -- boxed
+instance Arrayed a => Arrayed (Mat a) where
+  type Arr (Mat a) = V.Vector -- boxed
 
-instance (Vectored a, Eq0 a) => Eq0 (Mat a) where
+instance (Arrayed a, Eq0 a) => Eq0 (Mat a) where
   isZero (Mat n _ _ _) = n == 0
   {-# INLINE isZero #-}
 
 -- * Construction
 
 -- | Build a sparse matrix.
-fromList :: Vectored a => [(Key, a)] -> Mat a
+fromList :: Arrayed a => [(Key, a)] -> Mat a
 fromList xs = _Mat # H.modify (Sort.sortBy (compare `on` fst)) (H.fromList xs)
 {-# INLINABLE fromList #-}
 
 -- | Transpose a matrix
-transpose :: Vectored a => Mat a -> Mat a
+transpose :: Arrayed a => Mat a -> Mat a
 transpose xs = xs & _Mat %~ H.modify (Sort.sortBy (compare `on` fst)) . H.map (first swap)
 {-# INLINE transpose #-}
 
 -- | @singleton@ makes a matrix with a singleton value at a given location
-singleton :: Vectored a => Key -> a -> Mat a
+singleton :: Arrayed a => Key -> a -> Mat a
 singleton k v = _Mat # H.singleton (k,v)
 {-# INLINE singleton #-}
 
@@ -236,7 +233,7 @@
 --
 -- >>> ident 4
 -- fromList [(Key 0 0,1),(Key 1 1,1),(Key 2 2,1),(Key 3 3,1)]
-ident :: (Vectored a, Num a) => Int -> Mat a
+ident :: (Arrayed a, Num a) => Int -> Mat a
 ident w = Mat w (U.generate w fromIntegral) (U.generate w fromIntegral) (G.replicate w 1)
 {-# INLINE ident #-}
 
@@ -244,7 +241,7 @@
 --
 -- >>> empty :: Mat Int
 -- fromList []
-empty :: Vectored a => Mat a
+empty :: Arrayed a => Mat a
 empty = Mat 0 U.empty U.empty G.empty
 {-# INLINE empty #-}
 
@@ -265,7 +262,7 @@
 null (Mat n _ _ _) = n == 0
 {-# INLINE null #-}
 
-instance (Vectored a, Eq0 a) => Num (Mat a) where
+instance (Arrayed a, Eq0 a) => Num (Mat a) where
   {-# SPECIALIZE instance Num (Mat Int) #-}
   {-# SPECIALIZE instance Num (Mat Double) #-}
   {-# SPECIALIZE instance Num (Mat (Complex Double)) #-}
@@ -297,7 +294,7 @@
     where m = l + div (h-l) 2
 {-# INLINE search #-}
 
-split1 :: Vectored a => Word -> Word -> Mat a -> (Mat a, Mat a)
+split1 :: Arrayed a => Word -> Word -> Mat a -> (Mat a, Mat a)
 split1 ai bi (Mat n xs ys vs) = (m0,m1)
   where
     !aibi = xor ai bi
@@ -309,7 +306,7 @@
     !m1 = Mat (n-k) xs1 ys1 vs1
 {-# INLINE split1 #-}
 
-split2 :: Vectored a => Word -> Word -> Mat a -> (Mat a, Mat a)
+split2 :: Arrayed a => Word -> Word -> Mat a -> (Mat a, Mat a)
 split2 aj bj (Mat n xs ys vs) = (m0,m1)
   where
     !ajbj = xor aj bj
@@ -323,32 +320,34 @@
 
 -- | Merge two matrices where the indices coincide into a new matrix. This provides for generalized
 -- addition, but where the summation of two non-zero entries is necessarily non-zero.
-addWith :: Vectored a => (a -> a -> a) -> Mat a -> Mat a -> Mat a
+addWith :: Arrayed a => (a -> a -> a) -> Mat a -> Mat a -> Mat a
 addWith f xs ys = _Mat # G.unstream (mergeStreamsWith f (G.stream (xs^._Mat)) (G.stream (ys^._Mat)))
 {-# INLINE addWith #-}
 
 -- | Merge two matrices where the indices coincide into a new matrix. This provides for generalized
 -- addition. Return 'Nothing' for zero.
-addWith0 :: Vectored a => (a -> a -> Maybe a) -> Mat a -> Mat a -> Mat a
+addWith0 :: Arrayed a => (a -> a -> Maybe a) -> Mat a -> Mat a -> Mat a
 addWith0 f xs ys = _Mat # G.unstream (mergeStreamsWith0 f (G.stream (xs^._Mat)) (G.stream (ys^._Mat)))
 {-# INLINE addWith0 #-}
 
 -- | Multiply two matrices using the specified multiplication and addition operation.
-multiplyWith :: Vectored a => (a -> a -> a) -> (Maybe (Heap a) -> Stream (Key, a)) -> Mat a -> Mat a -> Mat a
+multiplyWith :: Arrayed a => (a -> a -> a) -> (Maybe (Heap a) -> Stream (Key, a)) -> Mat a -> Mat a -> Mat a
 {-# INLINEABLE multiplyWith #-}
 multiplyWith times make x0 y0 = case compare (size x0) 1 of
   LT -> empty
-  EQ | size y0 == 1 -> _Mat # (G.unstream $ hint $ make $ go11 (lo x0) (head x0) (lo y0) (head y0))
-     | otherwise    -> _Mat # (G.unstream $ hint $ make $ go12 (lo x0) (head x0) (lo y0) y0 (hi y0))
+  EQ | size y0 == 1 -> hinted $ go11 (lo x0) (head x0) (lo y0) (head y0)
+     | otherwise    -> hinted $ go12 (lo x0) (head x0) (lo y0) y0 (hi y0)
   GT -> case compare (size y0) 1 of
       LT -> empty
-      EQ -> _Mat # (G.unstream $ hint $ make $ go21 (lo x0) x0 (hi x0) (lo y0) (head y0))
-      GT -> _Mat # (G.unstream $ hint $ make $ go22 (lo x0) x0 (hi x0) (lo y0) y0 (hi y0))
+      EQ -> hinted $ go21 (lo x0) x0 (hi x0) (lo y0) (head y0)
+      GT -> hinted $ go22 (lo x0) x0 (hi x0) (lo y0) y0 (hi y0)
   where
-    hint x = sized x $ Max (size x0 * size y0)
+    hinted x = _Mat # G.unstream (sized (make x) (Max (size x0 * size y0)))
+
     go11 (Key i j) a (Key j' k) b
        | j == j' = Just $ Heap.singleton (Key i k) (times a b)
        | otherwise = Nothing
+    {-# INLINE go11 #-}
 
     -- internal cases in go22
     go22L0 xa x ya y yb
@@ -387,24 +386,53 @@
         xiyj = xi .|. yj
         ykxj = yk .|. xj
 
-    go21 _ mx _ yb b = Heap.timesSingleton times (G.stream (mx^._Mat)) yb b -- linear scan. use tree and fast rejects?
-    go12 xa a _ my _ = Heap.singletonTimes times xa a (G.stream (my^._Mat))
+    -- internal cases in go21
+    go21L0 xa x yb b
+      | size x == 1 = go11 xa (head x) yb b
+      | otherwise    = go21 xa x (hi x) yb b
+    {-# INLINE go21L0 #-}
 
+    go21L1 x xb yb b
+      | size x == 1 = go11 xb (head x) yb b
+      | otherwise    = go21 (lo x) x xb yb b
+    {-# INLINE go21L1 #-}
+
+    go21 xa@(Key xai xaj) x xb@(Key xbi xbj) yb@(Key ybj _ybk) b
+      | gts (xor xaj ybj) (xi.|.xj) = Nothing
+      | ges xi xj = case split1 xai xbi x of (m0,m1) -> go21L0 xa m0 yb b `mfby` go21L1 m1 xb yb b -- we can split on i, fby
+      | otherwise = case split2 xaj xbj x of (m0,m1) -> go21L0 xa m0 yb b `madd` go21L1 m1 xb yb b -- we split on j, mix
+      where
+        xi = xor xai xbi
+        xj = xor xaj xbj
+
+    go12R0 xa a ya y
+      | size y == 1 = go11 xa a ya (head y)
+      | otherwise   = go12 xa a ya y (hi y)
+    {-# INLINE go12R0 #-}
+
+    go12R1 xa a y yb
+      | size y == 1 = go11 xa a yb (head y)
+      | otherwise   = go12 xa a (lo y) y yb
+    {-# INLINE go12R1 #-}
+
+    go12 xa@(Key _xai xaj) a ya@(Key yaj yak) y yb@(Key ybj ybk)
+      | gts (xor xaj yaj) (yj.|.yk) = Nothing
+      | ges yj yk = case split1 yaj ybj y of (m0,m1) -> go12R0 xa a ya m0 `madd` go12R1 xa a m1 yb -- we had to split on j, mix
+      | otherwise = case split2 yak ybk y of (m0,m1) -> go12R0 xa a ya m0 `mfby` go12R1 xa a m1 yb -- we can split on k, fby
+      where
+        yj = xor yaj ybj
+        yk = xor yak ybk
+
     madd Nothing xs = xs
     madd xs Nothing = xs
     madd (Just x) (Just y) = Just (mix x y)
-    {-# INLINE madd #-}
 
     mfby Nothing xs = xs
     mfby xs Nothing = xs
     mfby (Just x) (Just y) = Just (fby x y)
-    {-# INLINE mfby #-}
 
     lo (Mat _ xs ys _) = Key (U.head xs) (U.head ys)
-    {-# INLINE lo #-}
 
     hi (Mat _ xs ys _) = Key (U.last xs) (U.last ys)
-    {-# INLINE hi #-}
 
     head (Mat _ _ _ vs) = G.head vs
-    {-# INLINE head #-}
diff --git a/src/Sparse/Matrix/Internal/Array.hs b/src/Sparse/Matrix/Internal/Array.hs
new file mode 100644
--- /dev/null
+++ b/src/Sparse/Matrix/Internal/Array.hs
@@ -0,0 +1,210 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module Sparse.Matrix.Internal.Array
+  ( Arrayed(..)
+  , Array
+  -- * Internals
+  , V_Complex(V_Complex)
+  , MV_Complex(MV_Complex)
+  , V_Pair(V_Pair)
+  , MV_Pair(MV_Pair)
+  ) where
+
+import Control.Monad
+import Data.Complex
+import Data.Int
+import Data.Monoid
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Fusion.Stream as Stream
+import qualified Data.Vector as B
+import Data.Word
+import qualified Sparse.Matrix.Internal.Key as Key
+import Text.Read
+
+-- * Data types that know how to store themselves in a Vector optimally, maximizing the level of unboxing provided.
+
+type Array a = Arr a a
+
+class (G.Vector (Arr a) a, Monoid (Arr a a)) => Arrayed a where
+  type Arr a :: * -> *
+  type Arr a = U.Vector
+
+-- * Unboxed vectors
+
+instance Arrayed ()
+instance Arrayed Double
+instance Arrayed Float
+instance Arrayed Int
+instance Arrayed Int8
+instance Arrayed Int16
+instance Arrayed Int32
+instance Arrayed Int64
+instance Arrayed Key.Key
+instance Arrayed Word
+instance Arrayed Word8
+instance Arrayed Word16
+instance Arrayed Word32
+instance Arrayed Word64
+
+-- * Boxed vectors
+
+instance Arrayed Integer where
+  type Arr Integer = B.Vector
+
+-- * Pairs are boxed or unboxed based on their components
+
+#ifndef HLINT
+data MV_Pair :: * -> * -> * where
+  MV_Pair:: {-# UNPACK #-} !Int -> !(G.Mutable (Arr a) s a) -> !(G.Mutable (Arr b) s b) -> MV_Pair s (a, b)
+
+data V_Pair :: * -> * where
+  V_Pair :: {-# UNPACK #-} !Int -> !(Array a) -> !(Array b) -> V_Pair (a, b)
+#endif
+
+type instance G.Mutable V_Pair = MV_Pair
+
+instance (Arrayed a, Arrayed b) => GM.MVector MV_Pair (a, b) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeReplicate #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  basicLength (MV_Pair l _ _) = l
+  basicUnsafeSlice i n (MV_Pair _ u v)                   = MV_Pair n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v)
+  basicOverlaps (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)      = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2
+  basicUnsafeNew n                                       = liftM2 (MV_Pair n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n)
+  basicUnsafeReplicate n (x, y)                          = liftM2 (MV_Pair n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y)
+  basicUnsafeRead (MV_Pair _ u v) i                      = liftM2 (,) (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i)
+  basicUnsafeWrite (MV_Pair _ u v) i (x, y)              = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y
+  basicClear (MV_Pair _ u v)                             = GM.basicClear u >> GM.basicClear v
+  basicSet (MV_Pair _ u v) (x, y)                        = GM.basicSet u x >> GM.basicSet v y
+  basicUnsafeCopy (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)    = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2
+  basicUnsafeMove (MV_Pair _ u1 v1) (MV_Pair _ u2 v2)    = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2
+  basicUnsafeGrow (MV_Pair _ u v) n                      = liftM2 (MV_Pair n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n)
+
+instance (Arrayed a, Arrayed b) => G.Vector V_Pair (a, b) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq #-}
+  basicLength (V_Pair v _ _) = v
+  basicUnsafeFreeze (MV_Pair n u v)                = liftM2 (V_Pair n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v)
+  basicUnsafeThaw (V_Pair n u v)                   = liftM2 (MV_Pair n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v)
+  basicUnsafeSlice i n (V_Pair _ u v)              = V_Pair n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v)
+  basicUnsafeIndexM (V_Pair _ u v) i               = liftM2 (,) (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i)
+  basicUnsafeCopy (MV_Pair _ mu mv) (V_Pair _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v
+  elemseq _ (x, y) z  = G.elemseq (undefined :: Array a) x
+                       $ G.elemseq (undefined :: Array b) y z
+
+instance (Arrayed a, Arrayed b, Show a, Show b, c ~ (a, b)) => Show (V_Pair c) where
+  showsPrec = G.showsPrec
+
+instance (Arrayed a, Arrayed b, Read a, Read b, c ~ (a, b)) => Read (V_Pair c) where
+  readPrec = G.readPrec
+  readListPrec = readListPrecDefault
+
+instance (Arrayed a, Arrayed b, Eq a, Eq b, c ~ (a, b)) => Eq (V_Pair c) where
+  xs == ys = Stream.eq (G.stream xs) (G.stream ys)
+  {-# INLINE (==) #-}
+
+instance (Arrayed a, Arrayed b, c ~ (a, b)) => Monoid (V_Pair c) where
+  mappend = (G.++)
+  {-# INLINE mappend #-}
+  mempty = G.empty
+  {-# INLINE mempty #-}
+  mconcat = G.concat
+  {-# INLINE mconcat #-}
+
+instance (Arrayed a, Arrayed b) => Arrayed (a, b) where
+  type Arr (a, b) = V_Pair
+
+-- * Complex numbers are boxed or unboxed based on their components
+
+#ifndef HLINT
+data MV_Complex :: * -> * -> * where
+  MV_Complex :: {-# UNPACK #-} !Int -> !(G.Mutable (Arr a) s a) -> !(G.Mutable (Arr a) s a) -> MV_Complex s (Complex a)
+
+data V_Complex :: * -> * where
+  V_Complex :: {-# UNPACK #-} !Int -> !(Array a) -> !(Array a) -> V_Complex (Complex a)
+#endif
+
+type instance G.Mutable V_Complex = MV_Complex
+
+instance (Arrayed a, RealFloat a) => GM.MVector MV_Complex (Complex a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeReplicate #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  basicLength (MV_Complex l _ _) = l
+  basicUnsafeSlice i n (MV_Complex _ u v)                   = MV_Complex n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v)
+  basicOverlaps (MV_Complex _ u1 v1) (MV_Complex _ u2 v2)   = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2
+  basicUnsafeNew n                                          = liftM2 (MV_Complex n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n)
+  basicUnsafeReplicate n (x :+ y)                           = liftM2 (MV_Complex n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y)
+  basicUnsafeRead (MV_Complex _ u v) i                      = liftM2 (:+) (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i)
+  basicUnsafeWrite (MV_Complex _ u v) i (x :+ y)            = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y
+  basicClear (MV_Complex _ u v)                             = GM.basicClear u >> GM.basicClear v
+  basicSet (MV_Complex _ u v) (x :+ y)                      = GM.basicSet u x >> GM.basicSet v y
+  basicUnsafeCopy (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2
+  basicUnsafeMove (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2
+  basicUnsafeGrow (MV_Complex _ u v) n                      = liftM2 (MV_Complex n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n)
+
+instance (Arrayed a, RealFloat a) => G.Vector V_Complex (Complex a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq #-}
+  basicLength (V_Complex v _ _) = v
+  basicUnsafeFreeze (MV_Complex n u v)                   = liftM2 (V_Complex n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v)
+  basicUnsafeThaw (V_Complex n u v)                      = liftM2 (MV_Complex n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v)
+  basicUnsafeSlice i n (V_Complex _ u v)                 = V_Complex n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v)
+  basicUnsafeIndexM (V_Complex _ u v) i                  = liftM2 (:+) (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i)
+  basicUnsafeCopy (MV_Complex _ mu mv) (V_Complex _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v
+  elemseq _ (x :+ y) z = G.elemseq (undefined :: Arr a a) x
+                       $ G.elemseq (undefined :: Arr a a) y z
+
+instance (Arrayed a, RealFloat a, Show a, b ~ Complex a) => Show (V_Complex b) where
+  showsPrec = G.showsPrec
+
+instance (Arrayed a, RealFloat a, Read a, b ~ Complex a) => Read (V_Complex b) where
+  readPrec = G.readPrec
+  readListPrec = readListPrecDefault
+
+instance (Arrayed a, RealFloat a, Eq a, b ~ Complex a) => Eq (V_Complex b) where
+  xs == ys = Stream.eq (G.stream xs) (G.stream ys)
+  {-# INLINE (==) #-}
+
+instance (Arrayed a, RealFloat a, b ~ Complex a) => Monoid (V_Complex b) where
+  mappend = (G.++)
+  {-# INLINE mappend #-}
+  mempty = G.empty
+  {-# INLINE mempty #-}
+  mconcat = G.concat
+  {-# INLINE mconcat #-}
+
+instance (Arrayed a, RealFloat a) => Arrayed (Complex a) where
+  type Arr (Complex a) = V_Complex
diff --git a/src/Sparse/Matrix/Internal/Heap.hs b/src/Sparse/Matrix/Internal/Heap.hs
--- a/src/Sparse/Matrix/Internal/Heap.hs
+++ b/src/Sparse/Matrix/Internal/Heap.hs
@@ -26,8 +26,6 @@
   , fromAscList
   , streamHeapWith
   , streamHeapWith0
-  , timesSingleton
-  , singletonTimes
   ) where
 
 import Control.Applicative
@@ -36,7 +34,6 @@
 import Data.Monoid
 import Data.Vector.Fusion.Stream.Monadic hiding (singleton, fromList, head, tail)
 import Data.Vector.Fusion.Stream.Size
-import Data.Vector.Fusion.Util
 import Sparse.Matrix.Internal.Key
 import Prelude hiding (head, tail)
 
@@ -97,14 +94,20 @@
 
 pops :: [Heap a] -> [Heap a] -> [Heap a] -> [Heap a]
 pops xs     []     [] = xs
-pops (x:xs) ls     rs = [fbys (Prelude.foldl mix x xs) ls rs]
+pops (x:xs) ls     rs = [fbys (merge x xs) ls rs]
 pops []     (l:ls) rs = [fbys l ls rs]
 pops []     []     rs = case reverse rs of
   f:fs -> [fbys f fs []]
   _    -> [] -- caught above by the 'go as [] []' case
 
+merge :: Heap a -> [Heap a] -> Heap a
+merge x (y:ys) = case ys of
+  (z:zs) -> mix x y `mix` merge z zs
+  []     -> mix x y
+merge x [] = x
+
 pop :: [Heap a] -> [Heap a] -> [Heap a] -> Maybe (Heap a)
-pop (x:xs) ls     rs = Just $ fbys (Prelude.foldl mix x xs) ls rs
+pop (x:xs) ls     rs = Just $ fbys (merge x xs) ls rs
 pop []     (l:ls) rs = Just $ fbys l ls rs
 pop []     []     rs = case reverse rs of
   f:fs -> Just (fbys f fs [])
@@ -175,37 +178,3 @@
   step Finished = return Done
   {-# INLINE [1] step #-}
 {-# INLINE [0] streamHeapWith0 #-}
-
--- | This is an internal 'Heap' fusion combinator used to multiply on the right by a singleton 'Key'/value pair.
-timesSingleton :: (a -> b -> c) -> Stream Id (Key, a) -> Key -> b -> Maybe (Heap c)
-timesSingleton f (Stream stepa sa0 _) (Key j k) b = start sa0 where
-  start sa = case unId (stepa sa) of
-    Yield (Key i j', a) sa'
-      | j == j'         -> Just $ run (singleton (Key i k) (f a b)) sa'
-      | otherwise       -> start sa'
-    Skip sa' -> start sa'
-    Done     -> Nothing
-  run h sa = case unId (stepa sa) of
-    Yield (Key i j', a) sa'
-      | j == j'   -> run (h `mix` singleton (Key i k) (f a b)) sa'
-      | otherwise -> run h sa'
-    Skip sa' -> run h sa'
-    Done     -> h
-{-# INLINE timesSingleton #-}
-
--- | This is an internal 'Heap' fusion combinator used to multiply on the right by a singleton 'Key'/value pair.
-singletonTimes :: (a -> b -> c) -> Key -> a -> Stream Id (Key, b) -> Maybe (Heap c)
-singletonTimes f (Key i j) a (Stream stepb sb0 _) = start sb0 where
-  start sb = case unId (stepb sb) of
-    Yield (Key j' k, b) sb'
-      | j == j'   -> Just $ run (singleton (Key i k) (f a b)) sb'
-      | otherwise -> start sb'
-    Skip sb' -> start sb'
-    Done     -> Nothing
-  run h sb = case unId (stepb sb) of
-    Yield (Key j' k, b) sb'
-      | j == j'   -> run (h `mix` singleton (Key i k) (f a b)) sb'
-      | otherwise -> run h sb'
-    Skip sb' -> run h sb'
-    Done     -> h
-{-# INLINE singletonTimes #-}
diff --git a/src/Sparse/Matrix/Internal/Vectored.hs b/src/Sparse/Matrix/Internal/Vectored.hs
deleted file mode 100644
--- a/src/Sparse/Matrix/Internal/Vectored.hs
+++ /dev/null
@@ -1,133 +0,0 @@
-{-# LANGUAGE CPP #-}
-{-# LANGUAGE GADTs #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE KindSignatures #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-module Sparse.Matrix.Internal.Vectored
-  ( Vectored(..)
-  , Vector
-  -- * Internals
-  , V_Complex(V_Complex)
-  , MV_Complex(MV_Complex)
-  ) where
-
-import Control.Monad
-import Data.Complex
-import Data.Int
-import Data.Monoid
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Generic.Mutable as GM
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Fusion.Stream as Stream
-import qualified Data.Vector as B
-import Data.Word
-import qualified Sparse.Matrix.Internal.Key as Key
-import Text.Read
-
--- * Data types that know how to store themselves in a Vector optimally, maximizing the level of unboxing provided.
-
-type Vector a = Vec a a
-
-class (G.Vector (Vec a) a, Monoid (Vec a a)) => Vectored a where
-  type Vec a :: * -> *
-  type Vec a = U.Vector
-
--- * Unboxed vectors
-
-instance Vectored ()
-instance Vectored Double
-instance Vectored Float
-instance Vectored Int
-instance Vectored Int8
-instance Vectored Int16
-instance Vectored Int32
-instance Vectored Int64
-instance Vectored Key.Key
-instance Vectored Word
-instance Vectored Word8
-instance Vectored Word16
-instance Vectored Word32
-instance Vectored Word64
-
--- * Boxed vectors
-
-instance Vectored Integer where
-  type Vec Integer = B.Vector
-
--- * Complex numbers are boxed or unboxed based on their components
-
-#ifndef HLINT
-data MV_Complex :: * -> * -> * where
-  MV_Complex :: {-# UNPACK #-} !Int -> !(G.Mutable (Vec a) s a) -> !(G.Mutable (Vec a) s a) -> MV_Complex s (Complex a)
-
-data V_Complex :: * -> * where
-  V_Complex :: {-# UNPACK #-} !Int -> !(Vector a) -> !(Vector a) -> V_Complex (Complex a)
-#endif
-
-type instance G.Mutable V_Complex = MV_Complex
-
-instance (Vectored a, RealFloat a) => GM.MVector MV_Complex (Complex a) where
-  {-# INLINE basicLength #-}
-  {-# INLINE basicUnsafeSlice #-}
-  {-# INLINE basicOverlaps #-}
-  {-# INLINE basicUnsafeNew #-}
-  {-# INLINE basicUnsafeReplicate #-}
-  {-# INLINE basicUnsafeRead #-}
-  {-# INLINE basicUnsafeWrite #-}
-  {-# INLINE basicClear #-}
-  {-# INLINE basicSet #-}
-  {-# INLINE basicUnsafeCopy #-}
-  {-# INLINE basicUnsafeGrow #-}
-  basicLength (MV_Complex l _ _) = l
-  basicUnsafeSlice i n (MV_Complex _ u v)                   = MV_Complex n (GM.basicUnsafeSlice i n u) (GM.basicUnsafeSlice i n v)
-  basicOverlaps (MV_Complex _ u1 v1) (MV_Complex _ u2 v2)   = GM.basicOverlaps u1 u2 || GM.basicOverlaps v1 v2
-  basicUnsafeNew n                                          = liftM2 (MV_Complex n) (GM.basicUnsafeNew n) (GM.basicUnsafeNew n)
-  basicUnsafeReplicate n (x :+ y)                           = liftM2 (MV_Complex n) (GM.basicUnsafeReplicate n x) (GM.basicUnsafeReplicate n y)
-  basicUnsafeRead (MV_Complex _ u v) i                      = liftM2 (:+) (GM.basicUnsafeRead u i) (GM.basicUnsafeRead v i)
-  basicUnsafeWrite (MV_Complex _ u v) i (x :+ y)            = GM.basicUnsafeWrite u i x >> GM.basicUnsafeWrite v i y
-  basicClear (MV_Complex _ u v)                             = GM.basicClear u >> GM.basicClear v
-  basicSet (MV_Complex _ u v) (x :+ y)                      = GM.basicSet u x >> GM.basicSet v y
-  basicUnsafeCopy (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeCopy u1 u2 >> GM.basicUnsafeCopy v1 v2
-  basicUnsafeMove (MV_Complex _ u1 v1) (MV_Complex _ u2 v2) = GM.basicUnsafeMove u1 u2 >> GM.basicUnsafeMove v1 v2
-  basicUnsafeGrow (MV_Complex _ u v) n                      = liftM2 (MV_Complex n) (GM.basicUnsafeGrow u n) (GM.basicUnsafeGrow v n)
-
-instance (Vectored a, RealFloat a) => G.Vector V_Complex (Complex a) where
-  {-# INLINE basicLength #-}
-  {-# INLINE basicUnsafeFreeze #-}
-  {-# INLINE basicUnsafeThaw #-}
-  {-# INLINE basicUnsafeSlice #-}
-  {-# INLINE basicUnsafeIndexM #-}
-  {-# INLINE elemseq #-}
-  basicLength (V_Complex v _ _) = v
-  basicUnsafeFreeze (MV_Complex n u v)                   = liftM2 (V_Complex n) (G.basicUnsafeFreeze u) (G.basicUnsafeFreeze v)
-  basicUnsafeThaw (V_Complex n u v)                      = liftM2 (MV_Complex n) (G.basicUnsafeThaw u) (G.basicUnsafeThaw v)
-  basicUnsafeSlice i n (V_Complex _ u v)                 = V_Complex n (G.basicUnsafeSlice i n u) (G.basicUnsafeSlice i n v)
-  basicUnsafeIndexM (V_Complex _ u v) i                  = liftM2 (:+) (G.basicUnsafeIndexM u i) (G.basicUnsafeIndexM v i)
-  basicUnsafeCopy (MV_Complex _ mu mv) (V_Complex _ u v) = G.basicUnsafeCopy mu u >> G.basicUnsafeCopy mv v
-  elemseq _ (x :+ y) z = G.elemseq (undefined :: Vec a a) x
-                       $ G.elemseq (undefined :: Vec a a) y z
-
-instance (Vectored a, RealFloat a, Show a, b ~ Complex a) => Show (V_Complex b) where
-  showsPrec = G.showsPrec
-
-instance (Vectored a, RealFloat a, Read a, b ~ Complex a) => Read (V_Complex b) where
-  readPrec = G.readPrec
-  readListPrec = readListPrecDefault
-
-instance (Vectored a, RealFloat a, Eq a, b ~ Complex a) => Eq (V_Complex b) where
-  xs == ys = Stream.eq (G.stream xs) (G.stream ys)
-  {-# INLINE (==) #-}
-
-instance (Vectored a, RealFloat a, b ~ Complex a) => Monoid (V_Complex b) where
-  mappend = (G.++)
-  {-# INLINE mappend #-}
-  mempty = G.empty
-  {-# INLINE mempty #-}
-  mconcat = G.concat
-  {-# INLINE mconcat #-}
-
-instance (Vectored a, RealFloat a) => Vectored (Complex a) where
-  type Vec (Complex a) = V_Complex