diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,13 @@
 # Change log for the `patch-image` package
 
+## 0.3.2:
+
+ * Use package `dsp` instead of `hmatrix`
+   for the solution of the linear least squares problem
+   for determining absolute coordinates from image pair displacements.
+   This removes dependency from LAPACK and GSL
+   and makes the code a bit simpler.
+
 ## 0.3.1:
 
  * Speed up computation by moving more stuff to Knead/LLVM.
diff --git a/patch-image.cabal b/patch-image.cabal
--- a/patch-image.cabal
+++ b/patch-image.cabal
@@ -1,5 +1,5 @@
 Name:           patch-image
-Version:        0.3.1
+Version:        0.3.2
 License:        BSD3
 License-File:   LICENSE
 Author:         Henning Thielemann <haskell@henning-thielemann.de>
@@ -41,7 +41,7 @@
   README.md
 
 Source-Repository this
-  Tag:         0.3.1
+  Tag:         0.3.2
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/patch-image/
 
@@ -83,28 +83,32 @@
 
   If flag(llvm)
     Build-Depends:
-      knead >=0.2.2 && <0.3,
-      llvm-extra >=0.7 && <0.8,
+      knead >=0.3 && <0.4,
+      llvm-extra >=0.8 && <0.9,
       llvm-tf >=3.1 && <3.2,
       tfp >=1.0 && <1.1,
       JuicyPixels >=2.0 && <3.3,
-      hmatrix >=0.15 && <0.16,
+      dsp >=0.2.4 && <0.3,
       vector >=0.10 && <0.13,
-      pqueue >=1.2 && <1.4,
+      pqueue >=1.2 && <1.5,
       enumset >=0.0.4 && <0.1,
-      containers >=0.4.2 && <0.6,
+      containers >=0.4.2 && <0.7,
+      semigroups >=0.1 && <1.0,
       fft >=0.1.7 && <0.2,
+      storable-complex >=0.2.2 && <0.3,
       storable-tuple >=0.0.3 && <0.1,
+      bool8 >=0.0 && <0.1,
       carray >=0.1.5 && <0.2,
       array >=0.4 && <0.6,
-      cassava >=0.4.5 && <0.5,
-      unordered-containers >=0.2.5 && <0.2.9,
+      cassava >=0.4.5 && <0.5.2,
+      unordered-containers >=0.2.5 && <0.2.10,
       bytestring >=0.9.2 && <0.11,
       explicit-exception >=0.1.7 && <0.2,
-      Cabal >=1.18 && <2,
-      filepath >=1.3 && <1.4,
+      Cabal >=1.18 && <3,
+      filepath >=1.3 && <1.5,
       non-empty >=0.2 && <0.4,
       utility-ht >=0.0.13 && <0.1,
+      prelude-compat ==0.0.*,
       base >=4 && <5
   Else
     Buildable: False
@@ -126,19 +130,27 @@
 
   If flag(cuda)
     Build-Depends:
-      accelerate-fourier >=0.0 && <0.1,
-      accelerate-arithmetic >=0.1 && <0.2,
-      accelerate-utility >=0.1 && <0.2,
-      accelerate-cufft >=0.0 && <0.1,
-      accelerate-cuda >=0.15 && <0.17,
-      accelerate-io >=0.15 && <0.16,
-      accelerate >=0.15 && <0.16,
+      accelerate-fourier >=1.0 && <1.1,
+      accelerate-arithmetic >=1.0 && <1.1,
+      accelerate-utility >=1.0 && <1.1,
+      accelerate-cufft >=1.0 && <1.1,
+      accelerate-llvm-ptx >=1.1 && <1.2,
+      accelerate-io >=1.0 && <1.1,
+      accelerate >=1.1 && <1.2,
       JuicyPixels >=2.0 && <3.3,
-      hmatrix >=0.15 && <0.16,
+      cassava >=0.4.5 && <0.5.2,
+      dsp >=0.2.4 && <0.3,
       gnuplot >=0.5 && <0.6,
+      containers >=0.4.2 && <0.7,
+      array >=0.4 && <0.6,
       vector >=0.10 && <0.13,
-      Cabal >=1.18 && <2,
-      filepath >=1.3 && <1.4,
+      unordered-containers >=0.2.5 && <0.2.10,
+      bytestring >=0.9.2 && <0.11,
+      enumset >=0.0.4 && <0.1,
+      explicit-exception >=0.1.7 && <0.2,
+      Cabal >=1.18 && <3,
+      filepath >=1.3 && <1.5,
+      non-empty >=0.2 && <0.4,
       utility-ht >=0.0.1 && <0.1,
       base >=4 && <5
   Else
diff --git a/src/Accelerate.hs b/src/Accelerate.hs
--- a/src/Accelerate.hs
+++ b/src/Accelerate.hs
@@ -1,11 +1,11 @@
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
 module Main where
 
 import qualified Option
 import qualified State
 
 import qualified Arithmetic as Arith
-import qualified Complex as Komplex
 import qualified Degree
 import LinearAlgebra (
    absolutePositionsFromPairDisplacements, fixAtLeastOnePosition,
@@ -31,8 +31,8 @@
 import qualified Data.Array.Accelerate.Fourier.Real as FourierReal
 import qualified Data.Array.Accelerate.CUFFT.Single as CUFFT
 import qualified Data.Array.Accelerate.Data.Complex as AComplex
-import qualified Data.Array.Accelerate.CUDA.Foreign as CUDAForeign
-import qualified Data.Array.Accelerate.CUDA as CUDA
+import qualified Data.Array.Accelerate.Data.Bits as ABits
+import qualified Data.Array.Accelerate.LLVM.PTX as CUDA
 import qualified Data.Array.Accelerate.IO as AIO
 import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg
 import qualified Data.Array.Accelerate.Utility.Lift.Run as Run
@@ -46,8 +46,7 @@
 import Data.Array.Accelerate.Utility.Ord (argmaximum)
 import Data.Array.Accelerate
           (Acc, Array, Exp, DIM1, DIM2, DIM3,
-           (:.)((:.)), Z(Z), Any(Any), All(All),
-           (<*), (<=*), (>=*), (==*), (&&*), (||*), (?), (!), )
+           (:.)((:.)), Z(Z), Any(Any), All(All), (?), (!), )
 
 import qualified Graphics.Gnuplot.Advanced as GP
 import qualified Graphics.Gnuplot.LineSpecification as LineSpec
@@ -106,8 +105,7 @@
                      (SV.length dat)
                return $
                   AIO.fromVectors
-                     (Z :. Pic.imageHeight pic :. Pic.imageWidth pic :. 3)
-                     ((), dat)
+                     (Z :. Pic.imageHeight pic :. Pic.imageWidth pic :. 3) dat
             _ -> ioError $ userError "unsupported image type"
 
 writeImage :: Int -> FilePath -> ColorImage8 -> IO ()
@@ -117,7 +115,7 @@
       Pic.Image {
          Pic.imageWidth = width,
          Pic.imageHeight = height,
-         Pic.imageData = snd $ AIO.toVectors arr
+         Pic.imageData = AIO.toVectors arr
       }
 
 writeGrey :: Int -> FilePath -> Array DIM2 Word8 -> IO ()
@@ -127,7 +125,7 @@
       Pic.Image {
          Pic.imageWidth = width,
          Pic.imageHeight = height,
-         Pic.imageData = snd $ AIO.toVectors arr
+         Pic.imageData = AIO.toVectors arr
       }
 
 colorImageExtent :: ColorImage8 -> (Int, Int)
@@ -135,12 +133,12 @@
    case A.arrayShape pic of Z:.height:.width:._chans -> (width, height)
 
 imageFloatFromByte ::
-   (A.Shape sh, A.Elt a, A.IsFloating a) =>
+   (A.Shape sh, A.Floating a, A.FromIntegral Word8 a) =>
    Acc (Array sh Word8) -> Acc (Array sh a)
 imageFloatFromByte = A.map ((/255) . A.fromIntegral)
 
 imageByteFromFloat ::
-   (A.Shape sh, A.Elt a, A.IsFloating a) =>
+   (A.Shape sh, A.RealFloat a) =>
    Acc (Array sh a) -> Acc (Array sh Word8)
 imageByteFromFloat = A.map (fastRound . (255*) . max 0 . min 1)
 
@@ -172,21 +170,30 @@
       arr
 
 
