diff --git a/accelerate-arithmetic.cabal b/accelerate-arithmetic.cabal
--- a/accelerate-arithmetic.cabal
+++ b/accelerate-arithmetic.cabal
@@ -1,10 +1,10 @@
 Name:             accelerate-arithmetic
-Version:          0.0.1
+Version:          0.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
 Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
-Homepage:         http://code.haskell.org/~thielema/accelerate-arithmetic/
+Homepage:         http://hub.darcs.net/thielema/accelerate-arithmetic/
 Category:         Math
 Synopsis:         Linear algebra and interpolation using the Accelerate framework
 Description:
@@ -17,66 +17,48 @@
 Build-Type:       Simple
 
 Source-Repository this
-  Tag:         0.0.1
+  Tag:         0.1
   Type:        darcs
-  Location:    http://code.haskell.org/~thielema/accelerate-arithmetic/
+  Location:    http://hub.darcs.net/thielema/accelerate-arithmetic/
 
 Source-Repository head
   Type:        darcs
-  Location:    http://code.haskell.org/~thielema/accelerate-arithmetic/
+  Location:    http://hub.darcs.net/thielema/accelerate-arithmetic/
 
 Library
   Build-Depends:
     accelerate-utility >=0.1 && <0.2,
     accelerate >=0.15 && <0.16,
     utility-ht >=0.0.8 && <0.1,
-    QuickCheck >=2.4 && <2.8,
-    base >=4.5 && <4.8
+    QuickCheck >=2.4 && <3,
+    base >=4.5 && <4.10
 
   GHC-Options:      -Wall -fwarn-missing-import-lists
-  Hs-Source-Dirs:   src
+  Hs-Source-Dirs:   src, private
   Default-Language: Haskell98
   Exposed-Modules:
-    Data.Array.Accelerate.Arithmetic.LinearAlgebra
-    Data.Array.Accelerate.Arithmetic.Sparse
     Data.Array.Accelerate.Arithmetic.Interpolation
+    Data.Array.Accelerate.LinearAlgebra
+    Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse
+    Data.Array.Accelerate.LinearAlgebra.Matrix.Banded
   Other-Modules:
     Data.Array.Accelerate.Arithmetic.Example
+    Data.Array.Accelerate.LinearAlgebra.Private
 
 Test-Suite test
   Type: exitcode-stdio-1.0
   Main-Is:          Test.hs
   GHC-Options:      -Wall -fwarn-missing-import-lists
-  Hs-Source-Dirs:   test
+  Hs-Source-Dirs:   test, private
   Default-Language: Haskell98
   Build-Depends:
     accelerate-arithmetic,
+    accelerate-utility,
     accelerate,
     QuickCheck,
     base
   Other-Modules:
+    Data.Array.Accelerate.LinearAlgebra.Private
     Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra
     Test.Data.Array.Accelerate.Arithmetic.Sparse
     Test.Data.Array.Accelerate.Arithmetic.Utility
