diff --git a/accelerate-cublas.cabal b/accelerate-cublas.cabal
--- a/accelerate-cublas.cabal
+++ b/accelerate-cublas.cabal
@@ -1,10 +1,10 @@
 Name:             accelerate-cublas
-Version:          0.0
+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-cublas/
+Homepage:         http://hub.darcs.net/thielema/accelerate-cublas/
 Category:         Math
 Synopsis:         Basic Linear Algebra using native CUBLAS library
 Description:
@@ -19,26 +19,26 @@
   default:     False
 
 Source-Repository this
-  Tag:         0.0
+  Tag:         0.1
   Type:        darcs
-  Location:    http://code.haskell.org/~thielema/accelerate-cublas/
+  Location:    http://hub.darcs.net/thielema/accelerate-cublas/
 
 Source-Repository head
   Type:        darcs
-  Location:    http://code.haskell.org/~thielema/accelerate-cublas/
+  Location:    http://hub.darcs.net/thielema/accelerate-cublas/
 
 Library
   Build-Depends:
-    accelerate-arithmetic >=0.0.1 && <0.1,
+    accelerate-arithmetic >=0.1 && <0.2,
     accelerate-utility >=0.1 && <0.2,
-    accelerate-cuda >=0.15 && <0.16,
+    accelerate-cuda >=0.16 && <0.17,
     accelerate-io >=0.15 && <0.16,
     accelerate >=0.15 && <0.16,
     cublas >=0.2.0.2 && <0.3,
-    cuda >=0.5 && <0.7,
-    vector >=0.10.11 && <0.11,
+    cuda >=0.5 && <0.8,
+    vector >=0.10.11 && <0.12,
     utility-ht >=0.0.8 && <0.1,
-    base >=4.5 && <4.8
+    base >=4.5 && <5
 
   GHC-Options:      -Wall -fwarn-missing-import-lists
   Hs-Source-Dirs:   src
@@ -60,6 +60,26 @@
     accelerate-arithmetic,
     accelerate-utility,
     accelerate,
-    cublas,
-    random >=1.0 && <1.1,
-    base >=4.5 && <4.8
+    random >=1.0 && <1.2,
+    base
+
+Benchmark newton-inverse
+  Type:             exitcode-stdio-1.0
+  Main-Is:          NewtonInverse.hs
+  Hs-Source-Dirs:   benchmark
+  Default-Language: Haskell98
+  GHC-Options:      -Wall -threaded
+  GHC-Prof-Options: -fprof-auto
+  Build-Depends:
+    accelerate-cublas,
+    accelerate-arithmetic,
+    accelerate-utility,
+    accelerate-cuda,
+    cuda >=0.5 && <0.8,
+    accelerate,
+    pooled-io >=0.0 && <0.1,
+    timeit >=1.0 && <1.1,
+    hmatrix >=0.15.2 && <0.16,
+    random >=1.0.1 && <1.2,
+    utility-ht,
+    base
diff --git a/benchmark/NewtonInverse.hs b/benchmark/NewtonInverse.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/NewtonInverse.hs
@@ -0,0 +1,225 @@
+{-# LANGUAGE FlexibleContexts #-}
+module Main where
+
+import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Batched
+
+import qualified Data.Array.Accelerate.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) =>
+   Batched.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) =>
+   Batched.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
+
+{-
+das wird durch mehrere Kerne gar nicht schneller,
+vielleicht st�rt das Serialisieren der Zufallszahlen?
+-}
+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 <- Batched.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 <- Batched.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) $ 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/example/Main.hs b/example/Main.hs
--- a/example/Main.hs
+++ b/example/Main.hs
@@ -2,7 +2,7 @@
 
 import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Batched
 
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg
 
 import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp
 import Data.Array.Accelerate.Utility.Lift.Exp (expr)
@@ -11,8 +11,6 @@
 import qualified Data.Array.Accelerate.CUDA as AC
 import qualified Data.Array.Accelerate as A
 
-import qualified Foreign.CUDA.Cublas as Cublas
-
 import System.Random (randomRs, mkStdGen)
 
 
@@ -25,7 +23,7 @@
    in  (A.generate (A.constant $ Z :. numMats :. 3 :. 4) f,
         A.generate (A.constant $ Z :. numMats :. 4 :. 2) f)
 