-fastRound ::
-   (A.Elt i, A.IsIntegral i, A.Elt a, A.IsFloating a) => Exp a -> Exp i
+fastRound :: (A.Elt i, A.IsIntegral i, A.RealFloat a) => Exp a -> Exp i
 fastRound x = A.floor (x+0.5)
 
 floatArray :: Acc (Array sh Float) -> Acc (Array sh Float)
 floatArray = id
 
 
-splitFraction :: (A.Elt a, A.IsFloating a) => Exp a -> (Exp Int, Exp a)
+splitFraction ::
+   (A.RealFloat a, A.FromIntegral Int a) => Exp a -> (Exp Int, Exp a)
 splitFraction x =
    let i = A.floor x
    in  (i, x - A.fromIntegral i)
 
 
+target :: CUDA.PTX
+target = unsafePerformIO CUFFT.getBestTarget
 
+cudaRun :: (A.Arrays a) => Acc a -> a
+cudaRun = CUDA.runWith target
+
+cudaRun1 :: (A.Arrays a, A.Arrays b) => (Acc a -> Acc b) -> a -> b
+cudaRun1 = CUDA.run1With target
+
+
 type Channel ix = Array (ix :. Int :. Int)
 type Plane = Channel Z
 
@@ -209,7 +216,7 @@
    in  arr ! A.lift (ix :. yc :. xc)
 
 indexFrac ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) =>
    Acc (Channel ix a) -> Exp ix :. Exp a :. Exp a -> Exp a
 indexFrac arr (ix:.y:.x) =
    let (xi,xf) = splitFraction x
@@ -230,7 +237,7 @@
 
 
 rotateStretchMoveCoords ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.Floating a, A.FromIntegral Int a) =>
    (Exp a, Exp a) ->
    (Exp a, Exp a) ->
    (Exp Int, Exp Int) ->
@@ -242,15 +249,15 @@
           in  A.lift $ trans (A.fromIntegral xdst, A.fromIntegral ydst)
 
 inBox ::
-   (A.Elt a, A.IsNum a, A.IsScalar a) =>
+   (A.Num a, A.Ord a) =>
    (Exp a, Exp a) ->
    (Exp a, Exp a) ->
    Exp Bool
 inBox (width,height) (x,y) =
-   0<=*x &&* x<*width &&* 0<=*y &&* y<*height
+   0 A.<= x  A.&&  x A.< width  A.&&  0 A.<= y  A.&&  y A.< height
 
 validCoords ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a) =>
    (Exp Int, Exp Int) ->
    Acc (Channel Z (a, a)) ->
    Acc (Channel Z Bool)
@@ -269,7 +276,7 @@
 and then moves the picture.
 -}
 rotateStretchMove ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) =>
    (Exp a, Exp a) ->
    (Exp a, Exp a) ->
    ExpDIM2 ix -> Acc (Channel ix a) ->
@@ -289,7 +296,7 @@
 
 
 rotateLeftTop ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) =>
    (Exp a, Exp a) -> Acc (Channel ix a) ->
    ((Acc (A.Scalar a), Acc (A.Scalar a)), Acc (Channel ix a))
 rotateLeftTop rot arr =
@@ -302,7 +309,7 @@
            (chans :. A.ceiling (bottom-top) :. A.ceiling (right-left)) arr)
 
 rotate ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.RealFloat a, A.FromIntegral Int a) =>
    (Exp a, Exp a) ->
    Acc (Channel ix a) -> Acc (Channel ix a)
 rotate rot arr = snd $ rotateLeftTop rot arr
@@ -321,7 +328,7 @@
    Degree Float -> ColorImage8 -> (ColorImage8, Array DIM1 Float)
 rotateHistogram =
    let rot =
-          Run.with CUDA.run1 $ \orient arr ->
+          Run.with cudaRun1 $ \orient arr ->
              let rotated =
                     rotate orient $
                     separateChannels $ imageFloatFromByte arr
@@ -357,9 +364,7 @@
 
 
 
-differentiate ::
-   (A.Elt a, A.IsNum a) =>
-   Acc (Array DIM1 a) -> Acc (Array DIM1 a)
+differentiate :: (A.Num a) => Acc (Array DIM1 a) -> Acc (Array DIM1 a)
 differentiate arr =
    let size = A.unindex1 $ A.shape arr
    in  A.generate (A.index1 (size-1)) $ \i ->
@@ -368,7 +373,7 @@
 scoreRotation :: Degree Float -> ColorImage8 -> Float
 scoreRotation =
    let rot =
-          Run.with CUDA.run1 $ \orient arr ->
+          Run.with cudaRun1 $ \orient arr ->
              A.sum $ A.map (^(2::Int)) $ differentiate $ rowHistogram $
              rotate orient $ separateChannels $ imageFloatFromByte arr
    in  \angle arr -> Acc.the $ rot (Degree.cis angle) arr
@@ -378,18 +383,16 @@
    Key.maximum (flip scoreRotation pic) angles
 
 
-magnitudeSqr :: (A.Elt a, A.IsNum a) => Exp (Complex a) -> Exp a
+magnitudeSqr :: (A.Num a) => Exp (Complex a) -> Exp a
 magnitudeSqr =
    Exp.modify (expr:+expr) $ \(r:+i) -> r*r+i*i
 
 fourierTransformationRun :: ColorImage8 -> IO (Array DIM2 Word8)
 fourierTransformationRun pic = do
    let (shape@(Z:.height:.width):._) = A.arrayShape pic
-   plan <-
-      CUDAForeign.inDefaultContext $
-      CUFFT.plan2D CUFFT.forwardReal shape
+   plan <- CUFFT.plan2D target CUFFT.forwardReal shape
    let trans =
-          Run.with CUDA.run1 $ \arr ->
+          Run.with cudaRun1 $ \arr ->
              imageByteFromFloat $
              A.map (1e-9*) $
              A.zipWith (*)
@@ -423,7 +426,7 @@
 thus I leave it as it is.
 -}
 scoreSlopes ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    (Exp Int, Exp Int) ->
    Acc (Channel Z (Complex a)) -> Acc (Array DIM1 a)
 scoreSlopes (minX, maxX) arr =
@@ -453,16 +456,14 @@
    (Degree Float, Degree Float) -> ColorImage8 -> IO (Degree Float)
 radonAngle (minAngle,maxAngle) pic = do
    let (shape@(Z :. height :. _width):._) = A.arrayShape pic
-   plan <-
-      CUDAForeign.inDefaultContext $
-      CUFFT.plan2D CUFFT.forwardReal shape
+   plan <- CUFFT.plan2D target CUFFT.forwardReal shape
    let height2 = fromIntegral (div height 2)
    let slope w = tan (Degree.toRadian w) * height2
    let minX = floor $ slope minAngle
    let maxX = ceiling $ slope maxAngle
    let angle s = Degree.fromRadian $ atan (s/height2)
    let trans =
-          Run.with CUDA.run1 $ \arr ->
+          Run.with cudaRun1 $ \arr ->
              A.map A.snd $ argmaximum $
              Arrange.mapWithIndex (\ix s -> A.lift (s, A.unindex1 ix)) $
              scoreSlopes (A.constant minX, A.constant maxX) $
@@ -475,7 +476,7 @@
 rotateManifest :: Degree Float -> ColorImage8 -> Array DIM3 Float
 rotateManifest =
    let rot =
-          Run.with CUDA.run1 $ \orient arr ->
+          Run.with cudaRun1 $ \orient arr ->
              rotate orient $ separateChannels $ imageFloatFromByte arr
    in  \angle arr -> rot (Degree.cis angle) arr
 
@@ -484,7 +485,7 @@
    Int -> (Degree Float, ColorImage8) -> ((Float,Float), Plane Float)
 prepareOverlapMatching =
    let rot =
-          Run.with CUDA.run1 $ \radius orient arr ->
+          Run.with cudaRun1 $ \radius orient arr ->
              rotateLeftTop orient $
              (if True
                 then highpass radius
@@ -497,7 +498,7 @@
 
 ceilingPow2 :: Exp Int -> Exp Int
 ceilingPow2 n =
-   A.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double)
+   ABits.setBit 0 $ A.ceiling $ logBase 2 (fromIntegral n :: Exp Double)
 
 pad ::
    (A.Elt a) =>
@@ -506,12 +507,12 @@
    let (height, width) = A.unlift $ A.unindex2 $ A.shape arr
    in  A.generate sh $ \p ->
           let (y, x) = A.unlift $ A.unindex2 p