-
-Benchmark newton-inverse
-  Type:             exitcode-stdio-1.0
-  Main-Is:          NewtonInverse.hs
-  Hs-Source-Dirs:   benchmark
-  Other-Modules:    CUBLASBatched
-  Default-Language: Haskell98
-  GHC-Options:      -Wall -threaded
-  GHC-Prof-Options: -fprof-auto
-  Build-Depends:
-    accelerate-arithmetic,
-    accelerate-utility,
-    accelerate-cuda >=0.15 && <0.16,
-    cublas >=0.2.0.2 && <0.3,
-    cuda >=0.5 && <0.7,
-    accelerate,
-    pooled-io >=0.0 && <0.1,
-    timeit >=1.0 && <1.1,
-    hmatrix >=0.15.2 && <0.16,
-    random >=1.0.1 && <1.1,
-    utility-ht,
-    base
diff --git a/benchmark/CUBLASBatched.hs b/benchmark/CUBLASBatched.hs
deleted file mode 100644
--- a/benchmark/CUBLASBatched.hs
+++ /dev/null
@@ -1,269 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ConstraintKinds #-}
-{-# LANGUAGE FlexibleContexts #-}
-module CUBLASBatched where
-
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg
-
-import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc
-import Data.Array.Accelerate.Utility.Lift.Acc (acc, expr)
-
-import Data.Array.Accelerate.Array.Sugar (EltRepr)
-import Data.Array.Accelerate (Array, DIM3, Acc, Z (..), (:.) (..), Exp)
-import qualified Data.Array.Accelerate.CUDA.Foreign as AF
-import qualified Data.Array.Accelerate.CUDA as AC
-import qualified Data.Array.Accelerate as A
-
-import qualified Foreign.CUDA.Cublas as Cublas
-import Foreign.CUDA.Ptr (DevicePtr, castDevPtr, advanceDevPtr)
-
-import Foreign.C.Types (CFloat, CDouble)
-import Foreign.Storable (Storable)
-
-import Data.Tuple.HT (uncurry3)
-
-
-type Matrix ix = Array (ix :. Int :. Int)
-type Vector ix = Array (ix :. Int)
-type Scalar ix = Array ix
-
-mul ::
-   (A.Shape ix, A.Slice ix, Eq ix, Element a, A.Elt a, A.IsNum a) =>
-   Cublas.Handle ->
-   Exp a ->
-   ALinAlg.Matrix ix a -> ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a
-mul handle alpha a b =
-   A.foreignAcc
-      (AF.CUDAForeignAcc "mul" $ uncurry3 $ mulPlain handle)
-      (Acc.modify (expr,acc,acc) $ \(alpha0, a0, b0) ->
-         A.map (alpha0 *) $
-         ALinAlg.multiplyMatrixMatrix a0 b0)
-   $
-   A.lift (A.unit alpha, a, b)
-
-mulPlain ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   A.Scalar a -> Matrix ix a -> Matrix ix a ->
-   AF.CIO (Matrix ix a)
-mulPlain handle alpha a b = do
-   let (aNumMatrices :. n  :. k) = A.arrayShape a
-   let (bNumMatrices :. _k :. m) = A.arrayShape b
-   let numMatrices =
-          if aNumMatrices == bNumMatrices
-            then aNumMatrices
-            else error "mul: mismatching shapes of matrix arrays"
-   c <- AF.allocateArray (numMatrices :. n :. m)
-   (pas, lda) <- arrayPtrs a
-   (pbs, ldb) <- arrayPtrs b
-   (pcs, ldc) <- arrayPtrs c
-   AF.liftIO $
-      Cublas.gemmBatched handle Cublas.N Cublas.N m n k
-         (storableFromScalar alpha)
-         pbs ldb
-         pas lda
-         0
-         pcs ldc
-         (A.arraySize numMatrices)
-   return c
-
-mac ::
-   (A.Shape ix, A.Slice ix, Eq ix, Element a, A.Elt a, A.IsNum a) =>
-   Cublas.Handle ->
-   Exp a -> ALinAlg.Matrix ix a -> ALinAlg.Matrix ix a ->
-   Exp a -> ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a
-mac handle alpha a b beta c =
-   A.foreignAcc
-      (AF.CUDAForeignAcc "mac" $
-       \((alpha0, a0, b0), (beta0, c0)) ->
-          macPlain handle alpha0 a0 b0 beta0 c0)
-      (Acc.modify ((expr,acc,acc),(expr,acc)) $
-       \((alpha0, a0, b0), (beta0, c0)) ->
-         A.zipWith (+)
-            (A.map (alpha0 *) $
-             ALinAlg.multiplyMatrixMatrix a0 b0)
-            (A.map (beta0 *) c0))
-   $
-   A.lift ((A.unit alpha, a, b), (A.unit beta, c))
-
-macPlain ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   A.Scalar a -> Matrix ix a -> Matrix ix a ->
-   A.Scalar a -> Matrix ix a ->
-   AF.CIO (Matrix ix a)
-macPlain handle alpha a b beta c = do
-   let (aNumMatrices :. n  :. k ) = A.arrayShape a
-   let (bNumMatrices :. _k :. m ) = A.arrayShape b
-   let (cNumMatrices :. n' :. m') = A.arrayShape c
-   let numMatrices =
-          if aNumMatrices == bNumMatrices
-             &&
-             aNumMatrices == cNumMatrices
-            then aNumMatrices
-            else error "mac: mismatching shapes of matrix arrays"
-   d <- AF.allocateArray (numMatrices :. n' :. m')
-   AF.copyArray c d
-   (pas, lda) <- arrayPtrs a
-   (pbs, ldb) <- arrayPtrs b
-   (pds, ldd) <- arrayPtrs d
-   AF.liftIO $
-      Cublas.gemmBatched handle Cublas.N Cublas.N m n k
-         (storableFromScalar alpha)
-         pbs ldb
-         pas lda
-         (storableFromScalar beta)
-         pds ldd
-         (A.arraySize numMatrices)
-   return d
-
-lu ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   ALinAlg.Matrix ix a ->
-   (ALinAlg.Matrix ix a, ALinAlg.Vector ix Int, ALinAlg.Scalar ix Int)
-lu handle =
-   A.unlift
-   .
-   A.foreignAcc
-      (AF.CUDAForeignAcc "lu" $ luPlain handle)
-      (error "Requires CUDA backend")
-
-luPlain ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   Matrix ix a ->
-   AF.CIO (Matrix ix a, Vector ix Int, Scalar ix Int)
-luPlain handle a = do
-   let sh@(numMatrices :. n  :. k) = A.arrayShape a
-   let size =
-          if n == k
-            then n
-            else error "lu: matrices must have square shape"
-   b <- AF.allocateArray sh
-   AF.copyArray a b
-   (pbs, ldb) <- arrayPtrs b
-
-   pivot <- AF.allocateArray (numMatrices :. size)
-   pivotPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray pivot
-
-   info <- AF.allocateArray numMatrices
-   infoPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray info
-
-   AF.liftIO $
-      Cublas.getrfBatched handle size
-         pbs ldb
-         pivotPtr infoPtr
-         (A.arraySize numMatrices)
-   return (b, pivot, info)
-
-
-luInv ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   (ALinAlg.Matrix ix a, ALinAlg.Vector ix Int, ALinAlg.Scalar ix Int) ->
-   ALinAlg.Matrix ix a
-luInv handle =
-   A.foreignAcc
-      (AF.CUDAForeignAcc "luInv" $ luInvPlain handle)
-      (error "Requires CUDA backend")
-   .
-   A.lift
-
-luInvPlain ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   (Matrix ix a, Vector ix Int, Scalar ix Int) ->
-   AF.CIO (Matrix ix a)
-luInvPlain handle (a, pivot, info) = do
-   let sh@(numMatrices :. n  :. k) = A.arrayShape a
-   let size =
-          if n == k
-            then n
-            else error "luInv: matrices must have square shape"
-   c <- AF.allocateArray sh
-   AF.copyArray a c
-   (pas, lda) <- arrayPtrs a
-   (pcs, ldc) <- arrayPtrs c
-
-   pivotPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray pivot
-   infoPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray info
-
-   AF.liftIO $
-      Cublas.getriBatched handle size
-         pas lda
-         pivotPtr
-         pcs ldc
-         infoPtr
-         (A.arraySize numMatrices)
-   return c
-
-
-inv ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
-   Cublas.Handle ->
-   ALinAlg.Matrix ix a ->
-   (ALinAlg.Matrix ix a, ALinAlg.Scalar ix Int)
-inv handle a =
-   let sol@(_,_,info) = lu handle a
-   in  (luInv handle sol, info)
-
-
-type Element a =
-        (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),
-         Fractional (StorableOf a),
-         Cublas.Cublas (StorableOf a),
-         Storable (StorableOf a),
-         Real a)
-
-type family StorableOf float
-type instance StorableOf Float = CFloat
-type instance StorableOf Double = CDouble
-
-storableFromScalar ::
-   (Real a, StorableOf a ~ b, Fractional b) => A.Scalar a -> b
-storableFromScalar x = realToFrac $ A.indexArray x Z
-
-arrayPtrs ::
-   (Storable a, StorableOf e ~ a,
-    A.Shape ix,
-    AF.DevicePtrs (EltRepr e) ~ ((), DevicePtr e)) =>
-   Array (ix :. Int :. Int) e -> AF.CIO ([DevicePtr a], Int)
-arrayPtrs arr = do
-   let (numMatrices :. n  :. k) = A.arrayShape arr
-   pa <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray arr
-   return (genPointers (n*k) pa (A.arraySize numMatrices), k)
-
-genPointers ::
-   (Storable a) =>
-   Int -> DevicePtr a -> Int -> [DevicePtr a]
-genPointers size p n =
-   take n $ iterate (flip advanceDevPtr size) p
-
-
-genMatrices :: (Acc (Array DIM3 Double), Acc (Array DIM3 Double))
-genMatrices = (a,b)
-   where
-   a = A.generate (A.constant sha) $ \ix ->
-      let (Z :. i :. j :. k) = unlift ix
-      in A.fromIntegral (i+j+k)
-   b = A.generate (A.constant shb) $ \ix ->
-      let (Z :. i :. j :. k) = unlift ix
-      in A.fromIntegral (i+j+k)
-   numMats = 100 :: Int
-   sha = Z :. numMats :. (3 :: Int) :. (4 :: Int)
-   shb = Z :. numMats :. (4 :: Int) :. (2 :: Int)
-   unlift :: Exp (Z :. Int :. Int :. Int)
-      -> Z :. Exp Int :. Exp Int :. Exp Int
-   unlift = A.unlift
-
-test :: IO ()
-test = do
-   handle <- Cublas.create
-   print genMatrices
-   print $ AC.run $
-      case genMatrices of
-         (a,b) -> mul handle 1 a b
diff --git a/benchmark/NewtonInverse.hs b/benchmark/NewtonInverse.hs
deleted file mode 100644
--- a/benchmark/NewtonInverse.hs
+++ /dev/null
@@ -1,224 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ConstraintKinds #-}
-{-# LANGUAGE FlexibleContexts #-}
-module Main where
-
-import qualified CUBLASBatched as Batched
-import qualified Foreign.CUDA.Cublas as Cublas
-
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg
-import qualified Data.Array.Accelerate.Utility.Loop as Loop
-import qualified Data.Array.Accelerate.CUDA as CUDA
-import qualified Data.Array.Accelerate as A
-import Data.Array.Accelerate (All(All), Z(Z), (:.)((:.)))
-
-import qualified Control.Concurrent.PooledIO.Independent as Pooled
-
-import qualified Data.Packed.Matrix as Matrix
-import qualified Data.Packed.Vector as Vector
-import qualified Numeric.Container as Container
-import qualified Numeric.LinearAlgebra.Algorithms as HMLinAlg
-
-import Numeric.Container (Container, (<>))
-import Data.Packed.Matrix (Matrix)
-import Data.Packed.Vector (Vector)
-
-import qualified System.Random as Rnd
-import System.TimeIt (timeIt)
-
-import Text.Printf (printf)
-
-import qualified Data.List.HT as ListHT
-import Data.Function.HT (nest)
-import Data.Tuple.HT (mapPair)
-
-
-
-newtonInverseStep ::
-   (Num a, Container Vector a, Container.Product a) =>
-   Matrix a -> Matrix a -> Matrix a
-newtonInverseStep a x =
-   Container.sub
-      (Container.scale 2 x)
-      (x <> a <> x)
-
-newtonInverse ::
-   (Num a, Container Vector a, Container.Product a) =>
-   Int -> Matrix a -> Matrix a -> Matrix a
-newtonInverse count start a =
-   nest count (newtonInverseStep a) start
-
-
-newtonInverseCUBLASStep, newtonInverseCUBLASStepMul ::
-   (A.Shape ix, A.Slice ix, Eq ix, Batched.Element a, A.IsNum a, A.Elt a) =>
-   Cublas.Handle ->
-   ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a
-newtonInverseCUBLASStep h a x =
-   Batched.mac h (-1) x (Batched.mul h 1 a x) 2 x
-
-newtonInverseCUBLASStepMul h a x =
-   A.zipWith (-) (A.map (2*) x) $
-   Batched.mul h 1 x $ Batched.mul h 1 a x
-
-newtonInverseCUBLAS ::
-   (A.Shape ix, A.Slice ix, Eq ix, Batched.Element a, A.IsNum a, A.Elt a) =>
-   Cublas.Handle ->
-   A.Exp Int ->
-   ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a ->
-   ALinAlg.Matrix ix a
-newtonInverseCUBLAS h n seed a =
-   Loop.nest n (newtonInverseCUBLASStep h a) seed
-
-
-randomMatrixInv :: Int -> (Matrix Double, Matrix Double)
-randomMatrixInv size =
-   let x =
-          Matrix.fromLists $ take size $ ListHT.sliceVertical size $
-          Rnd.randomRs (-1,1::Double) $ Rnd.mkStdGen 42
-   in  (x, HMLinAlg.inv x)
-
-
-parallel :: [a] -> (Int -> a -> IO ()) -> IO ()
-parallel xs f = Pooled.run $ zipWith f [0 ..] xs
-
-disturbedMatrices ::
-   (Container Vector a) =>
-   Matrix a -> [a] -> [Matrix a]
-disturbedMatrices x yelems =
-   let size = Matrix.rows x
-   in  map (Container.add x . Matrix.fromLists) $
-       ListHT.sliceVertical size $
-       ListHT.sliceVertical size $
-       yelems
-
-mainHMatrixDirect ::
-   (Show a, Container Vector a, HMLinAlg.Field a) =>
-   String ->
-   Int -> (Matrix a, Matrix a) -> [a] -> IO ()
-mainHMatrixDirect typ numberOfMatrices (x, _xinv) yelems = do
-   let yinvs = map HMLinAlg.inv $ disturbedMatrices x yelems
-   putStrLn $ "hmatrix-direct-" ++ typ
-   timeIt $ parallel (take numberOfMatrices yinvs) $ \ n y ->
-      writeFile (printf "/tmp/hmatrix-direct-%s%03d.txt" typ n) $ show y
-
-mainHMatrix ::
-   (Show a, Container Vector a, Container.Product a) =>
-   String ->
-   Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()
-mainHMatrix typ numberOfMatrices newtonIts (x, xinv) yelems = do
-   let yinvs = map (newtonInverse newtonIts xinv) $ disturbedMatrices x yelems
-   putStrLn $ "hmatrix-" ++ typ
-   timeIt $ parallel (take numberOfMatrices yinvs) $ \ n y ->
-      writeFile (printf "/tmp/hmatrix-%s%03d.txt" typ n) $ show y
-
-
-mainCUDA ::
-   (A.Elt a, A.IsNum a, Container.Element a) =>
-   String ->
-   Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()
-mainCUDA typ numberOfMatrices newtonIts (xm, xinvm) yelems = do
-   let size = Matrix.rows xm
-       matrixAccFromHM =
-          A.fromList (Z :. size :. size) .
-          Vector.toList . Matrix.flatten
-       xarr = matrixAccFromHM xm
-       xinvarr = matrixAccFromHM xinvm
-
-   let ysarr =
-          A.fromList (Z :. numberOfMatrices :. size :. size) yelems
-       rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)
-       yinvs =
-          CUDA.run1
-             (\args ->
-                case A.unlift args of
-                   (x, xinv, ys) ->
-                      ALinAlg.newtonInverse (A.constant newtonIts) (rep xinv) $
-                      A.zipWith (+) ys (rep x))
-             (xarr, xinvarr, ysarr)
-
-   putStrLn $ "cuda-" ++ typ
-   timeIt $ writeFile ("/tmp/cuda-"++typ++".txt") $ show yinvs
-
-
-mainCUBLASDirect ::
-   (Batched.Element a, Container.Element a, A.IsNum a, A.Elt a) =>
-   String ->
-   Int -> (Matrix a, Matrix a) -> [a] -> IO ()
-mainCUBLASDirect typ numberOfMatrices (xm, _xinvm) yelems = do
-   let size = Matrix.rows xm
-       matrixAccFromHM =
-          A.fromList (Z :. size :. size) .
-          Vector.toList . Matrix.flatten
-       xarr = matrixAccFromHM xm
-
-   handle <- Cublas.create
-   let ysarr =
-          A.fromList (Z :. numberOfMatrices :. size :. size) yelems
-       rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)
-       yinvs =
-          CUDA.run1
-             (\args ->
-                case A.unlift args of
-                   (x, ys) ->
-                      fst $ Batched.inv handle $ A.zipWith (+) ys (rep x))
-             (xarr, ysarr)
-
-   putStrLn $ "cublas-direct-" ++ typ
-   timeIt $ writeFile ("/tmp/cublas-direct-"++typ++".txt") $ show yinvs
-
-
-mainCUBLAS ::
-   (Batched.Element a, Container.Element a, A.IsNum a, A.Elt a) =>
-   String ->
-   Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()
-mainCUBLAS typ numberOfMatrices newtonIts (xm, xinvm) yelems = do
-   let size = Matrix.rows xm
-       matrixAccFromHM =
-          A.fromList (Z :. size :. size) .
-          Vector.toList . Matrix.flatten
-       xarr = matrixAccFromHM xm
-       xinvarr = matrixAccFromHM xinvm
-
-   handle <- Cublas.create
-   let ysarr =
-          A.fromList (Z :. numberOfMatrices :. size :. size) yelems
-       rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)
-       yinvs =
-          CUDA.run1
-             (\args ->
-                case A.unlift args of
-                   (x, xinv, ys) ->
-                      newtonInverseCUBLAS handle (A.constant newtonIts) (rep xinv) $
-                      A.zipWith (+) ys (rep x))
-             (xarr, xinvarr, ysarr)
-
-   putStrLn $ "cublas-" ++ typ
-   timeIt $ writeFile ("/tmp/cublas-"++typ++".txt") $ show yinvs
-
-
-main :: IO ()
-main = do
-   let n = 96
-   let sz = 50
-   let its = 20
-   let xmsDouble = randomMatrixInv sz
-       ysDouble = Rnd.randomRs (-0.01,0.01::Double) $ Rnd.mkStdGen 23
-   let xmsFloat =
-          mapPair
-             (Container.cmap realToFrac, Container.cmap realToFrac)
-             xmsDouble
-       ysFloat :: [Float]
-       ysFloat = map realToFrac ysDouble
-   mainHMatrixDirect "double" n xmsDouble ysDouble
---   mainHMatrixDirect "float" n xmsFloat ysFloat
-   mainHMatrix "double" n its xmsDouble ysDouble
-   mainHMatrix "float" n its xmsFloat ysFloat
-   mainCUBLASDirect "double" n xmsDouble ysDouble
-   mainCUBLASDirect "float" n xmsFloat ysFloat
-   mainCUBLAS "double" n its xmsDouble ysDouble
-   mainCUBLAS "float" n its xmsFloat ysFloat
-   mainCUDA "double" n its xmsDouble ysDouble
-   mainCUDA "float" n its xmsFloat ysFloat
diff --git a/private/Data/Array/Accelerate/LinearAlgebra/Private.hs b/private/Data/Array/Accelerate/LinearAlgebra/Private.hs
new file mode 100644
--- /dev/null
+++ b/private/Data/Array/Accelerate/LinearAlgebra/Private.hs
@@ -0,0 +1,254 @@
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
+module Data.Array.Accelerate.LinearAlgebra.Private where
+
+import qualified Data.Array.Accelerate.Utility.Loop as Loop
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate
+          (Acc, Array, Exp, Any(Any), All(All), Z(Z), (:.)((:.)))
+
+
+
+type Scalar ix a = Acc (Array ix a)
+type Vector ix a = Acc (Array (ix :. Int) a)
+type Matrix ix a = Acc (Array (ix :. Int :. Int) a)
+
+transpose ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Matrix ix a -> Matrix ix a
+transpose m =
+   A.backpermute
+      (A.lift $ swapIndex $ matrixShape m)
+      (A.lift . swapIndex . A.unlift)
+      m
+
+swapIndex ::
+   Exp ix :. Exp Int :. Exp Int ->
+   Exp ix :. Exp Int :. Exp Int
+swapIndex (ix :. r :. c) = (ix :. c :. r)
+
+
+numElems :: (A.Shape ix, A.Slice ix, A.Elt a) => Vector ix a -> Exp Int
+numElems m = case vectorShape m of _ix :. n -> n
+
+numRows :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int
+numRows m = case matrixShape m of _ix :. rows :. _cols -> rows
+
+numCols :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int
+numCols m = case matrixShape m of _ix :. _rows :. cols -> cols
+
+vectorShape ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Vector ix a -> Exp ix :. Exp Int
+vectorShape m = A.unlift $ A.shape m
+
+matrixShape ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Matrix ix a -> Exp ix :. Exp Int :. Exp Int
+matrixShape m = A.unlift $ A.shape m
+
+withVectorIndex ::
+   (A.Shape ix, A.Slice ix, A.Lift Exp a) =>
+   (Exp ix :. Exp Int -> a) ->
+   (Exp (ix :. Int) -> Exp (A.Plain a))
+withVectorIndex f = A.lift . f . A.unlift
+
+withMatrixIndex ::
+   (A.Shape ix, A.Slice ix, A.Lift Exp a) =>
+   (Exp ix :. Exp Int :. Exp Int -> a) ->
+   (Exp (ix :. Int :. Int) -> Exp (A.Plain a))
+withMatrixIndex f = A.lift . f . A.unlift
+
+
+outer ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Vector ix a -> Vector ix a -> Matrix ix a
+outer x y =
+   A.zipWith (*)
+      (A.replicate (A.lift $ Any :. All :. numElems y) x)
+      (A.replicate (A.lift $ Any :. numElems x :. All) y)
+
+multiplyMatrixVector ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Matrix ix a ->
+   Vector ix a ->
+   Vector ix a
+multiplyMatrixVector m v =
+   case matrixShape m of
+      (_ix :. rows :. _cols) ->
+         A.fold1 (+) $
+         A.zipWith (*) m
+            (A.replicate (A.lift $ Any :. rows :. All) v)
+
+multiplyMatrixMatrix ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Matrix ix a ->
+   Matrix ix a ->
+   Matrix ix a
+multiplyMatrixMatrix x y =
+   case (matrixShape x, matrixShape y) of
+      (_ :. rows :. _cols, _ :. _rows :. cols) ->
+         A.fold1 (+) $ transpose $
+         A.zipWith (*)
+            (A.replicate (A.lift $ Any :. All :. All :. cols) x)
+            (A.replicate (A.lift $ Any :. rows :. All :. All) y)
+
+newtonInverseStep ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Matrix ix a ->
+   Matrix ix a ->
+   Matrix ix a
+newtonInverseStep a x =
+   A.zipWith (-) (A.map (2*) x) $
+   multiplyMatrixMatrix x $ multiplyMatrixMatrix a x
+
+identity ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Exp (ix :. Int :. Int) -> Matrix ix a
+identity sh =
+   A.generate sh
+      (withMatrixIndex $
+       \(_ :. r :. c) -> A.fromIntegral $ A.boolToInt (r A.==* c))
+
+newtonInverse ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Exp Int ->
+   Matrix ix a ->
+   Matrix ix a ->
+   Matrix ix a
+newtonInverse n seed a =
+   Loop.nest n (newtonInverseStep a) seed
+
+
+
+scaleRows ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
+   Vector ix a -> Matrix ix a -> Matrix ix a
+scaleRows s x =
+   zipScalarVectorWith (*) s x
+
+
+
+zipScalarVectorWith ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   Scalar ix a -> Vector ix b -> Vector ix c
+zipScalarVectorWith f x ys =
+   case vectorShape ys of
+      _ix :. dim ->
+         A.zipWith f (A.replicate (A.lift (Any :. dim)) x) ys
+
+zipScalarMatrixWith ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   Scalar ix a -> Matrix ix b -> Matrix ix c
+zipScalarMatrixWith f x ys =
+   case matrixShape ys of
+      _ix :. rows :. cols ->
+         A.zipWith f
+            (A.replicate (A.lift (Any :. rows :. cols)) x) ys
+
+
+
+columnFromVector ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Vector ix a -> Matrix ix a
+columnFromVector a = A.reshape (Exp.indexCons (A.shape a) 1) a
+
+{- |
+input must be a matrix with exactly one column
+-}
+vectorFromColumn ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Matrix ix a -> Vector ix a
+vectorFromColumn a = A.reshape (A.indexTail $ A.shape a) a
+
+
+
+flattenMatrix, flattenMatrixReshape, flattenMatrixBackPermute ::
+   (A.Slice ix, A.Shape ix, A.Elt a) =>
+   Matrix ix a -> Vector ix a
+flattenMatrix = flattenMatrixBackPermute
+
+flattenMatrixReshape m =
+   case matrixShape m of
+      ix :. rows :. cols ->
+         A.reshape (A.lift $ ix :. rows*cols) m
+
+accDivMod :: Integral a => a -> a -> (a, a)
+accDivMod x y = (div x y, mod x y)
+
+flattenMatrixBackPermute m =
+   case matrixShape m of
+      ix :. rows :. cols ->
+         A.backpermute
+            (A.lift $ ix :. rows*cols)
+            (withVectorIndex $
+             \(vix :. n) -> case accDivMod n cols of (r,c) -> vix :. r :. c)
+            m
+
+
+restoreMatrix, restoreMatrixReshape, restoreMatrixBackPermute ::
+   (A.Slice ix, A.Shape ix, A.Elt a) =>
+   Exp Int -> Vector ix a -> Matrix ix a
+restoreMatrix = restoreMatrixBackPermute
+
+restoreMatrixReshape cols v =
+   case vectorShape v of
+      ix :. n ->
+         A.reshape (A.lift $ ix :. div n cols :. cols) v
+
+restoreMatrixBackPermute cols v =
+   case vectorShape v of
+      ix :. n ->
+         A.backpermute
+            (A.lift $ ix :. div n cols :. cols)
+            (withMatrixIndex $ \(vix :. k :. j) -> vix :. k*cols+j)
+            v
+
+
+
+extrudeVector ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Exp ix -> Vector Z a -> Vector ix a
+extrudeVector shape y =
+   -- A.replicate (A.lift $ shape :. All) y
+   A.backpermute
+      (A.lift $ shape :. numElems y)
+      (A.index1 . A.indexHead)
+      y
+
+extrudeMatrix ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Exp ix -> Matrix Z a -> Matrix ix a
+extrudeMatrix shape y =
+   A.backpermute
+      (A.lift $ shape :. numRows y :. numCols y)
+      (withMatrixIndex $ \(_:.r:.c) -> Z:.r:.c)
+      y
+
+zipExtrudedVectorWith ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   Vector Z a ->
+   Vector ix b ->
+   Vector ix c
+zipExtrudedVectorWith f x y =
+   A.zipWith f (extrudeVector (A.indexTail $ A.shape y) x) y
+
+zipExtrudedMatrixWith ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   Matrix Z a ->
+   Matrix ix b ->
+   Matrix ix c
+zipExtrudedMatrixWith f x y =
+   A.zipWith f (extrudeMatrix (A.indexTail $ A.indexTail $ A.shape y) x) y
+
+gatherFromVector ::
+   (A.Shape ix, A.Elt a) =>
+   Scalar ix Int -> Vector Z a -> Scalar ix a
+gatherFromVector indices =
+   Arrange.gather (A.map A.index1 indices)
diff --git a/src/Data/Array/Accelerate/Arithmetic/Example.hs b/src/Data/Array/Accelerate/Arithmetic/Example.hs
--- a/src/Data/Array/Accelerate/Arithmetic/Example.hs
+++ b/src/Data/Array/Accelerate/Arithmetic/Example.hs
@@ -1,8 +1,8 @@
 module Data.Array.Accelerate.Arithmetic.Example where
 
 import qualified Data.Array.Accelerate.Arithmetic.Interpolation as Ip