-mainMul :: Cublas.Handle -> IO ()
+mainMul :: Batched.Handle -> IO ()
 mainMul handle = do
    print factorMatrices
    print $ AC.run $
@@ -78,7 +76,7 @@
           let ix = A.index2 i j
           in  k==*0 ? (x!ix, y!ix)
 
-mainLU :: Cublas.Handle -> IO ()
+mainLU :: Batched.Handle -> IO ()
 mainLU handle = do
    print luMatrices
    let mat =
@@ -90,7 +88,7 @@
       Batched.mul handle 1 mat $
       A.replicate (A.lift $ Z :. (2::Int) :. A.All :. A.All) rhsMatrix
 
-mainInv :: Cublas.Handle -> IO ()
+mainInv :: Batched.Handle -> IO ()
 mainInv handle = do
    let dim = 4
        mat =
@@ -103,7 +101,7 @@
 
 main :: IO ()
 main = do
-   handle <- Cublas.create
+   handle <- Batched.create
    mainMul handle
    mainLU handle
    mainInv handle
diff --git a/src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs b/src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs
--- a/src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs
+++ b/src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs
@@ -1,8 +1,6 @@
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ConstraintKinds #-}
-{-# LANGUAGE FlexibleContexts #-}
 module Data.Array.Accelerate.CUBLAS.Level2.Batched (
+   Cublas.Handle,
+   Cublas.create,
    Level3.Element,
    mul,
    mac,
@@ -11,7 +9,7 @@
 import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Level3
 import Data.Array.Accelerate.CUBLAS.Level3.Batched (Element)
 
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg
 
 import qualified Data.Array.Accelerate as A
 import Data.Array.Accelerate (Exp)
diff --git a/src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs b/src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs
--- a/src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs
+++ b/src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs
@@ -1,8 +1,7 @@
 {-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ConstraintKinds #-}
-{-# LANGUAGE FlexibleContexts #-}
 module Data.Array.Accelerate.CUBLAS.Level3.Batched (
+   Cublas.Handle,
+   Cublas.create,
    Element,
    mul,
    mac,
@@ -18,7 +17,7 @@
 import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign as Foreign
 import Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign (Element, Vector)
 
-import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg
+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg
 
 import qualified Data.Array.Accelerate.Utility.Sliced1 as Sliced1
 import qualified Data.Array.Accelerate.Utility.Sliced as Sliced
@@ -54,7 +53,7 @@
    ALinAlg.Matrix ix a
 mul handle alpha a b =
    A.foreignAcc
-      (AF.CUDAForeignAcc "mul" $ uncurry3 $ Foreign.mul handle)
+      (AF.CUDAForeignAcc "mul" $ const $ uncurry3 $ Foreign.mul handle)
       (Acc.modify (expr,acc,acc) $ \(alpha0, a0, b0) ->
          A.map (alpha0 *) $
          ALinAlg.multiplyMatrixMatrix a0 b0)
@@ -69,7 +68,7 @@
    ALinAlg.Matrix ix a
 mac handle alpha a b beta c =
    A.foreignAcc
-      (AF.CUDAForeignAcc "mac" $
+      (AF.CUDAForeignAcc "mac" $ const $
        \((alpha0, a0, b0), (beta0, c0)) ->
           Foreign.mac handle alpha0 a0 b0 beta0 c0)
       (Acc.modify ((expr,acc,acc),(expr,acc)) $
@@ -262,5 +261,5 @@
    String -> (acc -> AF.CIO res) -> A.Acc acc -> A.Acc res
 cudaAcc name f =
    A.foreignAcc
-      (AF.CUDAForeignAcc name f)
+      (AF.CUDAForeignAcc name $ const f)
       (error $ name ++ ": requires CUDA backend")
diff --git a/src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs b/src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs
--- a/src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs
+++ b/src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs
@@ -1,6 +1,5 @@
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
 module Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign where
 
@@ -191,12 +190,17 @@
    return (b, info)
 
 
-type Element a =
-        (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),
-         Fractional (StorableOf a),
-         Cublas.Cublas (StorableOf a),
-         Storable (StorableOf a),
-         Real a)
+class
+   (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),
+    Fractional (StorableOf a),
+    Cublas.Cublas (StorableOf a),
+    Storable (StorableOf a),
+    Real a) =>
+       Element a where
+
+instance Element Float where
+instance Element Double where
+
 
 type family StorableOf float
 type instance StorableOf Float = CFloat