-          in  (y<*height &&* x<*width)
+          in  (y A.< height  A.&&  x A.< width)
               ?
               (arr ! A.index2 y x, a)
 
 mulConj ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp (Complex a) -> Exp (Complex a) -> Exp (Complex a)
 mulConj x y = x * AComplex.conjugate y
 
@@ -520,15 +521,14 @@
    (A.Elt e, CUFFT.Real e) =>
    CUFFT.Mode DIM2 e a b -> DIM2 -> CUFFT.Transform DIM2 a b
 fft2DGen mode sh =
-   CUFFT.transform $ unsafePerformIO $
-      CUDAForeign.inDefaultContext $ CUFFT.plan2D mode sh
+   CUFFT.transform $ unsafePerformIO $ CUFFT.plan2D target mode sh
 
 fft2DPlain ::
    (A.Elt e, CUFFT.Real e, A.Elt a, A.Elt b) =>
    CUFFT.Mode DIM2 e a b ->
    Channel Z a -> Acc (Channel Z b)
 fft2DPlain mode arr =
-   A.use $ CUDA.run1 (fft2DGen mode $ A.arrayShape arr) arr
+   A.use $ cudaRun1 (fft2DGen mode $ A.arrayShape arr) arr
 
 fft2D ::
    (A.Elt e, CUFFT.Real e, A.Elt a, A.Elt b) =>
@@ -547,13 +547,14 @@
        width  = ceilingPow2 $ widthx  + widthy
        height = ceilingPow2 $ heightx + heighty
        sh = A.index2 height width
-       forward z = fft2DPlain CUFFT.forwardReal $ CUDA.run $ pad 0 sh z
-   in  fft2DPlain CUFFT.inverseReal $ CUDA.run $
+       forward z = fft2DPlain CUFFT.forwardReal $ cudaRun $ pad 0 sh z
+   in  fft2DPlain CUFFT.inverseReal $ cudaRun $
        A.zipWith mulConj (forward x) (forward y)
 
 
 removeDCOffset ::
-   (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z a)
+   (A.Floating a, A.FromIntegral Int a) =>
+   Acc (Channel Z a) -> Acc (Channel Z a)
 removeDCOffset arr =
    let sh = A.shape arr
        (_z :. height :. width) = unliftDim2 sh
@@ -567,20 +568,19 @@
 because we already padded the images with zeros.
 -}
 clearDCCoefficient ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Acc (Array DIM2 (Complex a)) -> Acc (Array DIM2 (Complex a))
 clearDCCoefficient arr =
    A.generate (A.shape arr) $ \p ->
       let (_z:.y:.x) = unliftDim2 p
-      in  x==*0 ||* y==*0 ? (0, arr!p)
+      in  x A.== 0  A.||  y A.== 0 ? (0, arr!p)
 
 
 lowpass, highpass ::
-   (A.Elt a, A.IsFloating a) =>
-   Exp Int -> Acc (Channel Z a) -> Acc (Channel Z a)
+   (A.Floating a) => Exp Int -> Acc (Channel Z a) -> Acc (Channel Z a)
 lowpass count =
    Loop.nest count $
-      A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.Clamp
+      A.stencil (\(a,m,b) -> smooth3 (smooth3 a, smooth3 m, smooth3 b)) A.clamp
 
 highpass count arr =
    A.zipWith (-) arr $ lowpass count arr
@@ -619,7 +619,7 @@
 
 
 wrap :: Exp Int -> Exp Int -> Exp Int -> Exp Int
-wrap size split c = c<*split ? (c, c-size)
+wrap size split c = c A.< split ? (c, c-size)
 
 displacementMap ::
    Exp Int -> Exp Int -> Exp DIM2 -> Acc (Channel Z (Int, Int))
@@ -643,7 +643,7 @@
 that are actually digitalization artifacts.
 -}
 minimumOverlapScores ::
-   (A.Elt a, A.IsFloating a, A.IsScalar a) =>
+   (A.Floating a, A.IsScalar a) =>
    ((Exp Int, Exp Int) -> Exp a -> Exp a) ->
    Exp Int -> (Exp Int, Exp Int) -> (Exp Int, Exp Int) ->
    Acc (Channel Z (a, (Int, Int))) ->
@@ -653,7 +653,7 @@
        (Exp.modify (expr,(expr,expr)) $ \(v, dp@(dx,dy)) ->
           let clipWidth  = min widtha  (widthb  + dx) - max 0 dx
               clipHeight = min heighta (heightb + dy) - max 0 dy
-          in  ((clipWidth >=* minOverlap  &&*  clipHeight >=* minOverlap)
+          in  ((clipWidth A.>= minOverlap   A.&&   clipHeight A.>= minOverlap)
                ?
                (weight (clipWidth, clipHeight) v, 0),
                dp))
@@ -694,7 +694,7 @@
 allOverlapsRun ::
    DIM2 -> Float -> Plane Float -> Plane Float -> Plane Word8
 allOverlapsRun padExtent =
-   Run.with CUDA.run1 $ \minOverlap picA picB ->
+   Run.with cudaRun1 $ \minOverlap picA picB ->
       imageByteFromFloat $
       -- A.map (2*) $
       A.map (0.0001*) $
@@ -704,13 +704,13 @@
    DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int))
 optimalOverlap padExtent =
    let run =
-          Run.with CUDA.run1 $ \minimumOverlap a b ->
+          Run.with cudaRun1 $ \minimumOverlap a b ->
           argmaximum $ allOverlaps padExtent minimumOverlap a b
    in  \overlap a b -> Acc.the $ run overlap a b
 
 
 shrink ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.Floating a, A.FromIntegral Int a) =>
    GenDIM2 (Exp Int) -> Acc (Channel ix a) -> Acc (Channel ix a)
 shrink (_:.yk:.xk) arr =
    let (shape:.height:.width) = unliftDim2 $ A.shape arr
@@ -752,7 +752,7 @@
    DIM2 -> Float -> Plane Float -> Plane Float -> (Float, (Int, Int))
 optimalOverlapBig padExtent =
    let run =
-          Run.with CUDA.run1 $ \minimumOverlap a b ->
+          Run.with cudaRun1 $ \minimumOverlap a b ->
              let factors@(_z:.yk:.xk) =
                     shrinkFactors (A.floor, A.fromIntegral) padExtent
                        minimumOverlap
@@ -804,7 +804,7 @@
 optimalOverlapBigFine padExtent@(Z:.heightPad:.widthPad) =
    let overlaps = allOverlaps padExtent
        run =
-          Run.with CUDA.run1 $ \minimumOverlap a b ->
+          Run.with cudaRun1 $ \minimumOverlap a b ->
              let shapeA = A.unlift $ A.shape a
                  shapeB = A.unlift $ A.shape b
                  factors@(_z:.yk:.xk) =
@@ -848,20 +848,20 @@
    [(Float, (Int, Int), (Int, Int))]
 optimalOverlapBigMulti padExtent (Z:.heightStamp:.widthStamp) numCorrs =
    let overlapShrunk =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \minimumOverlap factors a b ->
              argmaximum $
              allOverlaps padExtent minimumOverlap
                 (shrink factors a) (shrink factors b)
        diffShrunk =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \shrunkd factors a b ->
              overlapDifference shrunkd
                 (shrink factors a) (shrink factors b)
 
        allOverlapsFine = allOverlaps (Z :. 2*heightStamp :. 2*widthStamp)
        overlapFine =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \minimumOverlap a b anchorA@(leftA, topA) anchorB@(leftB, topB)
                 extent@(width,height) ->
              let addCoarsePos =
@@ -922,7 +922,7 @@
 
 
 overlapDifference ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.Floating a, A.FromIntegral Int a) =>
    (Exp Int, Exp Int) ->
    Acc (Channel ix a) -> Acc (Channel ix a) -> Acc (A.Scalar a)
 overlapDifference (dx,dy) a b =
@@ -947,7 +947,7 @@
    (Int, Int) ->
    Plane Float -> Plane Float -> Float
 overlapDifferenceRun =
-   let diff = Run.with CUDA.run1 overlapDifference
+   let diff = Run.with cudaRun1 overlapDifference
    in  \d a b -> Acc.the $ diff d a b
 
 