-import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse
-import Data.Array.Accelerate.Arithmetic.LinearAlgebra (Vector, )
+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse as Sparse
+import Data.Array.Accelerate.LinearAlgebra (Vector, )
 
 import qualified Data.Array.Accelerate.Interpreter as AI
 import qualified Data.Array.Accelerate as A
@@ -11,9 +11,9 @@
 
 exampleSparseColumnMatrix :: IO ()
 exampleSparseColumnMatrix = do
-   let m :: Sparse.ColumnMatrix Z Double
+   let m :: Sparse.Columns Z Double
        m =
-          Sparse.ColumnMatrix (A.lift (3::Int)) $
+          Sparse.Columns (A.lift (3::Int)) $
           A.use $ A.fromList (Z :. 2 :. 5) $
           (0,1) : (2,2) : (1,3) : (0,4) : (2,5) :
           (1,6) : (2,7) : (0,8) : (2,9) : (1,10) :
@@ -22,13 +22,13 @@
        v :: Vector Z Double
        v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000]
 
-   print $ AI.run $ Sparse.multiplyColumnMatrixVector m v
+   print $ AI.run $ Sparse.multiplyColumnsVector m v
 
 exampleSparseRowMatrix :: IO ()
 exampleSparseRowMatrix = do
-   let m :: Sparse.RowMatrix Z Double
+   let m :: Sparse.Rows Z Double
        m =
