diff --git a/accelerate-arithmetic.cabal b/accelerate-arithmetic.cabal
--- a/accelerate-arithmetic.cabal
+++ b/accelerate-arithmetic.cabal
@@ -1,5 +1,5 @@
 Name:             accelerate-arithmetic
-Version:          0.0
+Version:          0.0.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -17,7 +17,7 @@
 Build-Type:       Simple
 
 Source-Repository this
-  Tag:         0.0
+  Tag:         0.0.1
   Type:        darcs
   Location:    http://code.haskell.org/~thielema/accelerate-arithmetic/
 
@@ -27,7 +27,7 @@
 
 Library
   Build-Depends:
-    accelerate-utility >=0.0 && <0.1,
+    accelerate-utility >=0.1 && <0.2,
     accelerate >=0.15 && <0.16,
     utility-ht >=0.0.8 && <0.1,
     QuickCheck >=2.4 && <2.8,
diff --git a/benchmark/CUBLASBatched.hs b/benchmark/CUBLASBatched.hs
--- a/benchmark/CUBLASBatched.hs
+++ b/benchmark/CUBLASBatched.hs
@@ -6,6 +6,9 @@
 
 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
@@ -26,7 +29,7 @@
 type Scalar ix = Array ix
 
 mul ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
+   (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 ->
@@ -34,7 +37,9 @@
 mul handle alpha a b =
    A.foreignAcc
       (AF.CUDAForeignAcc "mul" $ uncurry3 $ mulPlain handle)
-      (error "Requires CUDA backend")
+      (Acc.modify (expr,acc,acc) $ \(alpha0, a0, b0) ->
+         A.map (alpha0 *) $
+         ALinAlg.multiplyMatrixMatrix a0 b0)
    $
    A.lift (A.unit alpha, a, b)
 
@@ -54,7 +59,7 @@
    (pas, lda) <- arrayPtrs a
    (pbs, ldb) <- arrayPtrs b
    (pcs, ldc) <- arrayPtrs c
-   AF.liftIO $ do
+   AF.liftIO $
       Cublas.gemmBatched handle Cublas.N Cublas.N m n k
          (storableFromScalar alpha)
          pbs ldb
@@ -65,7 +70,7 @@
    return c
 
 mac ::
-   (A.Shape ix, Eq ix, Element a, A.Elt a) =>
+   (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 ->
@@ -73,11 +78,16 @@
 mac handle alpha a b beta c =
    A.foreignAcc
       (AF.CUDAForeignAcc "mac" $
-       \(aalpha, aa, ab, abeta, ac) ->
-          macPlain handle aalpha aa ab abeta ac)
-      (error "Requires CUDA backend")
+       \((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)
+   A.lift ((A.unit alpha, a, b), (A.unit beta, c))
 
 macPlain ::
    (A.Shape ix, Eq ix, Element a, A.Elt a) =>
@@ -100,7 +110,7 @@
    (pas, lda) <- arrayPtrs a
    (pbs, ldb) <- arrayPtrs b
    (pds, ldd) <- arrayPtrs d
-   AF.liftIO $ do
+   AF.liftIO $
       Cublas.gemmBatched handle Cublas.N Cublas.N m n k
          (storableFromScalar alpha)
          pbs ldb
diff --git a/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs b/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
--- a/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
+++ b/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs
@@ -3,6 +3,7 @@
 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
@@ -149,6 +150,23 @@
          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
@@ -210,6 +228,24 @@
       (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) =>
diff --git a/src/Data/Array/Accelerate/Arithmetic/Sparse.hs b/src/Data/Array/Accelerate/Arithmetic/Sparse.hs
--- a/src/Data/Array/Accelerate/Arithmetic/Sparse.hs
+++ b/src/Data/Array/Accelerate/Arithmetic/Sparse.hs
@@ -5,7 +5,7 @@
 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 (atom, )
+import Data.Array.Accelerate.Utility.Lift.Exp (expr, )
 
 import Data.Array.Accelerate.Arithmetic.LinearAlgebra
           (Matrix, Vector, matrixShape, )
@@ -101,7 +101,7 @@
       (_ :. xRows :. _xCols, _ :. _yRows :. yCols) ->
          -- it must be xCols == yRows
          A.zipWith
-            (Exp.modify2 (atom,atom) (atom,atom) $
+            (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)