@@ -983,7 +983,7 @@
    let rotat (rot,pic) =
           rotate rot $ separateChannels $ imageFloatFromByte pic
    in  (\f d (a,b) -> f d (mapFst Degree.cis a, mapFst Degree.cis b)) $
-       Run.with CUDA.run1 $
+       Run.with cudaRun1 $
        \(dx,dy) (a,b) ->
           imageByteFromFloat $ interleaveChannels $
           overlap2 (dx, dy) (rotat a, rotat b)
@@ -994,14 +994,14 @@
    ix :. Int :. Int ->
    (Plane Int, Channel ix Float)
 emptyCountCanvas =
-   Run.with CUDA.run1 $ \sh ->
+   Run.with cudaRun1 $ \sh ->
       let (_ix :. height :. width) = unliftDim2 sh
       in  (A.fill (A.lift $ Z:.height:.width) 0,
            A.fill sh 0)
 
 
 addToCountCanvas ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
+   (A.Slice ix, A.Shape ix, A.Num a, A.FromIntegral Int a) =>
    (Acc (Plane Bool), Acc (Channel ix a)) ->
    (Acc (Plane Int),  Acc (Channel ix a)) ->
    (Acc (Plane Int),  Acc (Channel ix a))
@@ -1016,7 +1016,7 @@
    (Plane Int, Channel DIM1 Float) ->
    (Plane Int, Channel DIM1 Float)
 updateCountCanvas =
-   Run.with CUDA.run1 $
+   Run.with cudaRun1 $
    \(rot, mov, pic) (count,canvas) ->
       addToCountCanvas
          (rotateStretchMove rot mov (unliftDim2 $ A.shape canvas) $
@@ -1025,7 +1025,7 @@
 
 finalizeCountCanvas :: (Plane Int, Channel DIM1 Float) -> ColorImage8
 finalizeCountCanvas =
-   Run.with CUDA.run1 $
+   Run.with cudaRun1 $
    \(count, canvas) ->
       imageByteFromFloat $ interleaveChannels $
       A.zipWith (/) canvas $
@@ -1044,30 +1044,30 @@
    in  A.cond xb (A.lift (True, A.cond yb (f xv yv) xv)) y
 
 maskedMinimum ::
-   (A.Shape ix, A.Elt a, A.IsScalar a) =>
+   (A.Shape ix, A.Ord a) =>
    LinAlg.Vector ix (Bool, a) ->
    LinAlg.Scalar ix (Bool, a)
 maskedMinimum = A.fold1 (maybePlus min)
 
 maskedMaximum ::
-   (A.Shape ix, A.Elt a, A.IsScalar a) =>
+   (A.Shape ix, A.Ord a) =>
    LinAlg.Vector ix (Bool, a) ->
    LinAlg.Scalar ix (Bool, a)
 maskedMaximum = A.fold1 (maybePlus max)
 
 
 project ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Point2 (Exp a) ->
    (Point2 (Exp a), Point2 (Exp a)) ->
    (Exp Bool, Point2 (Exp a))
 project x ab =
    let (r, y) = projectPerp x ab
-   in  (0<=*r &&* r<=*1, y)
+   in  (0 A.<= r  A.&&  r A.<= 1, y)
 
 
 distanceMapEdges ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp DIM2 -> Acc (Array DIM1 ((a,a),(a,a))) -> Acc (Channel Z a)
 distanceMapEdges sh edges =
    A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $
@@ -1081,13 +1081,13 @@
 distanceMapEdgesRun ::
    DIM2 -> Array DIM1 ((Float,Float),(Float,Float)) -> Plane Word8
 distanceMapEdgesRun =
-   Run.with CUDA.run1 $ \sh ->
+   Run.with cudaRun1 $ \sh ->
       imageByteFromFloat . A.map (0.01*) . distanceMapEdges sh
 
 type Geometry a = Arith.Geometry Int a
 
 distanceMapBox ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp DIM2 ->
    Exp (Geometry a) ->
    Acc (Channel Z (Bool, (((a,(a,a)), (a,(a,a))), ((a,(a,a)), (a,(a,a))))))
@@ -1138,7 +1138,7 @@
 
 distanceMapBoxRun :: DIM2 -> Geometry Float -> Plane Word8
 distanceMapBoxRun =
-   Run.with CUDA.run1 $ \sh geom ->
+   Run.with cudaRun1 $ \sh geom ->
       scaleDistanceMapGeom geom $
       A.map (Exp.modify (expr,expr) $ \(valid, dist) -> valid ? (dist, 0)) $
       maskedMinimum $
@@ -1161,12 +1161,12 @@
 
 
 containedAnywhere ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a) =>
    Acc (Array DIM1 (Geometry a)) ->
    Acc (Array DIM3 (a,a)) ->
    Acc (Array DIM3 Bool)
 containedAnywhere geoms arr =
-   A.fold1 (||*) $
+   A.fold1 (A.||) $
    breakFusion $
    outerVector
       (Exp.modify2 (expr,expr) ((expr,expr),(expr,expr),(expr,expr)) $
@@ -1177,7 +1177,7 @@
 
 
 distanceMapContained ::
-   (A.IsFloating a, A.Elt a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp DIM2 ->
    Exp (Geometry a) ->
    Acc (Array DIM1 (Geometry a)) ->
@@ -1189,19 +1189,19 @@
        maskedMinimum $
        A.zipWith
           (Exp.modify2 expr (expr,(expr,expr)) $ \c (b,(dist,_)) ->
-             (c&&*b, dist))
+             (c A.&& b, dist))
           contained distMap
 
 distanceMapContainedRun ::
    DIM2 -> Geometry Float -> [Geometry Float] -> Plane Word8
 distanceMapContainedRun =
    let distances =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \sh this -> scaleDistanceMapGeom this . distanceMapContained sh this
    in  \sh this others -> distances sh this $ array1FromList others
 
 scaleDistanceMapGeom ::
-   (A.IsFloating a, A.Elt a, A.Elt b, A.Shape ix) =>
+   (A.RealFloat a, A.FromIntegral Int a, A.Elt b, A.Shape ix) =>
    Exp (Geometry b) -> Acc (Array ix a) -> Acc (Array ix Word8)
 scaleDistanceMapGeom this =
    let scale = (4/) $ A.fromIntegral $ A.uncurry min $ Exp.thd3 this
@@ -1209,14 +1209,13 @@
 
 
 pixelCoordinates ::
-   (A.Elt a, A.IsFloating a) =>
-   Exp DIM2 -> Acc (Channel Z (a,a))
+   (A.RealFloat a, A.FromIntegral Int a) => Exp DIM2 -> Acc (Channel Z (a,a))
 pixelCoordinates sh =
    A.generate sh $ Exp.modify (expr:.expr:.expr) $ \(_z:.y:.x) ->
       (A.fromIntegral x, A.fromIntegral y)
 
 distanceMapPoints ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a) =>
+   (A.Slice ix, A.Shape ix, A.RealFloat a) =>
    Acc (Array ix (a,a)) ->
    Acc (Array DIM1 (a,a)) ->
    Acc (Array ix a)
@@ -1229,12 +1228,13 @@
 distanceMapPointsRun :: DIM2 -> [Point2 Float] -> Plane Word8
 distanceMapPointsRun =
    let distances =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \sh -> scaleDistanceMap . distanceMapPoints (pixelCoordinates sh)
    in  \sh points -> distances sh $ array1FromList points
 
 scaleDistanceMap ::
-   (A.Elt a, A.IsFloating a) => Acc (Channel Z a) -> Acc (Channel Z Word8)
+   (A.RealFloat a, A.FromIntegral Int a) =>
+   Acc (Channel Z a) -> Acc (Channel Z Word8)
 scaleDistanceMap arr =
    let scale =
          case Exp.unlift (expr:.expr:.expr) (A.shape arr) of
@@ -1258,7 +1258,7 @@
 and chose the minimal distance.
 -}
 distanceMap ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp DIM2 ->
    Exp (Geometry a) ->
    Acc (Array DIM1 (Geometry a)) ->
@@ -1277,7 +1277,7 @@
    Plane Word8
 distanceMapRun =
    let distances =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \sh this others -> scaleDistanceMap . distanceMap sh this others
    in  \sh this others points ->
           distances sh this
@@ -1286,7 +1286,7 @@
 
 
 distanceMapGamma ::
-   (A.Elt a, A.IsFloating a) =>
+   (A.RealFloat a, A.FromIntegral Int a) =>
    Exp a ->
    Exp DIM2 ->
    Exp (Geometry a) ->