-          Sparse.RowMatrix (A.lift (5::Int)) $
+          Sparse.Rows (A.lift (5::Int)) $
           A.use $ A.fromList (Z :. 3 :. 2) $
           (0,1) : (0,2) :
           (3,3) : (1,4) :
@@ -38,7 +38,7 @@
        v :: Vector Z Double
        v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000]
 
-   print $ AI.run $ Sparse.multiplyRowMatrixVector m v
+   print $ AI.run $ Sparse.multiplyRowsVector m v
 
 exampleLookup :: IO ()
 exampleLookup = do
diff --git a/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs b/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs
--- a/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs
+++ b/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs
@@ -4,12 +4,12 @@
    Interpolator13, sampleBasisFunctions13,
    ) where
 
-import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse as Sparse
+import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg
 import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
 import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
 import qualified Data.Array.Accelerate.Utility.Loop as Loop
-import Data.Array.Accelerate.Arithmetic.LinearAlgebra
+import Data.Array.Accelerate.LinearAlgebra
           (Scalar, Vector, numElems, extrudeVector, )
 
 import qualified Data.Array.Accelerate as A
@@ -68,9 +68,9 @@
 sampleBasisFunctions13 ::
    (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a, Num a) =>
    Interpolator13 (Exp a) ->
-   Vector Z a -> Vector ix a -> Sparse.RowMatrix ix a
+   Vector Z a -> Vector ix a -> Sparse.Rows ix a
 sampleBasisFunctions13 interpolate nodes zs =