@@ -1302,14 +1302,14 @@
    ix :. Int :. Int ->
    (Plane Float, Channel ix Float)
 emptyWeightedCanvas =
-   Run.with CUDA.run1 $ \sh ->
+   Run.with cudaRun1 $ \sh ->
       let (_ix :. height :. width) = unliftDim2 sh
       in  (A.fill (A.lift $ Z:.height:.width) 0,
            A.fill sh 0)
 
 
 addToWeightedCanvas ::
-   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>
+   (A.Slice ix, A.Shape ix, A.Num a) =>
    (Acc (Channel Z a), Acc (Channel ix a)) ->
    (Acc (Channel Z a), Acc (Channel ix a)) ->
    (Acc (Channel Z a), Acc (Channel ix a))
@@ -1328,7 +1328,7 @@
    (Plane Float, Channel DIM1 Float)
 updateWeightedCanvasMerged =
    let update =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \this others points pic (weightSum,canvas) ->
              let (rot, mov, _) =
                     Exp.unlift ((expr,expr), (expr,expr), expr) this
@@ -1352,9 +1352,9 @@
    (Plane Float, Channel DIM1 Float) ->
    (Plane Float, Channel DIM1 Float)
 updateWeightedCanvas =
-   let distances = Run.with CUDA.run1 distanceMapGamma
+   let distances = Run.with cudaRun1 distanceMapGamma
        update =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \this pic dist (weightSum,canvas) ->
              let (rot, mov, _) =
                     Exp.unlift ((expr,expr), (expr,expr), expr) this
@@ -1379,10 +1379,10 @@
    (Plane Float, Channel DIM1 Float) ->
    (Plane Float, Channel DIM1 Float)
 updateWeightedCanvasSplit =
-   let update = Run.with CUDA.run1 addToWeightedCanvas
-       distances = Run.with CUDA.run1 distanceMap
+   let update = Run.with cudaRun1 addToWeightedCanvas
+       distances = Run.with cudaRun1 distanceMap
        rotated =
-          Run.with CUDA.run1 $
+          Run.with cudaRun1 $
           \sh rot mov pic ->
              snd $ rotateStretchMove rot mov (unliftDim2 sh) $
              separateChannels $ imageFloatFromByte pic
@@ -1398,7 +1398,7 @@
 finalizeWeightedCanvas ::
    (Plane Float, Channel DIM1 Float) -> ColorImage8
 finalizeWeightedCanvas =
-   Run.with CUDA.run1 $
+   Run.with cudaRun1 $
    \(weightSum, canvas) ->
       imageByteFromFloat $ interleaveChannels $
       A.zipWith (/) canvas $
@@ -1505,8 +1505,7 @@
    let (errdx,errdy) =
           mapPair (maximum0, maximum0) $ unzip $
           zipWith
-             (\(dpx,dpy) (dx,dy) ->
-                (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy))
+             (\(dpx,dpy) (dx,dy) -> (abs $ dpx - dx, abs $ dpy - dy))
              dps (map snd overlaps)
 
    info $
@@ -1516,10 +1515,7 @@
       ++
       printf "maximum vertical error: %f\n" errdy
 
-   return
-      (map picPath pics,
-       map picColored pics,
-       map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss)
+   return (map picPath pics, map picColored pics, map (flip (,) 1) poss)
 
 
 processOverlapRotate ::
@@ -1610,13 +1606,7 @@
                dpx dpy xa ya xb yb)
          dps overlaps
 
-   return
-      (map picPath pics,
-       map picColored pics,
-       map
-         (mapPair
-            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))
-         posRots)
+   return (map picPath pics, map picColored pics, posRots)
 
 
 processRotation ::
@@ -1659,14 +1649,14 @@
       let pic0 : pic1 : _ = map snd rotated
           size = (Z:.512:.1024 :: DIM2)
       writeGrey (Option.quality opt) "/tmp/padded.jpeg" $
-         CUDA.run1 (imageByteFromFloat . pad 0 (A.lift size)) pic0
+         cudaRun1 (imageByteFromFloat . pad 0 (A.lift size)) pic0
       writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" $
-         CUDA.run $ imageByteFromFloat $ A.map AComplex.magnitude $
+         cudaRun $ imageByteFromFloat $ A.map AComplex.magnitude $
          fft2DPlain CUFFT.forwardReal $
-         CUDA.run1 (pad 0 (A.lift size)) $
+         cudaRun1 (pad 0 (A.lift size)) $
          pic0
       writeGrey (Option.quality opt) "/tmp/convolution.jpeg" $
-         CUDA.run $ imageByteFromFloat $ A.map (0.000001*) $
+         cudaRun $ imageByteFromFloat $ A.map (0.000001*) $
          correlatePadded size (A.use pic0) (A.use pic1)
 
    return $
diff --git a/src/Degree.hs b/src/Degree.hs
--- a/src/Degree.hs
+++ b/src/Degree.hs
@@ -4,14 +4,18 @@
 
 import qualified Data.Csv as Csv
 import Data.Monoid (Monoid, mempty, mappend)
+import Data.Semigroup (Semigroup, (<>))
 
 
 newtype Degree a = Degree {getDegree :: a}
    deriving (Eq, Show)
 
+instance (Num a) => Semigroup (Degree a) where
+   Degree x <> Degree y = Degree $ x+y
+
 instance (Num a) => Monoid (Degree a) where
    mempty = Degree 0
-   mappend (Degree x) (Degree y) = Degree $ x+y
+   mappend = (<>)
 
 instance Functor Degree where
    fmap f (Degree x) = Degree $ f x
diff --git a/src/Knead.hs b/src/Knead.hs
--- a/src/Knead.hs
+++ b/src/Knead.hs
@@ -27,11 +27,11 @@
 import qualified Data.Array.Knead.Simple.Physical as Phys
 import qualified Data.Array.Knead.Simple.ShapeDependent as ShapeDep
 import qualified Data.Array.Knead.Simple.Symbolic as Symb
-import qualified Data.Array.Knead.Index.Nested.Shape as Shape
+import qualified Data.Array.Knead.Shape.Nested as Shape
 import qualified Data.Array.Knead.Expression as Expr
 import Data.Array.Knead.Simple.Symbolic ((!))
 import Data.Array.Knead.Expression
-         (Exp, (==*), (/=*), (<*), (<=*), (>=*), (&&*))
+         (Exp, (==*), (<*), (<=*), (>=*), (||*), (&&*))
 
 import Data.Array.IArray (amap)
 import Data.Array.CArray (CArray)
@@ -52,6 +52,7 @@
 
 import qualified Data.Vector.Storable as SV
 import Foreign.ForeignPtr (ForeignPtr, castForeignPtr)
+import Foreign.Storable.Complex ()
 
 import qualified System.FilePath as FilePath
 import qualified System.IO as IO
@@ -68,7 +69,7 @@
 import qualified Data.Foldable as Fold
 import qualified Data.List as List
 import qualified Data.Map as Map
-import Data.Function.HT (Id)
+import qualified Data.Bool8 as Bool8
 import Data.Monoid ((<>))
 import Data.Maybe.HT (toMaybe)
 import Data.Maybe (mapMaybe, isJust, isNothing)
@@ -79,9 +80,13 @@
 import Data.Tuple.HT
          (mapPair, mapFst, mapSnd, mapTriple, swap, mapThd3, fst3, uncurry3)
 import Data.Word (Word8, Word32)
+import Data.Bool8 (Bool8)
 
+import Prelude2010
+import Prelude ()
 
 
+
 type SmallSize = Word32
 
 type Plane = Phys.Array Dim2
@@ -411,17 +416,17 @@
    (Exp a, Exp a) ->
    Exp Bool
 inBox (width,height) (x,y) =
-   Expr.liftM2 MultiValue.and (inRange width x) (inRange height y)
+   inRange width x &&* inRange height y
 
 validCoords ::
    (MultiValue.NativeFloating a ar,
     MultiValue.Field a, MultiValue.Real a,
     MultiValue.RationalConstant a) =>
    (Exp Size, Exp Size) ->
-   SymbPlane (a, a) -> SymbPlane MaskBool
+   SymbPlane (a, a) -> SymbPlane Bool8
 validCoords (width,height) =
    Symb.map $ Expr.modify (atom,atom) $ \(x,y) ->
-      maskFromBool $ inBox (width,height) (fastRound x, fastRound y)
+      Expr.bool8FromP $ inBox (width,height) (fastRound x, fastRound y)
 
 {- |
 @rotateStretchMove rot mov@
@@ -438,7 +443,7 @@
    Exp (a, a) ->
    Exp Dim2 ->
    SymbPlane v ->
-   SymbPlane (MaskBool, v)
+   SymbPlane (Bool8, v)
 rotateStretchMove vec rot mov sh img =
    let coords = rotateStretchMoveCoords rot mov sh
        (Vec2 heightSrc widthSrc) = Expr.decompose atomDim2 $ Symb.shape img
@@ -501,16 +506,8 @@
    (Symb.C array, MultiValue.Additive a) => array Dim1 a -> array Dim1 a
 differentiate xs = Symb.zipWith Expr.sub (tailArr xs) xs
 
-the :: Symb.Array () a -> Exp a
-the = Symb.the
-
-fold1All ::
-   (Shape.C sh, MultiValue.C a) =>
-   (Exp a -> Exp a -> Exp a) -> Symb.Array sh a -> Exp a
-fold1All f = the . Symb.fold1All f
-
 scoreHistogram :: (MultiValue.PseudoRing a) => Symb.Array Dim1 a -> Exp a
-scoreHistogram = fold1All Expr.add . Symb.map Expr.sqr . differentiate
+scoreHistogram = Symb.fold1All Expr.add . Symb.map Expr.sqr . differentiate
 
 
 runScoreRotation :: IO (Degree Float -> ColorImage8 -> IO Float)
@@ -558,7 +555,7 @@
 highpassMulti :: IO (Int -> Plane Float -> IO (Plane Float))
 highpassMulti = do
    lp <- lowpassMulti
-   sub <- RenderP.run $ Symb.zipWith Expr.sub . fixArray
+   sub <- RenderP.run $ Symb.zipWith Expr.sub . Symb.fix
    return $ \n img -> sub img =<< lp n img
 
 
@@ -623,13 +620,10 @@
 liftCArray f a = arrayKneadFromC . f <$> arrayCFromKnead a
 
 
-fixArray :: Id (Symb.Array sh a)
-fixArray = id
-
 prepareOverlapMatching ::
    IO (Int -> (Degree Float, ColorImage8) -> IO ((Float, Float), Plane Float))
 prepareOverlapMatching = do
-   bright <- RenderP.run $ brightnessPlane . colorImageFloatFromByte . fixArray
+   bright <- RenderP.run $ brightnessPlane . colorImageFloatFromByte . Symb.fix
    hp <- highpassMulti
    rotat <- RenderP.run $ rotate Arith.vecScalar
    return $ \radius (angle, img) ->
@@ -737,7 +731,7 @@
    (Shape.C sh,
     MultiValue.Comparison a, MultiValue.Select a, MultiValue.Select b) =>
    Symb.Array sh (a, b) -> Exp (a, b)
-argmaximum = fold1All argmax
+argmaximum = Symb.fold1All argmax
 
 optimalOverlap ::
    Dim2 -> IO (Float -> Plane Float -> Plane Float -> IO (Float, (Size, Size)))
@@ -961,7 +955,7 @@
        extentOverlap = (widthOverlap,heightOverlap)
    in  Expr.sqrt $
        (/(fromInt widthOverlap * fromInt heightOverlap)) $
-       fold1All (+) $
+       Symb.fold1All (+) $
        Symb.map Expr.sqr $
        Symb.zipWith (-)
           (clip (leftOverlap,topOverlap) extentOverlap a)
@@ -1020,31 +1014,20 @@
    RenderP.run $ \sh -> Symb.fill sh (Expr.zip 0 $ Expr.zip3 0 0 0)
 
 
-type MaskBool = Word8
-
-maskFromBool :: Exp Bool -> Exp MaskBool
-maskFromBool = Expr.liftM $ MultiValue.liftM $ LLVM.zext
-
-boolFromMask :: Exp MaskBool -> Exp Bool
-boolFromMask = (/=* 0)
-
-intFromBool :: Exp MaskBool -> Exp Word32
-intFromBool = Expr.liftM $ MultiValue.liftM $ LLVM.ext
-
 type RotatedImage = ((Float,Float), (Float,Float), ColorImage8)
 
 addToCountCanvas ::
    (MultiValue.PseudoRing a, MultiValue.NativeFloating a ar) =>
    VecExp a v ->
-   SymbPlane (MaskBool, v) ->
+   SymbPlane (Bool8, v) ->
    SymbPlane (Word32, v) ->
    SymbPlane (Word32, v)
 addToCountCanvas vec =
    Symb.zipWith
       (Expr.modify2 (atom,atom) (atom,atom) $ \(mask, pic) (count, canvas) ->
-         (Expr.add (intFromBool mask) count,
+         (Expr.add (Expr.intFromBool8 mask) count,
           Arith.vecAdd vec canvas $
-          Arith.vecScale vec (fromInt $ intFromBool mask) pic))
+          Arith.vecScale vec (Expr.floatFromBool8 mask) pic))
 
 updateCountCanvas ::
    IO (RotatedImage -> Plane (Word32, YUV Float) ->
@@ -1063,14 +1046,14 @@
       Symb.map
          (Expr.modify (atom,atom) $ \(count, pixel) ->
             Arith.vecScale vecYUV (recip $ fromInt count) pixel) .
-      fixArray
+      Symb.fix
 
 
 diffAbs :: (MultiValue.Real a) => Exp a -> Exp a -> Exp a
 diffAbs = Expr.liftM2 $ \x y -> MultiValue.abs =<< MultiValue.sub x y
 
 diffWithCanvas ::
-   IO (RotatedImage -> Plane (YUV Float) -> IO (Plane (MaskBool, Float)))
+   IO (RotatedImage -> Plane (YUV Float) -> IO (Plane (Bool8, Float)))
 diffWithCanvas =
    RenderP.run $ \(rot, mov, pic) avg ->
       Symb.zipWith
@@ -1088,20 +1071,20 @@
          (Expr.modify (atom,atom) $ \(count, pixel) ->
             Arith.vecScale vecYUV (recip $ fromInt count) pixel)
       .
-      fixArray
+      Symb.fix
 
 emptyCanvas :: IO (Dim2 -> IO ColorImage8)
 emptyCanvas = RenderP.run $ \sh -> Symb.fill sh (Expr.zip3 0 0 0)
 
 addMaskedToCanvas ::
    IO (RotatedImage ->
-       Plane MaskBool ->
+       Plane Bool8 ->
        Plane (YUV Word8) ->
        IO (Plane (YUV Word8)))
 addMaskedToCanvas =
    RenderP.run $ \(rot, mov, pic) mask canvas ->
       Symb.zipWith3 Expr.ifThenElse
-         (Symb.map boolFromMask mask)
+         (Symb.map Expr.boolPFrom8 mask)
          (Symb.map (yuvByteFromFloat . Expr.snd) $
           rotateStretchMove vecYUV rot mov (Symb.shape canvas) $
           colorImageFloatFromByte pic)
@@ -1116,7 +1099,8 @@
    RenderP.run $ \(rot, mov, pic) shape weightCanvas ->
       addToWeightedCanvas vecYUV
          (Symb.zipWith
-            (Expr.modify2 atom (atom,atom) $ \s (b,x) -> (fromInt b * s, x))
+            (Expr.modify2 atom (atom,atom) $ \s (b,x) ->
+               (Expr.floatFromBool8 b * s, x))
             shape $
           rotateStretchMove vecYUV rot mov (Symb.shape weightCanvas) $
           colorImageFloatFromByte pic)
@@ -1231,7 +1215,7 @@
    array sh (a,a) ->
    array sh Bool
 containedAnywhere geoms array =
-   Symb.fold1 (Expr.liftM2 MultiValue.or) $
+   Symb.fold1 (||*) $
    outerProduct
       (Expr.modify2 (atom,atom) ((atom,atom),(atom,atom),(atom,atom)) $
        \(xdst,ydst) (rot, mov, extent) ->
@@ -1257,7 +1241,7 @@
        maskedMinimum $
        Symb.zipWith
           (Expr.modify2 atom (atom,(atom,atom)) $ \c (b,(dist,_)) ->
-             (Expr.liftM2 MultiValue.and c b, dist))
+             (c &&* b, dist))
           contained distMap
 
 distanceMapContainedRun ::
@@ -1436,7 +1420,7 @@
       Symb.map
          (Expr.modify (atom,atom) $ \(weightSum, pixel) ->
             Arith.vecScale vecYUV (recip weightSum) pixel) .
-      fixArray
+      Symb.fix
 
 
 data
@@ -1544,8 +1528,7 @@
    let (errdx,errdy) =
           mapPair (maximum0, maximum0) $ unzip $
           zipWith
-             (\(dpx,dpy) (dx,dy) ->
-                (abs $ dpx - realToFrac dx, abs $ dpy - realToFrac dy))
+             (\(dpx,dpy) (dx,dy) -> (abs $ dpx - dx, abs $ dpy - dy))
              dps (map snd overlaps)
 
    info $
@@ -1555,10 +1538,7 @@
       ++
       printf "maximum vertical error: %f\n" errdy
 
-   return
-      (map picPath pics,
-       map picColored pics,
-       map (flip (,) 1) $ map (mapPair (realToFrac, realToFrac)) poss)
+   return (map picPath pics, map picColored pics, map (flip (,) 1) poss)
 
 
 processOverlapRotate ::
@@ -1649,13 +1629,7 @@
                dpx dpy xa ya xb yb)
          dps overlaps
 