-   Sparse.RowMatrix (numElems nodes) $
+   Sparse.Rows (numElems nodes) $
    let indices = lookupInterval (extrudeVector (A.shape zs) nodes) zs
        minIx = 1
        maxIx = numElems nodes - 3
diff --git a/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs b/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
deleted file mode 100644
--- a/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
+++ /dev/null
@@ -1,254 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE FlexibleContexts #-}
-module Data.Array.Accelerate.Arithmetic.LinearAlgebra where
-
-import qualified Data.Array.Accelerate.Utility.Loop as Loop
-import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
-import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
-import qualified Data.Array.Accelerate as A
-import Data.Array.Accelerate
-          (Acc, Array, Exp, Any(Any), All(All), Z(Z), (:.)((:.)))
-
-
-
-type Scalar ix a = Acc (Array ix a)
-type Vector ix a = Acc (Array (ix :. Int) a)
-type Matrix ix a = Acc (Array (ix :. Int :. Int) a)
-
-transpose ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Matrix ix a -> Matrix ix a
-transpose m =
-   A.backpermute
-      (A.lift $ swapIndex $ matrixShape m)
-      (A.lift . swapIndex . A.unlift)
-      m
-
-swapIndex ::
-   Exp ix :. Exp Int :. Exp Int ->
-   Exp ix :. Exp Int :. Exp Int
-swapIndex (ix :. r :. c) = (ix :. c :. r)
-
-
-numElems :: (A.Shape ix, A.Slice ix, A.Elt a) => Vector ix a -> Exp Int
-numElems m = case vectorShape m of _ix :. n -> n
-
-numRows :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int
-numRows m = case matrixShape m of _ix :. rows :. _cols -> rows
-
-numCols :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int
-numCols m = case matrixShape m of _ix :. _rows :. cols -> cols
-
-vectorShape ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Vector ix a -> Exp ix :. Exp Int
-vectorShape m = A.unlift $ A.shape m
-
-matrixShape ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Matrix ix a -> Exp ix :. Exp Int :. Exp Int
-matrixShape m = A.unlift $ A.shape m
-
-withVectorIndex ::
-   (A.Shape ix, A.Slice ix, A.Lift Exp a) =>
-   (Exp ix :. Exp Int -> a) ->
-   (Exp (ix :. Int) -> Exp (A.Plain a))
-withVectorIndex f = A.lift . f . A.unlift
-
-withMatrixIndex ::
-   (A.Shape ix, A.Slice ix, A.Lift Exp a) =>
-   (Exp ix :. Exp Int :. Exp Int -> a) ->
-   (Exp (ix :. Int :. Int) -> Exp (A.Plain a))
-withMatrixIndex f = A.lift . f . A.unlift
-
-
-outer ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Vector ix a -> Vector ix a -> Matrix ix a
-outer x y =
-   A.zipWith (*)
-      (A.replicate (A.lift $ Any :. All :. numElems y) x)
-      (A.replicate (A.lift $ Any :. numElems x :. All) y)
-
-multiplyMatrixVector ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Matrix ix a ->
-   Vector ix a ->
-   Vector ix a
-multiplyMatrixVector m v =
-   case matrixShape m of
-      (_ix :. rows :. _cols) ->
-         A.fold1 (+) $
-         A.zipWith (*) m
-            (A.replicate (A.lift $ Any :. rows :. All) v)
-
-multiplyMatrixMatrix ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Matrix ix a ->
-   Matrix ix a ->
-   Matrix ix a
-multiplyMatrixMatrix x y =
-   case (matrixShape x, matrixShape y) of
-      (_ :. rows :. _cols, _ :. _rows :. cols) ->
-         A.fold1 (+) $ transpose $
-         A.zipWith (*)
-            (A.replicate (A.lift $ Any :. All :. All :. cols) x)
-            (A.replicate (A.lift $ Any :. rows :. All :. All) y)
-
-newtonInverseStep ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Matrix ix a ->
-   Matrix ix a ->
-   Matrix ix a
-newtonInverseStep a x =
-   A.zipWith (-) (A.map (2*) x) $
-   multiplyMatrixMatrix x $ multiplyMatrixMatrix a x
-
-identity ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Exp (ix :. Int :. Int) -> Matrix ix a
-identity sh =
-   A.generate sh
-      (withMatrixIndex $
-       \(_ :. r :. c) -> A.fromIntegral $ A.boolToInt (r A.==* c))
-
-newtonInverse ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Exp Int ->
-   Matrix ix a ->
-   Matrix ix a ->
-   Matrix ix a
-newtonInverse n seed a =
-   Loop.nest n (newtonInverseStep a) seed
-
-
-
-scaleRows ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
-   Vector ix a -> Matrix ix a -> Matrix ix a
-scaleRows s x =
-   zipScalarVectorWith (*) s x
-
-
-
-zipScalarVectorWith ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
-   (Exp a -> Exp b -> Exp c) ->
-   Scalar ix a -> Vector ix b -> Vector ix c
-zipScalarVectorWith f x ys =
-   case vectorShape ys of
-      _ix :. dim ->
-         A.zipWith f (A.replicate (A.lift (Any :. dim)) x) ys
-
-zipScalarMatrixWith ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
-   (Exp a -> Exp b -> Exp c) ->
-   Scalar ix a -> Matrix ix b -> Matrix ix c
-zipScalarMatrixWith f x ys =
-   case matrixShape ys of
-      _ix :. rows :. cols ->
-         A.zipWith f
-            (A.replicate (A.lift (Any :. rows :. cols)) x) ys
-
-
-
-columnFromVector ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Vector ix a -> Matrix ix a
-columnFromVector a = A.reshape (Exp.indexCons (A.shape a) 1) a
-
-{- |
-input must be a matrix with exactly one column
--}
-vectorFromColumn ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Matrix ix a -> Vector ix a
-vectorFromColumn a = A.reshape (A.indexTail $ A.shape a) a
-
-
-
-flattenMatrix, flattenMatrixReshape, flattenMatrixBackPermute ::
-   (A.Slice ix, A.Shape ix, A.Elt a) =>
-   Matrix ix a -> Vector ix a
-flattenMatrix = flattenMatrixBackPermute
-
-flattenMatrixReshape m =
-   case matrixShape m of
-      ix :. rows :. cols ->
-         A.reshape (A.lift $ ix :. rows*cols) m
-
-accDivMod :: Integral a => a -> a -> (a, a)
-accDivMod x y = (div x y, mod x y)
-
-flattenMatrixBackPermute m =
-   case matrixShape m of
-      ix :. rows :. cols ->
-         A.backpermute
-            (A.lift $ ix :. rows*cols)
-            (withVectorIndex $
-             \(vix :. n) -> case accDivMod n cols of (r,c) -> vix :. r :. c)
-            m
-
-
-restoreMatrix, restoreMatrixReshape, restoreMatrixBackPermute ::
-   (A.Slice ix, A.Shape ix, A.Elt a) =>
-   Exp Int -> Vector ix a -> Matrix ix a
-restoreMatrix = restoreMatrixBackPermute
-
-restoreMatrixReshape cols v =
-   case vectorShape v of
-      ix :. n ->
-         A.reshape (A.lift $ ix :. div n cols :. cols) v
-
-restoreMatrixBackPermute cols v =
-   case vectorShape v of
-      ix :. n ->
-         A.backpermute
-            (A.lift $ ix :. div n cols :. cols)
-            (withMatrixIndex $ \(vix :. k :. j) -> vix :. k*cols+j)
-            v
-
-
-
-extrudeVector ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Exp ix -> Vector Z a -> Vector ix a
-extrudeVector shape y =
-   -- A.replicate (A.lift $ shape :. All) y
-   A.backpermute
-      (A.lift $ shape :. numElems y)
-      (A.index1 . A.indexHead)
-      y
-
-extrudeMatrix ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Exp ix -> Matrix Z a -> Matrix ix a
-extrudeMatrix shape y =
-   A.backpermute
-      (A.lift $ shape :. numRows y :. numCols y)
-      (withMatrixIndex $ \(_:.r:.c) -> Z:.r:.c)
-      y
-
-zipExtrudedVectorWith ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
-   (Exp a -> Exp b -> Exp c) ->
-   Vector Z a ->
-   Vector ix b ->
-   Vector ix c
-zipExtrudedVectorWith f x y =
-   A.zipWith f (extrudeVector (A.indexTail $ A.shape y) x) y
-
-zipExtrudedMatrixWith ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>
-   (Exp a -> Exp b -> Exp c) ->
-   Matrix Z a ->
-   Matrix ix b ->
-   Matrix ix c
-zipExtrudedMatrixWith f x y =
-   A.zipWith f (extrudeMatrix (A.indexTail $ A.indexTail $ A.shape y) x) y
-
-gatherFromVector ::
-   (A.Shape ix, A.Elt a) =>
-   Scalar ix Int -> Vector Z a -> Scalar ix a
-gatherFromVector indices =
-   Arrange.gather (A.map A.index1 indices)
diff --git a/src/Data/Array/Accelerate/Arithmetic/Sparse.hs b/src/Data/Array/Accelerate/Arithmetic/Sparse.hs
deleted file mode 100644
--- a/src/Data/Array/Accelerate/Arithmetic/Sparse.hs
+++ /dev/null
@@ -1,115 +0,0 @@
-{-# LANGUAGE TypeOperators #-}
-module Data.Array.Accelerate.Arithmetic.Sparse where
-
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
-import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
-import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
-import qualified Data.Array.Accelerate as A
-import Data.Array.Accelerate.Utility.Lift.Exp (expr, )
-
-import Data.Array.Accelerate.Arithmetic.LinearAlgebra
-          (Matrix, Vector, matrixShape, )
-import Data.Array.Accelerate
-          (Exp, Any(Any), All(All), (:.)((:.)), )
-
-
-{- |
-Sparse matrix with a definite number of non-zero entries per column.
--}
-data ColumnMatrix ix a =
-        ColumnMatrix {numRows :: Exp Int, columnMatrix :: Matrix ix (Int, a)}
-
-realIndex ::
-   (A.Shape ix, A.Slice ix, A.Elt a) =>
-   Matrix ix (Int, a) ->
-   Matrix ix (ix :. Int)
-realIndex m =
-   A.zipWith Exp.indexCons
-      (A.generate (A.shape m) (A.indexTail . A.indexTail))
-      (A.map A.fst m)
-
-multiplyColumnMatrixVector ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   ColumnMatrix ix a ->
-   Vector ix a ->
-   Vector ix a
-multiplyColumnMatrixVector (ColumnMatrix rows m) v =
-   Arrange.scatter (+)
-      (realIndex m)
-      (case matrixShape m of
-          sh :. _rows :. _cols -> A.fill (A.lift $ sh :. rows) 0) $
-   A.zipWith (*)
-      (A.map A.snd m)
-      (A.replicate (A.lift $ Any :. LinAlg.numRows m :. All) v)
-
-transposeColumnMatrix ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   ColumnMatrix ix a ->
-   RowMatrix ix a
-transposeColumnMatrix (ColumnMatrix n x) =
-   RowMatrix n $ LinAlg.transpose x
-
-
-{- |
-Sparse matrix with a definite number of non-zero entries per row.
--}
-data RowMatrix ix a =
-        RowMatrix {numCols :: Exp Int, rowMatrix :: Matrix ix (Int, a)}
-
-multiplyRowMatrixVector ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   RowMatrix ix a ->
-   Vector ix a ->
-   Vector ix a
-multiplyRowMatrixVector (RowMatrix _cols m) v =
-   A.fold1 (+) $
-   A.zipWith (*) (A.map A.snd m) $
-   Arrange.gather (realIndex m) v
-
-transposeRowMatrix ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   RowMatrix ix a ->
-   ColumnMatrix ix a
-transposeRowMatrix (RowMatrix n x) =
-   (ColumnMatrix n $ LinAlg.transpose x)
-
-multiplyMatrixMatrix ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   ColumnMatrix ix a ->
-   RowMatrix ix a ->
-   Matrix ix a
-multiplyMatrixMatrix
-      (ColumnMatrix rows x) (RowMatrix cols y) =
-   case matchMatrices x y of
-      m ->
-         let global = A.indexTail . A.indexTail . A.indexTail
-         in  Arrange.scatter (+)
-                (Arrange.mapWithIndex
-                   (\mix tix ->
-                      A.lift $ global mix :. A.fst tix :. A.snd tix) $
-                 A.map A.fst m)
-                (A.fill (A.lift $ global (A.shape m) :. rows :. cols) 0)
-                (A.map A.snd m)
-
-matchMatrices ::
-   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
-   Matrix ix (Int, a) ->
-   Matrix ix (Int, a) ->
-   Matrix (ix :. Int) ((Int, Int), a)
-matchMatrices x y =
-   case (matrixShape x, matrixShape y) of
-      (_ :. xRows :. _xCols, _ :. _yRows :. yCols) ->
-         -- it must be xCols == yRows
-         A.zipWith
-            (Exp.modify2 (expr,expr) (expr,expr) $
-             \(n,xi) (m,yi) -> ((n, m), xi*yi))
-            (A.replicate (A.lift $ Any :. All :. All :. yCols) x)
-            (A.replicate (A.lift $ Any :. xRows :. All :. All) y)
-
-
-scaleRowRows ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
-   Vector ix a -> RowMatrix ix a -> RowMatrix ix a
-scaleRowRows s (RowMatrix n x) =
-   RowMatrix n $
-   LinAlg.zipScalarVectorWith (\si xi -> Exp.mapSnd (si*) xi) s x
diff --git a/src/Data/Array/Accelerate/LinearAlgebra.hs b/src/Data/Array/Accelerate/LinearAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/LinearAlgebra.hs
@@ -0,0 +1,35 @@
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
+module Data.Array.Accelerate.LinearAlgebra (
+   LinAlg.Scalar,
+   LinAlg.Vector,
+   LinAlg.Matrix,
+   LinAlg.transpose,
+   LinAlg.numElems,
+   LinAlg.numRows,
+   LinAlg.numCols,
+   LinAlg.vectorShape,
+   LinAlg.matrixShape,
+   LinAlg.withVectorIndex,
+   LinAlg.withMatrixIndex,
+   LinAlg.outer,
+   LinAlg.multiplyMatrixVector,
+   LinAlg.multiplyMatrixMatrix,
+   LinAlg.newtonInverse,
+   LinAlg.newtonInverseStep,
+   LinAlg.identity,
+   LinAlg.scaleRows,
+   LinAlg.zipScalarVectorWith,
+   LinAlg.zipScalarMatrixWith,
+   LinAlg.columnFromVector,
+   LinAlg.vectorFromColumn,
+   LinAlg.flattenMatrix,
+   LinAlg.restoreMatrix,
+   LinAlg.extrudeVector,
+   LinAlg.extrudeMatrix,
+   LinAlg.zipExtrudedVectorWith,
+   LinAlg.zipExtrudedMatrixWith,
+   LinAlg.gatherFromVector,
+   ) where
+
+import qualified Data.Array.Accelerate.LinearAlgebra.Private as LinAlg
diff --git a/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Banded.hs b/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Banded.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Banded.hs
@@ -0,0 +1,26 @@
+module Data.Array.Accelerate.LinearAlgebra.Matrix.Banded (
+   Symmetric(..),
+   flattenSymmetric,
+   ) where
+
+import Data.Array.Accelerate.LinearAlgebra (Matrix, matrixShape)
+
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Utility.Lift.Exp (expr)
+import Data.Array.Accelerate ((:.)((:.)), (>*), (!), (?))
+
+
+newtype Symmetric ix a = Symmetric (Matrix ix a)
+
+flattenSymmetric ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
+   Symmetric ix a -> Matrix ix a
+flattenSymmetric (Symmetric m) =
+   case matrixShape m of
+      (sh :. rows :. width) ->
+         A.generate (A.lift $ sh :. rows :. rows) $
+         Exp.modify (expr:.expr:.expr) $ \(ix:.k0:.j0) ->
+            let k = min k0 j0
+                j = max k0 j0 - k
+            in  width >* j ? (m ! A.lift(ix:.k:.j), 0)
diff --git a/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Sparse.hs b/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Accelerate/LinearAlgebra/Matrix/Sparse.hs
@@ -0,0 +1,148 @@
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse (
+   Columns(..),
+   multiplyColumnsVector,
+   transposeColumns,
+   Rows(..),
+   multiplyRowsVector,
+   transposeRows,
+   multiplyColumnsRows,
+   realBandedGramian,
+   scaleRowRows,
+   ) where
+
+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Banded as BandMatrix
+import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg
+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange
+import qualified Data.Array.Accelerate as A
+import Data.Array.Accelerate.Utility.Lift.Exp (expr, )
+
+import Data.Array.Accelerate.LinearAlgebra
+          (Matrix, Vector, matrixShape, )
+import Data.Array.Accelerate
+          (Exp, Any(Any), All(All), (:.)((:.)), (>*), (?), )
+
+
+{- |
+Sparse matrix with a definite number of non-zero entries per column.
+-}
+data Columns ix a =
+        Columns {numRows :: Exp Int, columnMatrix :: Matrix ix (Int, a)}
+
+realIndex ::
+   (A.Shape ix, A.Slice ix, A.Elt a) =>
+   Matrix ix (Int, a) ->
+   Matrix ix (ix :. Int)
+realIndex m =
+   A.zipWith Exp.indexCons
+      (A.generate (A.shape m) (A.indexTail . A.indexTail))
+      (A.map A.fst m)
+
+multiplyColumnsVector ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Columns ix a ->
+   Vector ix a ->
+   Vector ix a
+multiplyColumnsVector (Columns rows m) v =
+   Arrange.scatter (+)
+      (realIndex m)
+      (case matrixShape m of
+          sh :. _rows :. _cols -> A.fill (A.lift $ sh :. rows) 0) $
+   A.zipWith (*)
+      (A.map A.snd m)
+      (A.replicate (A.lift $ Any :. LinAlg.numRows m :. All) v)
+
+transposeColumns ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Columns ix a ->
+   Rows ix a
+transposeColumns (Columns n x) =
+   Rows n $ LinAlg.transpose x
+
+
+{- |
+Sparse matrix with a definite number of non-zero entries per row.
+-}
+data Rows ix a =
+        Rows {numCols :: Exp Int, rowMatrix :: Matrix ix (Int, a)}
+
+multiplyRowsVector ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Rows ix a ->
+   Vector ix a ->
+   Vector ix a
+multiplyRowsVector (Rows _cols m) v =
+   A.fold1 (+) $
+   A.zipWith (*) (A.map A.snd m) $
+   Arrange.gather (realIndex m) v
+
+transposeRows ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Rows ix a ->
+   Columns ix a
+transposeRows (Rows n x) =
+   (Columns n $ LinAlg.transpose x)
+
+multiplyColumnsRows ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Columns ix a ->
+   Rows ix a ->
+   Matrix ix a
+multiplyColumnsRows (Columns rows x) (Rows cols y) =
+   let (ixs,prods) = A.unzip $ matchMatrices x y
+       global = A.indexTail . A.indexTail . A.indexTail
+   in  Arrange.scatter (+)
+          (Arrange.mapWithIndex
+             (Exp.modify2 expr (expr,expr) $ \mix (k,j) ->
+                global mix :. k :. j) $
+           ixs)
+          (A.fill (A.lift $ global (A.shape prods) :. rows :. cols) 0)
+          prods
+
+{- |
+Compute x^T*x, given that it has a band structure.
+You must pass the band-width as parameter
+and you must make sure that the Gramian stays within this band.
+Otherwise you cause out-of-bounds array accesses.
+So far, only correct for real matrices.
+-}
+realBandedGramian ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Exp Int ->
+   Rows ix a ->
+   BandMatrix.Symmetric ix a
+realBandedGramian width (Rows cols y) =
+   let (ixs,prods) = A.unzip $ matchMatrices (LinAlg.transpose y) y
+       global = A.indexTail . A.indexTail . A.indexTail
+   in  BandMatrix.Symmetric $
+       Arrange.scatter (+)
+          (Arrange.mapWithIndex
+             (Exp.modify2 expr (expr,expr) $ \mix (k,j) ->
+                k>*j ? (A.ignore, A.lift $ global mix :. k :. j-k)) $
+           ixs)
+          (A.fill (A.lift $ global (A.shape prods) :. cols :. width) 0)
+          prods
+
+matchMatrices ::
+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>
+   Matrix ix (Int, a) ->
+   Matrix ix (Int, a) ->
+   Matrix (ix :. Int) ((Int, Int), a)
+matchMatrices x y =
+   case (matrixShape x, matrixShape y) of
+      (_ :. xRows :. _xCols, _ :. _yRows :. yCols) ->
+         -- it must be xCols == yRows
+         A.zipWith
+            (Exp.modify2 (expr,expr) (expr,expr) $
+             \(n,xi) (m,yi) -> ((n, m), xi*yi))
+            (A.replicate (A.lift $ Any :. All :. All :. yCols) x)
+            (A.replicate (A.lift $ Any :. xRows :. All :. All) y)
+
+
+scaleRowRows ::
+   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
+   Vector ix a -> Rows ix a -> Rows ix a
+scaleRowRows s (Rows n x) =
+   Rows n $
+   LinAlg.zipScalarVectorWith (\si xi -> Exp.mapSnd (si*) xi) s x
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -10,6 +10,7 @@
 test :: IO ()
 test = mapM_ (\(msg,act) -> putStr (msg++": ") >> act) $
    ("sparseMatrix", quickCheck (\(Mod.Blind x) -> Sparse.multiplication x)) :
+   ("bandedGramian", quickCheck (\(Mod.Blind x) -> Sparse.bandedGramian x)) :
    ("flattenMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.flattenMatrix x)) :
    ("restoreMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.restoreMatrix x)) :
    ("flattenRestoreMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.flattenRestoreMatrix x)) :
diff --git a/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs b/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
--- a/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
+++ b/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
@@ -1,9 +1,9 @@
 module Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra where
 
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra.Private as LinAlg
 import qualified Data.Array.Accelerate as A
 
-import Data.Array.Accelerate.Arithmetic.LinearAlgebra (Matrix, numCols, )
+import Data.Array.Accelerate.LinearAlgebra.Private (Matrix, numCols, )
 import Data.Array.Accelerate (Z(Z), (:.)((:.)),)
 
 import Test.Data.Array.Accelerate.Arithmetic.Utility (arbitraryArray, (=!=), )
diff --git a/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs b/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs
--- a/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs
+++ b/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs
@@ -4,11 +4,12 @@
 
 import Test.Data.Array.Accelerate.Arithmetic.Utility (arbitraryArray, (=!=), )
 
-import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Banded as BandMatrix
+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse as Sparse
+import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg
 
 import qualified Data.Array.Accelerate as A
-import Data.Array.Accelerate (Z(Z), (:.)((:.)))
+import Data.Array.Accelerate (Exp, Z(Z), (:.)((:.)))
 
 import qualified Test.QuickCheck as QC
 
@@ -20,8 +21,8 @@
 data
    CRVTriple a =
       CRVTriple
-         (Sparse.ColumnMatrix Z a)
-         (Sparse.RowMatrix Z a)
+         (Sparse.Columns Z a)
+         (Sparse.Rows Z a)
          (LinAlg.Vector Z a)
 
 instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (CRVTriple a) where
@@ -40,13 +41,42 @@
       v <- arbitraryArray (Z :. nr) QC.arbitrary
       return $
          CRVTriple
-            (Sparse.ColumnMatrix (A.lift nc) (A.use mc))
-            (Sparse.RowMatrix (A.lift nr) (A.use mr))
+            (Sparse.Columns (A.lift nc) (A.use mc))
+            (Sparse.Rows (A.lift nr) (A.use mr))
             (A.use v)
 
 
 multiplication :: CRVTriple Word32 -> Bool
 multiplication (CRVTriple mc mr v) =
-   LinAlg.multiplyMatrixVector (Sparse.multiplyMatrixMatrix mc mr) v
+   LinAlg.multiplyMatrixVector (Sparse.multiplyColumnsRows mc mr) v
    =!=
-   Sparse.multiplyColumnMatrixVector mc (Sparse.multiplyRowMatrixVector mr v)
+   Sparse.multiplyColumnsVector mc (Sparse.multiplyRowsVector mr v)
+
+
+
+data BandGramian a = BandGramian (Exp Int) (Sparse.Rows Z a)
+
+instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (BandGramian a) where
+   arbitrary = do
+      width <- QC.choose (1,10)
+      rows <- QC.choose (1,100)
+      cols <- QC.choose (width,100)
+
+      m <-
+         fmap (A.fromList (Z :. rows :. width) . concat) $
+         QC.vectorOf rows $
+         liftM2
+            (\start row -> zip [start..] row)
+            (QC.choose (0,cols-width))
+            (QC.vectorOf width QC.arbitrary)
+
+      return $
+         BandGramian (A.lift width)
+            (Sparse.Rows (A.lift cols) (A.use m))
+
+
+bandedGramian :: BandGramian Word32 -> Bool
+bandedGramian (BandGramian width m) =
+   Sparse.multiplyColumnsRows (Sparse.transposeRows m) m
+   =!=
+   BandMatrix.flattenSymmetric (Sparse.realBandedGramian width m)