-   return
-      (map picPath pics,
-       map picColored pics,
-       map
-         (mapPair
-            (mapPair (realToFrac, realToFrac), Komplex.map realToFrac))
-         posRots)
+   return (map picPath pics, map picColored pics, posRots)
 
 
 processRotation ::
@@ -1695,14 +1669,14 @@
       let pic0 : pic1 : _ = map snd rotated
           size = Vec2 1024 768
       makeByteImage <-
-         RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . fixArray
+         RenderP.run $ \k -> imageByteFromFloat . Symb.map (k*) . Symb.fix
       runPad <- RenderP.run pad
       writeGrey (Option.quality opt) "/tmp/padded.jpeg" =<<
          (makeByteImage 1 =<< runPad 0 size pic0)
       runMagnitude <-
          RenderP.run $
          Symb.map (Expr.modify atomComplex $ \(r:+i) -> Expr.sqrt$ r*r+i*i)
-            . fixArray
+            . Symb.fix
       writeGrey (Option.quality opt) "/tmp/spectrum.jpeg" =<<
          (makeByteImage 0.1 =<< runMagnitude =<<
           liftCArray (FFT.dftRCN [0,1]) =<< runPad 0 size pic0)
@@ -1815,9 +1789,9 @@
       avg <- finalizeCanv sumImg
       diff <- diffWithCanvas
       picDiffs <- mapM (flip diff avg) rotMovPics
-      getSnd <- RenderP.run $ Symb.map Expr.snd . fixArray
+      getSnd <- RenderP.run $ Symb.map Expr.snd . Symb.fix
       lp <- lowpassMulti
-      masks <- map (amap ((0/=) . fst)) <$> mapM arrayCFromKnead picDiffs
+      masks <- map (amap (Bool8.toBool . fst)) <$> mapM arrayCFromKnead picDiffs
       let smoothRadius = Option.shapeSmooth opt
       smoothPicDiffs <-
          mapM (arrayCFromKnead <=< lp smoothRadius <=< getSnd) picDiffs
@@ -1840,7 +1814,7 @@
             foldM
                (\canvas (shape, rotMovPic) ->
                   addMasked rotMovPic
-                     (arrayKneadFromC $ amap (fromIntegral . fromEnum) shape)
+                     (arrayKneadFromC $ amap Bool8.fromBool shape)
                      canvas)
                emptyPlain (zip shapes rotMovPics)
 
@@ -1851,7 +1825,7 @@
                 amap (fromIntegral . fromEnum))
                shapes
          forM_ (Option.outputShape opt) $ \format -> do
-            makeByteImage <- RenderP.run $ imageByteFromFloat . fixArray
+            makeByteImage <- RenderP.run $ imageByteFromFloat . Symb.fix
             forM_ (zip names smoothShapes) $ \(name,shape) ->
                writeGrey (Option.quality opt) (printf format name)
                   =<< makeByteImage shape
diff --git a/src/Knead/Shape.hs b/src/Knead/Shape.hs
--- a/src/Knead/Shape.hs
+++ b/src/Knead/Shape.hs
@@ -2,11 +2,12 @@
 {-# LANGUAGE EmptyDataDecls #-}
 module Knead.Shape where
 
-import qualified Data.Array.Knead.Index.Nested.Shape as Shape
+import qualified Data.Array.Knead.Shape.Nested as Shape
 import qualified Data.Array.Knead.Expression as Expr
 
 import qualified LLVM.Extra.Multi.Value.Memory as MultiMem
 import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Iterator as Iter
 import qualified LLVM.Extra.Arithmetic as A
 import LLVM.Extra.Multi.Value (atom)
 
@@ -159,6 +160,11 @@
       case unzipShape nm of
          Vec2 n m ->
             Shape.loop (\i -> Shape.loop (\j -> code (zipShape i j)) m) n
+   iterator nm =
+      case unzipShape nm of
+         Vec2 n m ->
+            fmap (uncurry zipShape) $
+            Iter.cartesian (Shape.iterator n) (Shape.iterator m)
 
 
 instance (Expr.Compose n) => Expr.Compose (Vec2 tag n) where
diff --git a/src/LinearAlgebra.hs b/src/LinearAlgebra.hs
--- a/src/LinearAlgebra.hs
+++ b/src/LinearAlgebra.hs
@@ -1,19 +1,21 @@
 module LinearAlgebra where
 
-import qualified Data.Packed.Matrix as Matrix
-import qualified Data.Packed.Vector as Vector
-import qualified Data.Packed.ST as PackST
-import qualified Numeric.Container as Container
-import Numeric.Container ((<\>), (<>))
+import qualified Matrix.Sparse as Sparse
+import qualified Matrix.Vector as Vector
+import qualified Matrix.QR.Givens as QR
 
 import Data.Complex (Complex((:+)))
 
 import qualified Data.List.HT as ListHT
 import qualified Data.List as List
+import qualified Data.Array as Array
+import qualified Data.Map as Map
+import Data.Array (Array)
 import Data.Tuple.HT (mapPair, mapSnd)
-import Data.Maybe (isJust, fromMaybe)
+import Data.Maybe (isJust)
+import Data.Map (Map)
 
-import Control.Applicative ((<$>))
+import Control.Applicative ((<$>), (<$))
 
 
 fixAtLeastOne :: a -> [Maybe a] -> [Maybe a]
@@ -47,20 +49,33 @@
 parallel fs = uncurry zip . mapPair fs . unzip
 
 
-sparseMatrix :: Int -> Int -> [((Int, Int), Double)] -> Matrix.Matrix Double
-sparseMatrix numRows numCols xs =
-   PackST.runSTMatrix $ do
-      mat <- PackST.newMatrix 0 numRows numCols
-      mapM_ (\((r,c), x) -> PackST.writeMatrix mat r c x) xs
-      return mat
+type Matrix = Sparse.Matrix Int Int
+type Vector = Array Int
 
+sparseMatrix :: Int -> Int -> [((Int, Int), Float)] -> Matrix Float
+sparseMatrix numRows numCols =
+   Sparse.fromMap ((0,0), (numRows-1, numCols-1)) . Map.fromList
+
+-- every column must be non-empty
+sparseToColumns :: Matrix a -> [Map Int a]
+sparseToColumns = Map.elems . Sparse.toColumns
+
+sparseFromColumns :: (Int,Int) -> [Map Int a] -> Matrix a
+sparseFromColumns (m0,m1) cols =
+   Sparse.fromColumns ((m0,0), (m1, length cols-1)) $
+   Map.fromList $ zip [0..] cols
+
+addSparseColumn :: (Num a) => Vector a -> Map Int a -> Vector a
+addSparseColumn v col = Array.accum (+) v $ Map.toList col
+
+
 elm :: Int -> Int -> a -> ((Int, Int), a)
 elm row col x = ((row, col), x)
 
 
 absolutePositionsFromPairDisplacements ::
    [(Maybe Float, Maybe Float)] -> [((Int, Int), (Float, Float))] ->
-   ([(Double,Double)], [(Double,Double)])
+   ([(Float,Float)], [(Float,Float)])
 absolutePositionsFromPairDisplacements mxys displacements =
    let numPics = length mxys
        (mxs, mys) = unzip mxys
@@ -68,18 +83,14 @@
        matrix =
           sparseMatrix (length is) numPics $ concat $
           zipWith (\k (ia,ib) -> [elm k ia (-1), elm k ib 1]) [0..] is
-       solve ms ds =
-          leastSquaresSelected matrix
-             (map (fmap realToFrac) ms)
-             (Vector.fromList (map realToFrac ds))
+       solve ms ds = leastSquaresSelected matrix ms (Vector.fromList ds)
        (pxs, achievedDxs) = solve mxs dxs
        (pys, achievedDys) = solve mys dys
    in  (zip pxs pys, zip achievedDxs achievedDys)
 
 
 leastSquaresSelected ::
-   Matrix.Matrix Double -> [Maybe Double] -> Vector.Vector Double ->
-   ([Double], [Double])
+   Matrix Float -> [Maybe Float] -> Vector Float -> ([Float], [Float])
 leastSquaresSelected m mas rhs0 =
    let (lhsCols,rhsCols) =
           ListHT.unzipEithers $
@@ -87,13 +98,12 @@
              (\col ma ->
                 case ma of
                    Nothing -> Left col
-                   Just a -> Right $ Container.scale a col)
-             (Matrix.toColumns m) mas
-       lhs = Matrix.fromColumns lhsCols
-       rhs = foldl Container.add (Container.scale 0 rhs0) rhsCols
-       sol = lhs <\> Container.sub rhs0 rhs
-   in  if Vector.dim rhs0 == 0 then (map (fromMaybe 0) mas, []) else
-       (snd $
+                   Just a -> Right $ fmap (a*) col)
+             (sparseToColumns m) mas
+       lhs = sparseFromColumns (Array.bounds rhs0) lhsCols
+       rhs = foldl addSparseColumn (0 <$ rhs0) rhsCols
+       sol = QR.leastSquares lhs $ Vector.sub rhs0 rhs
+   in  (snd $
         List.mapAccumL
            (curry $ \x ->
                case x of
@@ -101,14 +111,11 @@
                   (a:as, Nothing) -> (as, a)
                   ([], Nothing) -> error "too few elements in solution vector")
            (Vector.toList sol) mas,
-        Vector.toList $
-        Container.add (lhs <> sol) rhs)
+        Vector.toList $ Vector.add (Sparse.mulVector lhs sol) rhs)
 
 
-zeroVector, _zeroVector :: Int -> Vector.Vector Double
+zeroVector :: Int -> Vector Float
 zeroVector n = Vector.fromList $ replicate n 0
--- fails for vectors of size 0
-_zeroVector n = Container.constant 0 n
 
 {-
 Approximate rotation from point correspondences.
@@ -123,58 +130,43 @@
 \0 1 y0  x0/   |c |   \y1/
                \s /
 
-We try to scale dx and dy down using 'weight'.
-Otherwise they are weighted much more than the rotation.
-However the weight will only influence the result
+Formerly, we tried to scale dx and dy down using a weight,
+such that dx and dy had magnitude similar to c and s.
+However the weight would only influence the result
 for under-constrained equation systems.
-This is usually not the case.
+This is usually not the case and
+the solver does not support rank-deficient matrices, anyway.
 -}
 layoutFromPairDisplacements ::
    [(Maybe (Float, Float), (Maybe Float, Maybe Float))] ->
    [((Int, Int), ((Float, Float), (Float, Float)))] ->
-   ([((Double,Double), Complex Double)],
-    [(Double,Double)])
+   ([((Float,Float), Complex Float)], [(Float,Float)])
 layoutFromPairDisplacements mrxys correspondences =
    let numPics = length mrxys
-       weight =
-          let xs =
-                 concatMap
-                    (\(_i, ((xai,yai),(xbi,ybi))) -> [xai, yai, xbi, ybi])
-                    correspondences
-          in  if null xs then 1 else realToFrac $ maximum xs - minimum xs
        matrix =
           sparseMatrix (2 * length correspondences) (4*numPics) $ concat $
           zipWith
-             (\k ((ia,ib), ((xai,yai),(xbi,ybi))) ->
-                let xa = realToFrac xai
-                    xb = realToFrac xbi
-                    ya = realToFrac yai
-                    yb = realToFrac ybi
-                in  elm (k+0) (4*ia+0) (-weight) :
-                    elm (k+1) (4*ia+1) (-weight) :
-                    elm (k+0) (4*ia+2) (-xa) :
-                    elm (k+0) (4*ia+3) ya :
-                    elm (k+1) (4*ia+2) (-ya) :
-                    elm (k+1) (4*ia+3) (-xa) :
-                    elm (k+0) (4*ib+0) weight :
-                    elm (k+1) (4*ib+1) weight :
-                    elm (k+0) (4*ib+2) xb :
-                    elm (k+0) (4*ib+3) (-yb) :
-                    elm (k+1) (4*ib+2) yb :
-                    elm (k+1) (4*ib+3) xb :
-                    [])
+             (\k ((ia,ib), ((xa,ya),(xb,yb))) ->
+                elm (k+0) (4*ia+0) (-1) :
+                elm (k+1) (4*ia+1) (-1) :
+                elm (k+0) (4*ia+2) (-xa) :
+                elm (k+0) (4*ia+3) ya :
+                elm (k+1) (4*ia+2) (-ya) :
+                elm (k+1) (4*ia+3) (-xa) :
+                elm (k+0) (4*ib+0) 1 :
+                elm (k+1) (4*ib+1) 1 :
+                elm (k+0) (4*ib+2) xb :
+                elm (k+0) (4*ib+3) (-yb) :
+                elm (k+1) (4*ib+2) yb :
+                elm (k+1) (4*ib+3) xb :
+                [])
              [0,2..] correspondences
        (solution, projection) =
           leastSquaresSelected matrix
-             (concatMap
-                (\(mr, (mx,my)) ->
-                   [(/weight) . realToFrac <$> mx,
-                    (/weight) . realToFrac <$> my,
-                    realToFrac . fst <$> mr,
-                    realToFrac . snd <$> mr]) $
+             (concatMap (\(mr, (mx,my)) -> [mx, my, fst <$> mr, snd <$> mr]) $
               mrxys)
              (zeroVector (2 * length correspondences))
-   in  (map (\[dx,dy,rx,ry] -> ((weight*dx,weight*dy), rx :+ ry)) $
+   in  (map (\[dx,dy,rx,ry] -> ((dx,dy), rx :+ ry)) $
         ListHT.sliceVertical 4 solution,
         map (\[x,y] -> (x,y)) $
         ListHT.sliceVertical 2 projection)
diff --git a/src/MatchImageBorders.hs b/src/MatchImageBorders.hs
--- a/src/MatchImageBorders.hs
+++ b/src/MatchImageBorders.hs
@@ -90,11 +90,15 @@
    //
    map (flip (,) locBorder) (Set.toList border)
 
+
+type
+   Queue i j =
+      MaxPQueue Float ((IOCArray (i,j) Location, CArray (i,j) Float), (i,j))
+
 prepareShaping ::
    (Ix i, Enum i, Ix j, Enum j) =>
    [(CArray (i,j) Bool, CArray (i,j) Float)] ->
-   IO ([IOCArray (i,j) Location],
-       MaxPQueue Float ((IOCArray (i,j) Location, CArray (i,j) Float), (i,j)))
+   IO ([IOCArray (i,j) Location], Queue i j)
 prepareShaping maskWeightss =
    fmap (mapSnd PQ.unions . unzip) $
    forM maskWeightss $ \(mask, weights) -> do
@@ -107,11 +111,10 @@
 shapeParts ::
    IOCArray (Int, Int) Int ->
    [IOCArray (Int, Int) Location] ->
-   MaxPQueue Float
-      ((IOCArray (Int, Int) Location, CArray (Int, Int) Float), (Int, Int)) ->
-   IO [CArray (Int, Int) Bool]
+   Queue Int Int -> IO [CArray (Int, Int) Bool]
 shapeParts count masks =
-   let loop queue =
+   let loop :: Queue Int Int -> IO [CArray (Int, Int) Bool]
+       loop queue =
          case PQ.maxView queue of
             Nothing -> mapM (fmap (amap (/=locOutside)) . freeze) masks
             Just (((locs, diffs), pos@(y,x)), remQueue) -> do