diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2014
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Makefile b/Makefile
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,5 @@
+run-test:
+	runhaskell Setup configure --user --enable-tests
+	runhaskell Setup build
+	runhaskell Setup haddock
+	./dist/build/knead-test/knead-test
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/knead.cabal b/knead.cabal
new file mode 100644
--- /dev/null
+++ b/knead.cabal
@@ -0,0 +1,119 @@
+Name:             knead
+Version:          1.0.1.1
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         https://hub.darcs.net/thielema/knead/
+Category:         Data Structures
+Synopsis:         Repa-like array processing using LLVM JIT
+Description:
+  This library processes arrays like @Repa@ and @Accelerate@,
+  but it uses the just-in-time compiler of @LLVM@
+  for generating the machine code.
+  That is, you get very efficient vectorised code
+  that can be run without a GPU.
+  You do not need to care about inlining and strictness annotations,
+  because the LLVM code is by default inlined and strict.
+  The package is intended as the basis
+  for an LLVM backend for the @Accelerate@ framework.
+  .
+  Highlights:
+  .
+  * Very flexible index handling,
+    even more flexible than the one of 'Data.Array'.
+    It is much more expressive and type-safe than that of @repa@ and @array@.
+  .
+  * Extensible element types, e.g. complex numbers.
+    (Maybe this is also possible with accelerate, e.g. with RGB type.)
+  .
+  * Every compilable program also runs.
+    In contrast to that, @accelerate@ may accept a program
+    that cannot be run by a particular backend, like @accelerate-cuda@.
+  .
+  Known deficiencies:
+  .
+  * The functions do not check array bounds.
+    (Of course, we can think about temporary bound checking
+    for debugging purposes.)
+  .
+  * The package does not try to distribute work across multiple processors.
+    It is certainly simpler, more efficient and more reliable
+    if you do that at a higher level.
+  .
+  The name of the package is inspired by the visualization of typical operations
+  like reshaping, collapsing a dimension and extruding another one.
+Tested-With:      GHC==8.4.4, GHC==8.6.5, GHC==8.10.7
+Tested-With:      GHC==9.0.2, GHC==9.2.8, GHC==9.4.6
+Cabal-Version:    >=1.10
+Build-Type:       Simple
+Extra-Source-Files:
+  Makefile
+
+Source-Repository this
+  Tag:         1.0.1.1
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/knead/
+
+Source-Repository head
+  Type:        darcs
+  Location:    https://hub.darcs.net/thielema/knead/
+
+Library
+  Build-Depends:
+    llvm-dsl >=0.1.1 && <0.2,
+    llvm-extra >=0.11 && <0.13,
+    llvm-tf >=9.0 && <17.1,
+    tfp >=1.0 && <1.1,
+    comfort-array >=0.5 && <0.6,
+    fixed-length >=0.2.1 && <0.3,
+    storable-record >=0.0.5 && <0.1,
+    storable-enum >=0.0 && <0.1,
+    bool8 >=0.0 && <0.1,
+    transformers >=0.3 && <0.7,
+    tagged >=0.7 && <0.9,
+    utility-ht >=0.0.15 && <0.1,
+    prelude-compat >=0.0 && <0.0.1,
+    base >=4 && <5
+
+  Default-Language: Haskell98
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src
+  Exposed-Modules:
+    Data.Array.Knead.Shape
+    Data.Array.Knead.Shape.Cubic
+    Data.Array.Knead.Shape.Cubic.Int
+    Data.Array.Knead.Expression
+    Data.Array.Knead.Symbolic
+    Data.Array.Knead.Symbolic.ShapeDependent
+    Data.Array.Knead.Symbolic.Physical
+    Data.Array.Knead.Symbolic.Slice
+    Data.Array.Knead.Symbolic.Fold
+    Data.Array.Knead.Symbolic.Render
+  Other-Modules:
+    Data.Array.Knead.Symbolic.RenderAlt
+    Data.Array.Knead.Symbolic.Render.Basic
+    Data.Array.Knead.Symbolic.Render.Argument
+    Data.Array.Knead.Symbolic.Private
+    Data.Array.Knead.Symbolic.PhysicalParametric
+    Data.Array.Knead.Symbolic.PhysicalPrivate
+    Data.Array.Knead.Code
+    Data.Array.Knead.Shape.Orphan
+
+Test-Suite knead-test
+  Type: exitcode-stdio-1.0
+  Build-Depends:
+    QuickCheck >=2 && <3,
+    knead,
+    comfort-array,
+    llvm-extra,
+    llvm-tf,
+    tfp,
+    utility-ht,
+    base
+  Default-Language: Haskell98
+  GHC-Options: -Wall
+  Hs-Source-Dirs: test
+  Main-Is: Main.hs
+  Other-Modules:
+    Test.Array
diff --git a/src/Data/Array/Knead/Code.hs b/src/Data/Array/Knead/Code.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Code.hs
@@ -0,0 +1,24 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Code where
+
+import qualified Data.Array.Knead.Shape as Shape
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value as MultiValue
+
+import qualified LLVM.Core as LLVM
+
+import Foreign.Ptr (Ptr)
+
+import Prelude2010
+import Prelude ()
+
+
+getElementPtr ::
+   (Shape.C sh, Shape.Index sh ~ ix, Storable.C a) =>
+   MultiValue.T sh -> LLVM.Value (Ptr a) ->
+   MultiValue.T ix ->
+   LLVM.CodeGenFunction r (LLVM.Value (Ptr a))
+getElementPtr sh ptr ix =
+   flip Storable.advancePtr ptr =<< LLVM.bitcast =<< Shape.offset sh ix
diff --git a/src/Data/Array/Knead/Expression.hs b/src/Data/Array/Knead/Expression.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Expression.hs
@@ -0,0 +1,91 @@
+module Data.Array.Knead.Expression (
+   Exp,
+   Value,
+   lift0,
+   lift1,
+   lift2,
+   lift3,
+   lift4,
+   liftM,
+   liftM2,
+   liftM3,
+   unliftM1,
+   unliftM2,
+   unliftM3,
+   liftReprM,
+   liftReprM2,
+   liftReprM3,
+   zip,
+   zip3,
+   zip4,
+   unzip,
+   unzip3,
+   unzip4,
+   fst,
+   snd,
+   mapFst,
+   mapSnd,
+   mapPair,
+   swap,
+   curry,
+   uncurry,
+   fst3,
+   snd3,
+   thd3,
+   mapFst3,
+   mapSnd3,
+   mapThd3,
+   mapTriple,
+   tuple,
+   untuple,
+   modifyMultiValue,
+   modifyMultiValue2,
+   modifyMultiValueM,
+   modifyMultiValueM2,
+   Compose(..),
+   Decompose(..),
+   modify,
+   modify2,
+   consComplex,
+   deconsComplex,
+   cons,
+   unit,
+   zero,
+   add,
+   sub,
+   mul,
+   sqr,
+   sqrt,
+   idiv,
+   irem,
+   shl,
+   shr,
+   fromInteger',
+   fromRational',
+   boolPFrom8,
+   bool8FromP,
+   intFromBool8,
+   floatFromBool8,
+   fromFastMath,
+   toFastMath,
+   minBound, maxBound,
+   cmp,
+   (==*), (/=*), (<*), (>=*), (>*), (<=*),
+   min, max,
+   true, false,
+   (&&*),
+   (||*),
+   not,
+   select,
+   ifThenElse,
+   complement,
+   (.&.*),
+   (.|.*),
+   xor,
+   toMaybe,
+   maybe,
+   ) where
+
+import LLVM.DSL.Expression
+
+import Prelude ()
diff --git a/src/Data/Array/Knead/Shape.hs b/src/Data/Array/Knead/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Shape.hs
@@ -0,0 +1,388 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE UndecidableInstances #-}
+module Data.Array.Knead.Shape (
+   C(..), Index,
+   Size,
+   value,
+   paramWith,
+   load,
+   intersect,
+   offset,
+
+   ZeroBased(ZeroBased), zeroBased, zeroBasedSize,
+
+   Range(Range), range, rangeFrom, rangeTo,
+   Shifted(Shifted), shifted, shiftedOffset, shiftedSize,
+   Cyclic(Cyclic), cyclic, cyclicSize,
+
+   Enumeration(Enumeration), EnumBounded(..),
+
+   Scalar(..),
+   Sequence(..),
+   ) where
+
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Shape.Orphan
+         (zeroBased, zeroBasedSize, cyclic, cyclicSize,
+          singletonRange, unzipRange, singletonShifted, unzipShifted)
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape
+         (Index, ZeroBased, Range(Range), Shifted(Shifted), Cyclic,
+          Enumeration(Enumeration))
+import Data.Ix (Ix)
+
+import qualified LLVM.DSL.Parameter as Param
+
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Multi.Iterator as IterMV
+import qualified LLVM.Extra.Tuple as Tuple
+import qualified LLVM.Extra.Memory as Memory
+import qualified LLVM.Extra.Iterator as Iter
+import qualified LLVM.Extra.ScalarOrVector as SoV
+import qualified LLVM.Extra.Arithmetic as A
+import LLVM.Extra.Multi.Value (atom)
+
+import qualified LLVM.Core as LLVM
+
+import qualified Data.Enum.Storable as Enum
+import Data.Tagged (Tagged)
+import Data.Tuple.HT (mapSnd)
+import Data.Word (Word8, Word16, Word32, Word64, Word)
+import Data.Int (Int8, Int16, Int32, Int64)
+
+import qualified Control.Monad.HT as Monad
+import Control.Applicative ((<$>))
+
+import Prelude2010
+import Prelude ()
+
+
+type Size = Word
+
+value :: (C sh, Expr.Value val) => sh -> val sh
+value = Expr.lift0 . MultiValue.cons
+
+paramWith ::
+   (Marshal.C b) =>
+   Param.T p b ->
+   (forall parameters.
+    (Marshal.C parameters) =>
+    (p -> parameters) ->
+    (forall val. (Expr.Value val) =>
+     MultiValue.T parameters -> val b) ->
+    a) ->
+   a
+paramWith p f =
+   Param.withMulti p (\get val -> f get (Expr.lift0 . val))
+
+load ::
+   (Marshal.C sh) =>
+   f sh -> LLVM.Value (LLVM.Ptr (Marshal.Struct sh)) ->
+   LLVM.CodeGenFunction r (MultiValue.T sh)
+load _ = Memory.load
+
+intersect :: (C sh) => Exp sh -> Exp sh -> Exp sh
+intersect = Expr.liftM2 intersectCode
+
+offset ::
+   (C sh) =>
+   MultiValue.T sh -> MultiValue.T (Index sh) ->
+   LLVM.CodeGenFunction r (LLVM.Value Size)
+offset sh ix = ($ ix) . snd =<< sizeOffset sh
+
+class (MultiValue.C sh, MultiValue.C (Index sh), Shape.Indexed sh) => C sh where
+   {-
+   It would be better to restrict zipWith to matching shapes
+   and turn shape intersection into a bound check.
+   -}
+   intersectCode ::
+      MultiValue.T sh -> MultiValue.T sh ->
+      LLVM.CodeGenFunction r (MultiValue.T sh)
+   size :: MultiValue.T sh -> LLVM.CodeGenFunction r (LLVM.Value Size)
+   {- |
+   Result is @(size, offset)@.
+   @size@ must equal the result of 'size'.
+   We use this for sharing intermediate results.
+   -}
+   sizeOffset ::
+      (Index sh ~ ix) =>
+      MultiValue.T sh ->
+      LLVM.CodeGenFunction r
+         (LLVM.Value Size,
+          MultiValue.T ix -> LLVM.CodeGenFunction r (LLVM.Value Size))
+   iterator :: (Index sh ~ ix) => MultiValue.T sh -> Iter.T r (MultiValue.T ix)
+   loop ::
+      (Index sh ~ ix, Tuple.Phi state) =>
+      (MultiValue.T ix -> state -> LLVM.CodeGenFunction r state) ->
+      MultiValue.T sh -> state -> LLVM.CodeGenFunction r state
+   loop f sh = Iter.mapState_ f (iterator sh)
+
+
+instance C () where
+   intersectCode _ _ = return $ MultiValue.cons ()
+   size _ = return A.one
+   sizeOffset _ = return (A.one, \_ -> return A.zero)
+   iterator = Iter.singleton
+   loop = id
+
+
+class C sh => Scalar sh where
+   scalar :: (Expr.Value val) => val sh
+   zeroIndex :: (Expr.Value val) => f sh -> val (Index sh)
+
+instance Scalar () where
+   scalar = Expr.lift0 $ MultiValue.Cons ()
+   zeroIndex _ = Expr.lift0 $ MultiValue.Cons ()
+
+
+class
+   (C sh,
+    MultiValue.IntegerConstant (Index sh),
+    MultiValue.Additive (Index sh)) =>
+      Sequence sh where
+   sequenceShapeFromIndex ::
+      MultiValue.T (Index sh) -> LLVM.CodeGenFunction r (MultiValue.T sh)
+
+
+class
+   (MultiValue.Additive n, MultiValue.Real n, MultiValue.IntegerConstant n) =>
+      ToSize n where
+   toSize :: MultiValue.T n -> LLVM.CodeGenFunction r (LLVM.Value Size)
+
+instance ToSize Word8  where toSize (MultiValue.Cons n) = LLVM.ext n
+instance ToSize Word16 where toSize (MultiValue.Cons n) = LLVM.ext n
+instance ToSize Word32 where toSize (MultiValue.Cons n) = LLVM.adapt n
+instance ToSize Word64 where toSize (MultiValue.Cons n) = LLVM.adapt n
+instance ToSize Word   where toSize (MultiValue.Cons n) = LLVM.adapt n
+instance ToSize Int8  where toSize (MultiValue.Cons n) = LLVM.zext n
+instance ToSize Int16 where toSize (MultiValue.Cons n) = LLVM.zext n
+instance ToSize Int32 where toSize (MultiValue.Cons n) = LLVM.zadapt n
+instance ToSize Int64 where toSize (MultiValue.Cons n) = LLVM.zadapt n
+instance ToSize Int   where toSize (MultiValue.Cons n) = LLVM.zadapt n
+
+
+{- |
+Array dimensions and indexes cannot be negative,
+but computations in indices may temporarily yield negative values
+or we want to add negative values to indices.
+
+So maybe, we would better have type Index (ZeroBased Word64) = Int64.
+This is not possible.
+Maybe we need an additional ZeroBased type for unsigned array sizes.
+-}
+instance
+      (Integral n, ToSize n, MultiValue.Comparison n) => C (ZeroBased n) where
+   intersectCode sha shb =
+      zeroBased <$> MultiValue.min (zeroBasedSize sha) (zeroBasedSize shb)
+   size = toSize . zeroBasedSize
+   sizeOffset sh = Monad.lift2 (,) (toSize $ zeroBasedSize sh) (return toSize)
+   iterator sh =
+      IterMV.take (zeroBasedSize sh) $
+      Iter.iterate MultiValue.inc MultiValue.zero
+
+instance
+   (Integral n, ToSize n, MultiValue.Comparison n) =>
+      Sequence (ZeroBased n) where
+   sequenceShapeFromIndex = return . zeroBased
+
+
+rangeSize ::
+   (ToSize n) =>
+   Range (MultiValue.T n) -> LLVM.CodeGenFunction r (LLVM.Value Size)
+rangeSize (Range from to) =
+   toSize =<< MultiValue.inc =<< MultiValue.sub to from
+
+
+rangeFrom :: (Expr.Value val) => val (Range n) -> val n
+rangeFrom = Expr.lift1 $ Shape.rangeFrom . unzipRange
+
+rangeTo :: (Expr.Value val) => val (Range n) -> val n
+rangeTo = Expr.lift1 $ Shape.rangeTo . unzipRange
+
+range :: (Expr.Value val) => val n -> val n -> val (Range n)
+range =
+   Expr.lift2 $
+      \(MultiValue.Cons from) (MultiValue.Cons to) ->
+         MultiValue.Cons (Range from to)
+
+instance (Ix n, ToSize n, MultiValue.Comparison n) => C (Range n) where
+   intersectCode =
+      MultiValue.modifyF2 (singletonRange atom) (singletonRange atom) $
+            \(Range fromN toN) (Range fromM toM) ->
+         Monad.lift2 Range (MultiValue.max fromN fromM) (MultiValue.min toN toM)
+   size = rangeSize . unzipRange
+   sizeOffset rngValue =
+      case unzipRange rngValue of
+         rng@(Range from _to) ->
+            Monad.lift2 (,) (rangeSize rng)
+               (return $ \i -> toSize =<< MultiValue.sub i from)
+   iterator rngValue =
+      case MultiValue.decompose (singletonRange atom) rngValue of
+         Range from to ->
+            IterMV.takeWhile (MultiValue.cmp LLVM.CmpGE to) $
+            Iter.iterate MultiValue.inc from
+
+
+
+shiftedOffset :: (Expr.Value val) => val (Shifted n) -> val n
+shiftedOffset = Expr.lift1 $ Shape.shiftedOffset . unzipShifted
+
+shiftedSize :: (Expr.Value val) => val (Shifted n) -> val n
+shiftedSize = Expr.lift1 $ Shape.shiftedSize . unzipShifted
+
+shifted :: (Expr.Value val) => val n -> val n -> val (Shifted n)
+shifted =
+   Expr.lift2 $
+      \(MultiValue.Cons from) (MultiValue.Cons to) ->
+         MultiValue.Cons (Shifted from to)
+
+
+instance (Integral n, ToSize n, MultiValue.Comparison n) => C (Shifted n) where
+   intersectCode =
+      MultiValue.modifyF2 (singletonShifted atom) (singletonShifted atom) $
+            \(Shifted startN lenN) (Shifted startM lenM) -> do
+         start <- MultiValue.max startN startM
+         endN <- MultiValue.add startN lenN
+         endM <- MultiValue.add startM lenM
+         end <- MultiValue.min endN endM
+         Shifted start <$> MultiValue.sub end start
+   size = toSize . shiftedSize
+   sizeOffset shapeValue =
+      case unzipShifted shapeValue of
+         Shifted start len ->
+            Monad.lift2 (,) (toSize len)
+               (return $ \i -> toSize =<< MultiValue.sub i start)
+   iterator rngValue =
+      case MultiValue.decompose (singletonShifted atom) rngValue of
+         Shifted from len ->
+            IterMV.take len $ Iter.iterate MultiValue.inc from
+
+
+instance
+      (Integral n, ToSize n, MultiValue.Comparison n) => C (Cyclic n) where
+   intersectCode sha shb =
+      cyclic <$> MultiValue.min (cyclicSize sha) (cyclicSize shb)
+   size = toSize . cyclicSize
+   sizeOffset sh = Monad.lift2 (,) (toSize $ cyclicSize sh) (return toSize)
+   iterator sh =
+      IterMV.take (cyclicSize sh) $
+      Iter.iterate MultiValue.inc MultiValue.zero
+
+
+class (IterMV.Enum enum, MultiValue.Bounded enum) => EnumBounded enum where
+   enumOffset :: MultiValue.T enum -> LLVM.CodeGenFunction r (LLVM.Value Size)
+
+instance
+   (ToSize w, MultiValue.Additive w,
+    LLVM.IsInteger w, SoV.IntegerConstant w, Num w,
+    MultiValue.Repr w ~ LLVM.Value w,
+    LLVM.CmpRet w, LLVM.IsPrimitive w,
+    Enum e, Bounded e) =>
+      EnumBounded (Enum.T w e) where
+   enumOffset ix =
+      toSize =<<
+      MultiValue.sub
+         (MultiValue.fromEnum ix)
+         (MultiValue.fromEnum $ MultiValue.minBound `asTypeOf` ix)
+
+instance
+      (Enum enum, Bounded enum, EnumBounded enum) => C (Enumeration enum) where
+   intersectCode _sha shb = return shb
+   size = return . A.fromInteger' . toInteger . Shape.size . plainEnumeration
+   sizeOffset sh = do
+      sz <- size sh
+      return (sz, enumOffset)
+   iterator _ = IterMV.enumFromTo MultiValue.minBound MultiValue.maxBound
+
+plainEnumeration :: val (Enumeration enum) -> Enumeration enum
+plainEnumeration _ = Enumeration
+
+
+instance (C sh) => C (Tagged tag sh) where
+   intersectCode = MultiValue.liftTaggedM2 intersectCode
+   size = size . MultiValue.untag
+   sizeOffset =
+      fmap (mapSnd (. MultiValue.untag)) . sizeOffset . MultiValue.untag
+   iterator = fmap MultiValue.tag . iterator . MultiValue.untag
+
+
+instance (C n, C m) => C (n,m) where
+   intersectCode a b =
+      case (MultiValue.unzip a, MultiValue.unzip b) of
+         ((an,am), (bn,bm)) ->
+            Monad.lift2 MultiValue.zip
+               (intersectCode an bn)
+               (intersectCode am bm)
+   size nm =
+      case MultiValue.unzip nm of
+         (n,m) -> Monad.liftJoin2 A.mul (size n) (size m)
+   sizeOffset nm =
+      case MultiValue.unzip nm of
+         (n,m) -> do
+            (ns, iOffset) <- sizeOffset n
+            (ms, jOffset) <- sizeOffset m
+            sz <- A.mul ns ms
+            return
+               (sz,
+                \ij ->
+                  case MultiValue.unzip ij of
+                     (i,j) -> do
+                        il <- iOffset i
+                        jl <- jOffset j
+                        A.add jl =<< A.mul ms il)
+   iterator nm =
+      case MultiValue.unzip nm of
+         (n,m) ->
+            uncurry MultiValue.zip <$>
+            Iter.cartesian (iterator n) (iterator m)
+   loop code nm =
+      case MultiValue.unzip nm of
+         (n,m) -> loop (\i -> loop (\j -> code (MultiValue.zip i j)) m) n
+
+instance (C n, C m, C l) => C (n,m,l) where
+   intersectCode a b =
+      case (MultiValue.unzip3 a, MultiValue.unzip3 b) of
+         ((ai,aj,ak), (bi,bj,bk)) ->
+            Monad.lift3 MultiValue.zip3
+               (intersectCode ai bi)
+               (intersectCode aj bj)
+               (intersectCode ak bk)
+   size nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            Monad.liftJoin2 A.mul (size n) $
+            Monad.liftJoin2 A.mul (size m) (size l)
+   sizeOffset nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) -> do
+            (ns, iOffset) <- sizeOffset n
+            (ms, jOffset) <- sizeOffset m
+            (ls, kOffset) <- sizeOffset l
+            sz <- A.mul ns =<< A.mul ms ls
+            return
+               (sz,
+                \ijk ->
+                  case MultiValue.unzip3 ijk of
+                     (i,j,k) -> do
+                        il <- iOffset i
+                        jl <- jOffset j
+                        kl <- kOffset k
+                        A.add kl =<< A.mul ls =<< A.add jl =<< A.mul ms il)
+   iterator nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            fmap (\(a,(b,c)) -> MultiValue.zip3 a b c) $
+            Iter.cartesian (iterator n) $
+            Iter.cartesian (iterator m) (iterator l)
+   loop code nml =
+      case MultiValue.unzip3 nml of
+         (n,m,l) ->
+            loop (\i -> loop (\j -> loop (\k ->
+               code (MultiValue.zip3 i j k))
+            l) m) n
diff --git a/src/Data/Array/Knead/Shape/Cubic.hs b/src/Data/Array/Knead/Shape/Cubic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Shape/Cubic.hs
@@ -0,0 +1,328 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE UndecidableInstances #-}
+module Data.Array.Knead.Shape.Cubic (
+   constant,
+   paramWith,
+   tunnel,
+
+   T(..),
+   Z(Z), z,
+   (:.)((:.)),
+   Shape,
+   Index,
+   cons, (#:.),
+   head,
+   tail,
+   switchR,
+   ) where
+
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Shape.Cubic.Int as Index
+
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified Data.Array.Comfort.Shape as ComfortShape
+import Data.Array.Comfort.Shape (ZeroBased(ZeroBased))
+
+import qualified LLVM.DSL.Parameter as Param
+
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Multi.Iterator as IterMV
+import qualified LLVM.Extra.Iterator as Iter
+import qualified LLVM.Extra.Arithmetic as A
+import qualified LLVM.Extra.Tuple as Tuple
+import qualified LLVM.Extra.Control as C
+import LLVM.Extra.Multi.Value (Atom)
+
+import qualified LLVM.Core as LLVM
+
+import qualified Foreign.Storable as St
+import Foreign.Storable.FixedArray (sizeOfArray, )
+import Foreign.Ptr (castPtr, )
+
+import qualified Type.Data.Num.Decimal as Dec
+import qualified Type.Data.Num.Unary as Unary
+import Type.Base.Proxy (Proxy(Proxy))
+
+import qualified Data.Traversable as Trav
+import qualified Data.Foldable as Fold
+import qualified Data.FixedLength as FixedLength
+import Data.FixedLength ((!:))
+
+import Control.Monad (liftM2, )
+import Control.Applicative (pure, (<$>), )
+
+import Prelude hiding (min, head, tail, )
+
+
+newtype T tag rank = Cons {decons :: FixedLength.T rank Index.Int}
+
+data ShapeTag
+data IndexTag
+
+type Shape = T ShapeTag
+type Index = T IndexTag
+
+
+paramWith ::
+   (Unary.Natural rank,
+    Dec.Natural (Dec.FromUnary rank),
+    Dec.Natural (Dec.FromUnary rank Dec.:*: LLVM.SizeOf Shape.Size)) =>
+   Param.T p (T tag rank) ->
+   (forall parameters.
+    (Marshal.C parameters) =>
+    (p -> parameters) ->
+    (forall val. (Expr.Value val) =>
+     MultiValue.T parameters -> val (T tag rank)) ->
+    a) ->
+   a
+paramWith p f =
+   case tunnel p of
+      Param.Tunnel get val -> f get (Expr.lift0 . val)
+
+tunnel ::
+   (Unary.Natural rank,
+    Dec.Natural (Dec.FromUnary rank),
+    Dec.Natural (Dec.FromUnary rank Dec.:*: LLVM.SizeOf Shape.Size)) =>
+   Param.T p (T tag rank) -> Param.Tunnel p (T tag rank)
+tunnel p = Param.tunnel MultiValue.cons p
+
+
+data Z = Z
+   deriving (Eq, Ord, Read, Show)
+
+
+infixl 3 :., #:.
+
+data tail :. head = !tail :. !head
+   deriving (Eq, Ord, Read, Show)
+
+
+(#:.) ::
+   (Expr.Value val) =>
+   val (T tag rank) -> val Index.Int -> val (T tag (Unary.Succ rank))
+(#:.) = cons
+
+cons ::
+   (Expr.Value val) =>
+   val (T tag rank) -> val Index.Int -> val (T tag (Unary.Succ rank))
+cons =
+   Expr.lift2 $
+      \(MultiValue.Cons t) (MultiValue.Cons h) -> MultiValue.Cons (h!:t)
+
+z :: (Expr.Value val) => val (T tag Unary.Zero)
+z = Expr.lift0 $ MultiValue.Cons FixedLength.end
+
+head ::
+   (Expr.Value val, Unary.Natural rank) =>
+   val (T tag (Unary.Succ rank)) -> val Index.Int
+head =
+   Expr.lift1 $ \(MultiValue.Cons sh) -> MultiValue.Cons $ FixedLength.head sh
+
+tail ::
+   (Expr.Value val, Unary.Natural rank) =>
+   val (T tag (Unary.Succ rank)) -> val (T tag rank)
+tail =
+   Expr.lift1 $ \(MultiValue.Cons sh) -> MultiValue.Cons $ FixedLength.tail sh
+
+switchR ::
+   (Unary.Natural rank) =>
+   Expr.Value val =>
+   (val (T tag rank) -> val Index.Int -> a) ->
+   val (T tag (Unary.Succ rank)) -> a
+switchR f ix = f (tail ix) (head ix)
+
+
+rank :: T tag rank -> Proxy rank
+rank (Cons _) = Proxy
+
+
+instance (tag ~ ShapeTag, rank ~ Unary.Zero) => Shape.Scalar (T tag rank) where
+   scalar = Expr.lift0 $ MultiValue.Cons FixedLength.end
+   zeroIndex _ = Expr.lift0 $ MultiValue.Cons FixedLength.end
+
+
+type family AtomRank sh
+type instance AtomRank (Atom (T tag rank)) = rank
+type instance AtomRank (sh:.s) = Unary.Succ (AtomRank s)
+
+type family AtomTag sh
+type instance AtomTag (Atom (T tag rank)) = tag
+type instance AtomTag (sh:.s) = AtomTag sh
+
+type instance MultiValue.PatternTuple (sh:.s) =
+   T (AtomTag sh) (Unary.Succ (AtomRank sh))
+
+type instance MultiValue.Decomposed f (sh:.s) =
+   MultiValue.Decomposed f sh :. f Index.Int
+
+instance
+   (Expr.Decompose sh, Expr.Decompose s,
+    MultiValue.Decomposed Exp s ~ Exp Index.Int,
+    MultiValue.PatternTuple s ~ Index.Int,
+    MultiValue.PatternTuple sh ~ T (AtomTag sh) (AtomRank sh),
+    Unary.Natural (AtomRank sh)) =>
+      Expr.Decompose (sh :. s) where
+   decompose (psh:.ps) x =
+      Expr.decompose psh (tail x) :. Expr.decompose ps (head x)
+
+
+type family Rank sh
+type instance Rank (T tag rank) = rank
+
+type family Tag sh
+type instance Tag (T tag rank) = tag
+
+instance
+   (Expr.Compose sh,
+    Expr.Composed sh ~ T (Tag (Expr.Composed sh)) (Rank (Expr.Composed sh)),
+    Expr.Compose s,
+    Expr.Composed s ~ Index.Int) =>
+      Expr.Compose (sh :. s) where
+   type Composed (sh :. s) =
+            T (Tag (Expr.Composed sh)) (Unary.Succ (Rank (Expr.Composed sh)))
+   compose (sh :. s) = cons (Expr.compose sh) (Expr.compose s)
+
+
+instance (Unary.Natural rank) => St.Storable (T tag rank) where
+   sizeOf sh = sizeOfArray (Unary.integralFromProxy $ rank sh) (0::Shape.Size)
+   alignment (Cons _sh) = St.alignment (0::Shape.Size)
+   poke ptr = St.poke (castPtr ptr) . fmap (\(Index.Int i) -> i) . decons
+   peek = fmap (Cons . fmap Index.Int) . St.peek . castPtr
+
+instance
+   (Unary.Natural rank,
+    Dec.Natural (Dec.FromUnary rank),
+    Dec.Natural (Dec.FromUnary rank Dec.:*: LLVM.SizeOf Shape.Size)) =>
+      Marshal.C (T tag rank) where
+   pack = LLVM.Array . map Marshal.pack . Fold.toList . decons
+   unpack (LLVM.Array sh) = Cons $ toFixedList $ map Marshal.unpack sh
+
+toFixedList :: (Unary.Natural n) => [a] -> FixedLength.T n a
+toFixedList xs = snd $ Trav.mapAccumL (\(y:ys) () -> (ys,y)) xs (pure ())
+
+
+instance (Unary.Natural rank) => MultiValue.C (T tag rank) where
+   type Repr (T tag rank) = FixedLength.T rank (MultiValue.Repr Index.Int)
+   cons = MultiValue.Cons . fmap (\(Index.Int i) -> LLVM.valueOf i) . decons
+   undef = constant $ MultiValue.undef
+   zero = constant $ MultiValue.zero
+   addPhi bb (MultiValue.Cons a) (MultiValue.Cons b) =
+      Tuple.addPhiFoldable bb a b
+   phi bb (MultiValue.Cons a) =
+      fmap MultiValue.Cons . Tuple.phiTraversable bb $ a
+
+constant ::
+   (Unary.Natural rank) => MultiValue.T Index.Int -> MultiValue.T (T tag rank)
+constant (MultiValue.Cons x) = MultiValue.Cons $ pure x
+
+instance
+   (tag ~ ShapeTag, Unary.Natural rank) =>
+      ComfortShape.C (T tag rank) where
+   size = Fold.product . fmap (ComfortShape.size . shapeFromInt) . decons
+
+instance
+   (tag ~ ShapeTag, Unary.Natural rank) =>
+      ComfortShape.Indexed (T tag rank) where
+   type Index (T tag rank) = Index rank
+   indices (Cons ix) =
+      map (Cons . fmap Index.Int) $
+      Trav.mapM (ComfortShape.indices . shapeFromInt) ix
+   inBounds (Cons sh) (Cons ix) =
+      Fold.and $
+      FixedLength.zipWith ComfortShape.inBounds
+         (shapeFromInt <$> sh) (indexFromInt <$> ix)
+   unifiedOffset (Cons sh) (Cons ix) =
+      Fold.foldlM
+         (\off (s,i) -> do
+            ioff <- ComfortShape.unifiedOffset s i
+            return $! off * ComfortShape.size s + ioff)
+         0 $
+      FixedLength.zipWith (,) (shapeFromInt <$> sh) (indexFromInt <$> ix)
+
+shapeFromInt :: Index.Int -> ZeroBased Shape.Size
+shapeFromInt (Index.Int i) = ZeroBased i
+
+indexFromInt :: Index.Int -> Shape.Size
+indexFromInt (Index.Int i) = i
+
+
+instance (tag ~ ShapeTag, Unary.Natural rank) => Shape.C (T tag rank) where
+   size (MultiValue.Cons sh) = Fold.foldlM A.mul A.one sh
+   intersectCode (MultiValue.Cons sh0) (MultiValue.Cons sh1) =
+      fmap MultiValue.Cons $ Trav.sequence $ FixedLength.zipWith A.min sh0 sh1
+   sizeOffset sh =
+      -- would a joint implementation be more efficient?
+      liftM2 (,) (Shape.size sh) (return $ offsetCode sh)
+   iterator = iterator
+   loop = loop
+
+
+offsetCode ::
+   (Unary.Natural rank) =>
+   MultiValue.T (Shape rank) -> MultiValue.T (Index rank) ->
+   LLVM.CodeGenFunction r (LLVM.Value Shape.Size)
+offsetCode (MultiValue.Cons sh) (MultiValue.Cons ix) =
+   Fold.foldlM (\off (s,i) -> A.mul off s >>= A.add i) A.zero $
+   FixedLength.zipWith (,) sh ix
+
+
+newtype Iterator r rank =
+   Iterator {
+      runIterator ::
+         MultiValue.T (Shape rank) -> Iter.T r (MultiValue.T (Index rank))
+   }
+
+iterator ::
+   (Unary.Natural rank) =>
+   MultiValue.T (Shape rank) -> Iter.T r (MultiValue.T (Index rank))
+iterator =
+   runIterator $
+   Unary.switchNat
+      (Iterator $ \ _z -> Iter.singleton z)
+      (Iterator $ switchR $ \sh n ->
+       fmap (\(ix,i) -> ix#:.i) $
+       Iter.cartesian
+         (iterator sh)
+         (IterMV.takeWhile (MultiValue.cmp LLVM.CmpGT n) $
+          Iter.iterate MultiValue.inc MultiValue.zero))
+
+
+newtype Loop r state rank =
+   Loop {
+      runLoop ::
+         (MultiValue.T (Index rank) ->
+          state ->
+          LLVM.CodeGenFunction r state) ->
+         MultiValue.T (Shape rank) ->
+         state ->
+         LLVM.CodeGenFunction r state
+   }
+
+loop ::
+   (Unary.Natural rank, Tuple.Phi state) =>
+   (MultiValue.T (Index rank) ->
+    state ->
+    LLVM.CodeGenFunction r state) ->
+   MultiValue.T (Shape rank) ->
+   state ->
+   LLVM.CodeGenFunction r state
+loop =
+   runLoop $
+   Unary.switchNat
+      (Loop $ \code _z -> code z)
+      (Loop $ \code -> switchR $ \sh (MultiValue.Cons n) ->
+         loop
+            (\ix ptrStart ->
+               fmap fst $
+               C.fixedLengthLoop n (ptrStart, A.zero) $ \(ptr, k) ->
+                  liftM2 (,)
+                     (code (ix #:. MultiValue.Cons k) ptr)
+                     (A.inc k))
+            sh)
diff --git a/src/Data/Array/Knead/Shape/Cubic/Int.hs b/src/Data/Array/Knead/Shape/Cubic/Int.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Shape/Cubic/Int.hs
@@ -0,0 +1,67 @@
+{-# LANGUAGE TypeFamilies #-}
+module Data.Array.Knead.Shape.Cubic.Int (
+   Single(..),
+   Int(Int), cons, decons,
+   ) where
+
+import qualified Data.Array.Knead.Expression as Expr
+
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Arithmetic as A
+
+import qualified LLVM.Core as LLVM
+
+import Data.Word (Word)
+
+import Prelude hiding (Int, head, tail, )
+
+
+newtype Int = Int Word
+
+cons :: (Expr.Value val) => val Word -> val Int
+cons = Expr.lift1 $ \(MultiValue.Cons x) -> MultiValue.Cons x
+
+decons :: (Expr.Value val) => val Int -> val Word
+decons = Expr.lift1 $ \(MultiValue.Cons x) -> MultiValue.Cons x
+
+
+class Single ix where
+   switchSingle :: f Int -> f ix
+
+instance Single Int where
+   switchSingle x = x
+
+
+instance MultiValue.C Int where
+   type Repr Int = LLVM.Value Word
+   cons (Int x) = MultiValue.consPrimitive x
+   undef = MultiValue.undefPrimitive
+   zero = MultiValue.zeroPrimitive
+   phi = MultiValue.phiPrimitive
+   addPhi = MultiValue.addPhiPrimitive
+
+instance MultiValue.Additive Int where
+   add = MultiValue.liftM2 A.add
+   sub = MultiValue.liftM2 A.sub
+   neg = MultiValue.liftM A.neg
+
+instance MultiValue.PseudoRing Int where
+   mul = MultiValue.liftM2 A.mul
+
+instance MultiValue.Real Int where
+   min = MultiValue.liftM2 A.min
+   max = MultiValue.liftM2 A.max
+   abs = MultiValue.liftM A.abs
+   signum = MultiValue.liftM A.signum
+
+instance MultiValue.IntegerConstant Int where
+   fromInteger' = cons . A.fromInteger'
+
+instance MultiValue.Comparison Int where
+   cmp mode = MultiValue.liftM2 $ A.cmp mode
+
+
+instance Marshal.C Int where
+   pack (Int i) = i
+   unpack = Int
diff --git a/src/Data/Array/Knead/Shape/Orphan.hs b/src/Data/Array/Knead/Shape/Orphan.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Shape/Orphan.hs
@@ -0,0 +1,281 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+module Data.Array.Knead.Shape.Orphan where
+
+import qualified Data.Array.Knead.Expression as Expr
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape
+         (ZeroBased(ZeroBased), Range(Range), Shifted(Shifted),
+          Cyclic(Cyclic),
+          Enumeration(Enumeration))
+
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Memory as Memory
+import qualified LLVM.Extra.Tuple as Tuple
+
+import qualified Control.Monad.HT as Monad
+import Control.Applicative ((<$>))
+
+import Prelude2010
+import Prelude ()
+
+
+
+unzipZeroBased :: MultiValue.T (ZeroBased n) -> ZeroBased (MultiValue.T n)
+unzipZeroBased (MultiValue.Cons (ZeroBased n)) = ZeroBased (MultiValue.Cons n)
+
+zeroBasedSize :: (Expr.Value val) => val (ZeroBased n) -> val n
+zeroBasedSize = Expr.lift1 $ Shape.zeroBasedSize . unzipZeroBased
+
+zeroBased :: (Expr.Value val) => val n -> val (ZeroBased n)
+zeroBased = Expr.lift1 $ \(MultiValue.Cons n) -> MultiValue.Cons (ZeroBased n)
+
+instance (Tuple.Undefined n) => Tuple.Undefined (ZeroBased n) where
+   undef = ZeroBased Tuple.undef
+
+instance (Tuple.Phi n) => Tuple.Phi (ZeroBased n) where
+   phi bb = fmap ZeroBased . Tuple.phi bb . Shape.zeroBasedSize
+   addPhi bb (Shape.ZeroBased a) (Shape.ZeroBased b) = Tuple.addPhi bb a b
+
+instance (Tuple.Value n) => Tuple.Value (ZeroBased n) where
+   type ValueOf (ZeroBased n) = ZeroBased (Tuple.ValueOf n)
+   valueOf (ZeroBased n) = ZeroBased $ Tuple.valueOf n
+
+instance (MultiValue.C n) => MultiValue.C (ZeroBased n) where
+   type Repr (ZeroBased n) = ZeroBased (MultiValue.Repr n)
+   cons (ZeroBased n) = zeroBased (MultiValue.cons n)
+   undef = zeroBased MultiValue.undef
+   zero = zeroBased MultiValue.zero
+   phi bb = Monad.lift zeroBased . MultiValue.phi bb . zeroBasedSize
+   addPhi bb a b = MultiValue.addPhi bb (zeroBasedSize a) (zeroBasedSize b)
+
+type instance
+   MultiValue.Decomposed f (ZeroBased pn) =
+      ZeroBased (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (ZeroBased pn) =
+      ZeroBased (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (ZeroBased n) where
+   type Composed (ZeroBased n) = ZeroBased (MultiValue.Composed n)
+   compose (ZeroBased n) = zeroBased (MultiValue.compose n)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (ZeroBased pn) where
+   decompose (ZeroBased p) sh =
+      MultiValue.decompose p <$> unzipZeroBased sh
+
+instance (Expr.Compose n) => Expr.Compose (ZeroBased n) where
+   type Composed (ZeroBased n) = ZeroBased (Expr.Composed n)
+   compose (ZeroBased n) = Expr.lift1 zeroBased (Expr.compose n)
+
+instance (Expr.Decompose pn) => Expr.Decompose (ZeroBased pn) where
+   decompose (ZeroBased p) = ZeroBased . Expr.decompose p . zeroBasedSize
+
+instance (Memory.C n) => Memory.C (ZeroBased n) where
+   type Struct (ZeroBased n) = Memory.Struct n
+   compose = Memory.compose . Shape.zeroBasedSize
+   decompose = fmap ZeroBased . Memory.decompose
+
+instance (Marshal.C n) => Marshal.C (ZeroBased n) where
+   pack = Marshal.pack . Shape.zeroBasedSize
+   unpack = Shape.ZeroBased . Marshal.unpack
+
+
+
+singletonRange :: n -> Range n
+singletonRange n = Range n n
+
+unzipRange :: MultiValue.T (Range n) -> Range (MultiValue.T n)
+unzipRange (MultiValue.Cons (Range from to)) =
+   Range (MultiValue.Cons from) (MultiValue.Cons to)
+
+zipRange :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Range n)
+zipRange (MultiValue.Cons from) (MultiValue.Cons to) =
+   MultiValue.Cons (Range from to)
+
+instance (Tuple.Undefined n) => Tuple.Undefined (Range n) where
+   undef = Range Tuple.undef Tuple.undef
+
+instance (Tuple.Value n) => Tuple.Value (Range n) where
+   type ValueOf (Range n) = Range (Tuple.ValueOf n)
+   valueOf (Range from to) = Range (Tuple.valueOf from) (Tuple.valueOf to)
+
+instance (MultiValue.C n) => MultiValue.C (Range n) where
+   type Repr (Range n) = Range (MultiValue.Repr n)
+   cons (Range from to) = zipRange (MultiValue.cons from) (MultiValue.cons to)
+   undef = MultiValue.compose $ singletonRange MultiValue.undef
+   zero = MultiValue.compose $ singletonRange MultiValue.zero
+   phi bb a =
+      case unzipRange a of
+         Range a0 a1 ->
+            Monad.lift2 zipRange (MultiValue.phi bb a0) (MultiValue.phi bb a1)
+   addPhi bb a b =
+      case (unzipRange a, unzipRange b) of
+         (Range a0 a1, Range b0 b1) ->
+            MultiValue.addPhi bb a0 b0 >>
+            MultiValue.addPhi bb a1 b1
+
+type instance
+   MultiValue.Decomposed f (Range pn) = Range (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (Range pn) = Range (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Range n) where
+   type Composed (Range n) = Range (MultiValue.Composed n)
+   compose (Range from to) =
+      zipRange (MultiValue.compose from) (MultiValue.compose to)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Range pn) where
+   decompose (Range pfrom pto) rng =
+      case unzipRange rng of
+         Range from to ->
+            Range
+               (MultiValue.decompose pfrom from)
+               (MultiValue.decompose pto to)
+
+
+
+singletonShifted :: n -> Shifted n
+singletonShifted n = Shifted n n
+
+unzipShifted :: MultiValue.T (Shifted n) -> Shifted (MultiValue.T n)
+unzipShifted (MultiValue.Cons (Shifted from to)) =
+   Shifted (MultiValue.Cons from) (MultiValue.Cons to)
+
+zipShifted :: MultiValue.T n -> MultiValue.T n -> MultiValue.T (Shifted n)
+zipShifted (MultiValue.Cons from) (MultiValue.Cons to) =
+   MultiValue.Cons (Shifted from to)
+
+instance (Tuple.Undefined n) => Tuple.Undefined (Shifted n) where
+   undef = Shifted Tuple.undef Tuple.undef
+
+instance (Tuple.Value n) => Tuple.Value (Shifted n) where
+   type ValueOf (Shifted n) = Shifted (Tuple.ValueOf n)
+   valueOf (Shifted start len) =
+      Shifted (Tuple.valueOf start) (Tuple.valueOf len)
+
+instance (MultiValue.C n) => MultiValue.C (Shifted n) where
+   type Repr (Shifted n) = Shifted (MultiValue.Repr n)
+   cons (Shifted start len) =
+      zipShifted (MultiValue.cons start) (MultiValue.cons len)
+   undef = MultiValue.compose $ singletonShifted MultiValue.undef
+   zero = MultiValue.compose $ singletonShifted MultiValue.zero
+   phi bb a =
+      case unzipShifted a of
+         Shifted a0 a1 ->
+            Monad.lift2 zipShifted
+               (MultiValue.phi bb a0) (MultiValue.phi bb a1)
+   addPhi bb a b =
+      case (unzipShifted a, unzipShifted b) of
+         (Shifted a0 a1, Shifted b0 b1) ->
+            MultiValue.addPhi bb a0 b0 >>
+            MultiValue.addPhi bb a1 b1
+
+type instance
+   MultiValue.Decomposed f (Shifted pn) =
+      Shifted (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (Shifted pn) =
+      Shifted (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Shifted n) where
+   type Composed (Shifted n) = Shifted (MultiValue.Composed n)
+   compose (Shifted start len) =
+      zipShifted (MultiValue.compose start) (MultiValue.compose len)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Shifted pn) where
+   decompose (Shifted pstart plen) rng =
+      case unzipShifted rng of
+         Shifted start len ->
+            Shifted
+               (MultiValue.decompose pstart start)
+               (MultiValue.decompose plen len)
+
+
+
+unzipCyclic :: MultiValue.T (Cyclic n) -> Cyclic (MultiValue.T n)
+unzipCyclic (MultiValue.Cons (Cyclic n)) = Cyclic (MultiValue.Cons n)
+
+cyclicSize :: (Expr.Value val) => val (Cyclic n) -> val n
+cyclicSize = Expr.lift1 $ Shape.cyclicSize . unzipCyclic
+
+cyclic :: (Expr.Value val) => val n -> val (Cyclic n)
+cyclic = Expr.lift1 $ \(MultiValue.Cons n) -> MultiValue.Cons (Cyclic n)
+
+instance (Tuple.Undefined n) => Tuple.Undefined (Cyclic n) where
+   undef = Cyclic Tuple.undef
+
+instance (Tuple.Phi n) => Tuple.Phi (Cyclic n) where
+   phi bb = fmap Cyclic . Tuple.phi bb . Shape.cyclicSize
+   addPhi bb (Shape.Cyclic a) (Shape.Cyclic b) = Tuple.addPhi bb a b
+
+instance (Tuple.Value n) => Tuple.Value (Cyclic n) where
+   type ValueOf (Cyclic n) = Cyclic (Tuple.ValueOf n)
+   valueOf (Cyclic n) = Cyclic $ Tuple.valueOf n
+
+instance (MultiValue.C n) => MultiValue.C (Cyclic n) where
+   type Repr (Cyclic n) = Cyclic (MultiValue.Repr n)
+   cons (Cyclic n) = cyclic (MultiValue.cons n)
+   undef = cyclic MultiValue.undef
+   zero = cyclic MultiValue.zero
+   phi bb = Monad.lift cyclic . MultiValue.phi bb . cyclicSize
+   addPhi bb a b = MultiValue.addPhi bb (cyclicSize a) (cyclicSize b)
+
+type instance
+   MultiValue.Decomposed f (Cyclic pn) = Cyclic (MultiValue.Decomposed f pn)
+type instance
+   MultiValue.PatternTuple (Cyclic pn) = Cyclic (MultiValue.PatternTuple pn)
+
+instance (MultiValue.Compose n) => MultiValue.Compose (Cyclic n) where
+   type Composed (Cyclic n) = Cyclic (MultiValue.Composed n)
+   compose (Cyclic n) = cyclic (MultiValue.compose n)
+
+instance (MultiValue.Decompose pn) => MultiValue.Decompose (Cyclic pn) where
+   decompose (Cyclic p) sh = MultiValue.decompose p <$> unzipCyclic sh
+
+instance (Expr.Compose n) => Expr.Compose (Cyclic n) where
+   type Composed (Cyclic n) = Cyclic (Expr.Composed n)
+   compose (Cyclic n) = Expr.lift1 cyclic (Expr.compose n)
+
+instance (Expr.Decompose pn) => Expr.Decompose (Cyclic pn) where
+   decompose (Cyclic p) = Cyclic . Expr.decompose p . cyclicSize
+
+instance (Memory.C n) => Memory.C (Cyclic n) where
+   type Struct (Cyclic n) = Memory.Struct n
+   compose = Memory.compose . Shape.cyclicSize
+   decompose = fmap Cyclic . Memory.decompose
+
+instance (Marshal.C n) => Marshal.C (Cyclic n) where
+   pack = Marshal.pack . Shape.cyclicSize
+   unpack = Shape.Cyclic . Marshal.unpack
+
+
+
+instance (Enum enum, Bounded enum) => MultiValue.C (Enumeration enum) where
+   type Repr (Enumeration enum) = ()
+   cons = MultiValue.consUnit
+   undef = MultiValue.undefUnit
+   zero = MultiValue.zeroUnit
+   phi = MultiValue.phiUnit
+   addPhi = MultiValue.addPhiUnit
+
+type instance MultiValue.Decomposed f (Enumeration enum) = Enumeration enum
+type instance MultiValue.PatternTuple (Enumeration enum) = Enumeration enum
+
+instance
+      (Enum enum, Bounded enum) => MultiValue.Compose (Enumeration enum) where
+   type Composed (Enumeration enum) = Enumeration enum
+   compose = MultiValue.cons
+
+instance MultiValue.Decompose (Enumeration enum) where
+   decompose Enumeration _ = Enumeration
+
+
+instance (Enum enum, Bounded enum) => Expr.Compose (Enumeration enum) where
+   type Composed (Enumeration enum) = Enumeration enum
+   compose = Expr.cons
+
+instance Expr.Decompose (Enumeration enum) where
+   decompose Enumeration _ = Enumeration
diff --git a/src/Data/Array/Knead/Symbolic.hs b/src/Data/Array/Knead/Symbolic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic.hs
@@ -0,0 +1,94 @@
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic (
+   Core.Array,
+   Core.C(..),
+   Exp,
+   fix,
+   shape,
+   (Core.!),
+   Core.the,
+   Core.fromScalar,
+   Core.fill,
+   gather,
+   backpermute,
+   Core.backpermute2,
+   Core.id,
+   Core.map,
+   Core.mapWithIndex,
+   zipWith,
+   zipWith3,
+   zipWith4,
+   zip,
+   zip3,
+   zip4,
+   Core.fold1,
+   Core.fold1All,
+   Core.findAll,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.ShapeDependent as ShapeDep
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import Data.Array.Knead.Symbolic.Private (Array, shape, gather, )
+
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified LLVM.Extra.Multi.Value as MultiValue
+
+import Data.Function.HT (Id)
+
+import Prelude hiding (zipWith, zipWith3, zip, zip3, replicate, )
+
+
+fix :: Id (Array sh a)
+fix = id
+
+backpermute ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    MultiValue.C a) =>
+   Exp sh1 ->
+   (Exp ix1 -> Exp ix0) ->
+   Array sh0 a ->
+   Array sh1 a
+backpermute sh1 f = gather (Core.map f (Core.id sh1))
+
+zipWith ::
+   (Core.C array, Shape.C sh) =>
+   (Exp a -> Exp b -> Exp c) ->
+   array sh a -> array sh b -> array sh c
+zipWith = ShapeDep.backpermute2 Shape.intersect id id
+
+zipWith3 ::
+   (Core.C array, Shape.C sh) =>
+   (Exp a -> Exp b -> Exp c -> Exp d) ->
+   array sh a -> array sh b -> array sh c -> array sh d
+zipWith3 f a b c =
+   zipWith (\ab -> uncurry f (Expr.unzip ab)) (zipWith Expr.zip a b) c
+
+zipWith4 ::
+   (Core.C array, Shape.C sh) =>
+   (Exp a -> Exp b -> Exp c -> Exp d -> Exp e) ->
+   array sh a -> array sh b -> array sh c -> array sh d -> array sh e
+zipWith4 f a b c d =
+   zipWith3 (\ab -> uncurry f (Expr.unzip ab)) (zipWith Expr.zip a b) c d
+
+
+zip ::
+   (Core.C array, Shape.C sh) =>
+   array sh a -> array sh b -> array sh (a,b)
+zip = zipWith (Expr.lift2 MultiValue.zip)
+
+zip3 ::
+   (Core.C array, Shape.C sh) =>
+   array sh a -> array sh b -> array sh c -> array sh (a,b,c)
+zip3 = zipWith3 (Expr.lift3 MultiValue.zip3)
+
+zip4 ::
+   (Core.C array, Shape.C sh) =>
+   array sh a -> array sh b -> array sh c -> array sh d ->
+   array sh (a,b,c,d)
+zip4 = zipWith4 (Expr.lift4 MultiValue.zip4)
diff --git a/src/Data/Array/Knead/Symbolic/Fold.hs b/src/Data/Array/Knead/Symbolic/Fold.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Fold.hs
@@ -0,0 +1,98 @@
+{- |
+Reduce selected dimensions.
+Alternatively you may reorder dimensions with 'ShapeDep.backpermute'
+and fold once along multiple dimensions.
+-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic.Fold (
+   T,
+   Cubic,
+   apply,
+   passAny,
+   pass,
+   fold,
+   (Core.$:.),
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import Data.Array.Knead.Symbolic.Private (Array(Array), Code, Val, )
+
+import qualified Data.Array.Knead.Shape.Cubic.Int as Index
+import qualified Data.Array.Knead.Shape.Cubic as Cubic
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Shape.Cubic ((#:.), (:.)((:.)), )
+
+import LLVM.DSL.Expression (Exp, unExp)
+
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import LLVM.Extra.Multi.Value (atom, )
+
+import qualified Type.Data.Num.Unary as Unary
+
+import Prelude hiding (zipWith, zipWith3, zip, zip3, replicate, )
+
+
+data T sh0 sh1 a =
+   forall ix0 ix1.
+   (Shape.Index sh0 ~ ix0, Shape.Index sh1 ~ ix1) =>
+   Cons
+      (Exp sh0 -> Exp sh1)
+      (forall r. Val sh0 -> (Val ix0 -> Code r a) -> (Val ix1 -> Code r a))
+
+
+apply ::
+   (Core.C array, Shape.C sh0, Shape.C sh1, MultiValue.C a) =>
+   T sh0 sh1 a ->
+   array sh0 a ->
+   array sh1 a
+apply (Cons fsh reduce) =
+   Core.lift1 $ \(Array sh code) ->
+      Array (fsh sh) (\ix -> do sh0 <- unExp sh; reduce sh0 code ix)
+
+
+type Cubic rank0 rank1 = T (Cubic.Shape rank0) (Cubic.Shape rank1)
+
+passAny :: Cubic rank rank a
+passAny = Cons id (const id)
+
+pass ::
+   (Unary.Natural rank0, Unary.Natural rank1, MultiValue.C a) =>
+   Cubic rank0 rank1 a ->
+   Cubic (Unary.Succ rank0) (Unary.Succ rank1) a
+pass (Cons fsh reduce) =
+   Cons
+      (Expr.modify (atom:.atom) $ \(sh:.s) -> fsh sh :. s)
+      (\sh code ->
+       Cubic.switchR $ \jx j ->
+          reduce (Cubic.tail sh) (\kx -> code (kx #:. j)) jx)
+
+
+fold1CodeLinear ::
+   (Unary.Natural rank, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Exp Index.Int ->
+   (Val (Cubic.Index (Unary.Succ rank)) -> Code r a) ->
+   (Val (Cubic.Index rank) -> Code r a)
+fold1CodeLinear f nc code ix =
+   Core.fold1Code f
+      (Expr.lift1 (MultiValue.compose . Shape.ZeroBased) $ Index.decons nc)
+      (\j -> code (ix #:. Index.cons j))
+
+fold ::
+   (Unary.Natural rank0, Unary.Natural rank1, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Cubic rank0 rank1 a ->
+   Cubic (Unary.Succ rank0) rank1 a
+fold f (Cons fsh reduce) =
+   Cons
+      (fsh . Cubic.tail)
+      (\sh code jx ->
+          reduce (Cubic.tail sh)
+             (fold1CodeLinear f (Expr.lift0 (Cubic.head sh)) code) jx)
+
+
+instance Core.Process (T sh0 sh1 a) where
diff --git a/src/Data/Array/Knead/Symbolic/Physical.hs b/src/Data/Array/Knead/Symbolic/Physical.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Physical.hs
@@ -0,0 +1,195 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+module Data.Array.Knead.Symbolic.Physical (
+   Array,
+   shape,
+   toList,
+   fromList,
+   vectorFromList,
+   with,
+   render,
+   scanl1,
+   mapAccumLSimple,
+   scatter,
+   scatterMaybe,
+   permute,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.PhysicalPrivate as Priv
+import qualified Data.Array.Knead.Symbolic.Private as Sym
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Symbolic.PhysicalPrivate (MarshalPtr)
+import Data.Array.Knead.Code (getElementPtr)
+
+import qualified LLVM.DSL.Execution as Code
+import LLVM.DSL.Expression (Exp, unExp)
+
+import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as MutArray
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+import qualified Data.Array.Comfort.Shape as ComfortShape
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Memory as Memory
+import qualified LLVM.Extra.Maybe as Maybe
+
+import qualified LLVM.Core as LLVM
+
+import Foreign.Storable (Storable, )
+import Foreign.ForeignPtr (withForeignPtr, mallocForeignPtrArray, )
+import Foreign.Ptr (FunPtr, Ptr, )
+
+import Control.Monad.HT (void, (<=<), )
+import Control.Applicative (liftA2, (<$>), )
+
+import Prelude2010 hiding (scanl1)
+import Prelude ()
+
+
+shape :: Array sh a -> sh
+shape = Array.shape
+
+toList ::
+   (Shape.C sh, Storable a) =>
+   Array sh a -> IO [a]
+toList = MutArray.toList <=< MutArray.unsafeThaw
+
+fromList ::
+   (Shape.C sh, Storable a) =>
+   sh -> [a] -> IO (Array sh a)
+fromList sh = MutArray.unsafeFreeze <=< MutArray.fromList sh
+
+vectorFromList ::
+   (Num n, Storable a) =>
+   [a] -> IO (Array (ComfortShape.ZeroBased n) a)
+vectorFromList xs =
+   Array.mapShape (\(Shape.ZeroBased n) -> Shape.ZeroBased $ fromIntegral n) <$>
+   (MutArray.unsafeFreeze =<< MutArray.vectorFromList xs)
+
+
+{- |
+The symbolic array is only valid inside the enclosed action.
+-}
+with ::
+   (Shape.C sh, Storable.C a) =>
+   (Sym.Array sh a -> IO b) ->
+   Array sh a -> IO b
+with f (Array sh fptr) =
+   withForeignPtr fptr $ \ptr ->
+      f $
+      Sym.Array
+         (Shape.value sh)
+         (\ix ->
+            Storable.load =<<
+               getElementPtr (Shape.value sh) (LLVM.valueOf ptr) ix)
+
+
+type Importer f = FunPtr f -> f
+
+foreign import ccall safe "dynamic" callShaper ::
+   Importer (LLVM.Ptr sh -> IO Shape.Size)
+
+foreign import ccall safe "dynamic" callRenderer ::
+   Importer (LLVM.Ptr sh -> Ptr a -> IO ())
+
+
+materialize ::
+   (Shape.C sh, Marshal.C sh, Storable.C a) =>
+   String ->
+   Exp sh ->
+   (LLVM.Value (MarshalPtr sh) ->
+    LLVM.Value (Ptr a) -> LLVM.CodeGenFunction () ()) ->
+   IO (Array sh a)
+materialize name esh code =
+   Marshal.alloca $ \lshptr -> do
+      (fsh, farr) <-
+         Code.compile name $
+         liftA2 (,)
+            (Code.createFunction callShaper "shape" $ \ptr -> do
+               sh <- unExp esh
+               Memory.store sh ptr
+               Shape.size sh)
+            (Code.createFunction callRenderer "fill" code)
+      n <- fsh lshptr
+      fptr <- mallocForeignPtrArray (fromIntegral n)
+      withForeignPtr fptr $ farr lshptr
+      sh <- Marshal.peek lshptr
+      return (Array sh fptr)
+
+render ::
+   (Shape.C sh, Marshal.C sh, Storable.C a) =>
+   Sym.Array sh a -> IO (Array sh a)
+render (Sym.Array esh code) =
+   materialize "render" esh $ \sptr ptr -> do
+      let step ix p = flip Storable.storeNext p =<< code ix
+      sh <- Shape.load esh sptr
+      void $ Shape.loop step sh ptr
+
+scanl1 ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    Storable.C a, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array (sh, n) a -> IO (Array (sh, n) a)
+scanl1 f (Sym.Array esh code) =
+   materialize "scanl1" esh $ \sptr ptr -> do
+      (sh, n) <- MultiValue.unzip <$> Shape.load esh sptr
+      let step ix ptrStart =
+             fmap fst $
+             (\body -> Shape.loop body n (ptrStart, Maybe.nothing)) $
+                   \k0 (ptr0, macc0) -> do
+                a <- code $ MultiValue.zip ix k0
+                acc1 <- Maybe.run macc0 (return a) (flip (Expr.unliftM2 f) a)
+                ptr1 <- Storable.storeNext acc1 ptr0
+                return (ptr1, Maybe.just acc1)
+      void $ Shape.loop step sh ptr
+
+mapAccumLSimple ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc, Storable.C x, Storable.C y) =>
+   (Exp acc -> Exp x -> Exp (acc,y)) ->
+   Sym.Array sh acc -> Sym.Array (sh, n) x -> IO (Array (sh, n) y)
+mapAccumLSimple f arrInit arrData =
+   materialize "mapAccumLSimple" (Sym.shape arrData) $
+      Priv.mapAccumLSimple f arrInit arrData
+
+scatterMaybe ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1, Marshal.C sh1,
+    Storable.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array sh1 a ->
+   Sym.Array sh0 (Maybe (ix1, a)) -> IO (Array sh1 a)
+scatterMaybe accum arrInit arrMap =
+   materialize "scatterMaybe" (Sym.shape arrInit) $
+      Priv.scatterMaybe accum arrInit arrMap
+
+scatter ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1, Marshal.C sh1,
+    Storable.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array sh1 a ->
+   Sym.Array sh0 (ix1, a) -> IO (Array sh1 a)
+scatter accum arrInit arrMap =
+   materialize "scatter" (Sym.shape arrInit) $
+      Priv.scatter accum arrInit arrMap
+
+permute ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1, Marshal.C sh1,
+    Storable.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array sh1 a ->
+   (Exp ix0 -> Exp ix1) ->
+   Sym.Array sh0 a ->
+   IO (Array sh1 a)
+permute accum deflt ixmap input =
+   scatter accum deflt
+      (Sym.mapWithIndex (Expr.lift2 MultiValue.zip . ixmap) input)
diff --git a/src/Data/Array/Knead/Symbolic/PhysicalParametric.hs b/src/Data/Array/Knead/Symbolic/PhysicalParametric.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/PhysicalParametric.hs
@@ -0,0 +1,455 @@
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+module Data.Array.Knead.Symbolic.PhysicalParametric (
+   the,
+   theMarshal,
+   render,
+   MapFilter(..),
+   mapFilter,
+   FilterOuter(..),
+   filterOuter,
+   Scatter(..),
+   scatter,
+   ScatterMaybe(..),
+   scatterMaybe,
+   MapAccumLSimple(..),
+   mapAccumLSimple,
+   MapAccumLSequence(..),
+   mapAccumLSequence,
+   MapAccumL(..),
+   mapAccumL,
+   FoldOuterL(..),
+   foldOuterL,
+   AddDimension(..),
+   addDimension,
+
+   Parametric,
+   Rendered,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.PhysicalPrivate as Priv
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Symbolic.PhysicalPrivate (MarshalPtr)
+
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+
+import qualified LLVM.DSL.Execution as Code
+import LLVM.DSL.Expression (Exp(Exp), unExp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Memory as Memory
+import qualified LLVM.Extra.Arithmetic as A
+
+import qualified LLVM.Core as LLVM
+
+import Foreign.Marshal.Array (allocaArray, )
+import Foreign.Marshal.Alloc (alloca, )
+import Foreign.Storable (Storable, peek, peekElemOff, )
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, mallocForeignPtrArray, )
+import Foreign.Ptr (FunPtr, Ptr, )
+
+import Control.Exception (finally)
+import Control.Monad.HT (void, )
+import Control.Applicative (liftA2, )
+
+
+mallocArray :: (Storable a) => Shape.Size -> IO (ForeignPtr a)
+mallocArray = mallocForeignPtrArray . fromIntegral
+
+
+type Importer f = FunPtr f -> f
+
+
+
+
+type Parametric p a = Exp p -> a
+type Rendered p a = IO (p, IO ()) -> IO a
+
+withManagedParam :: Monad m => (p -> IO a) -> m (Rendered p a)
+withManagedParam act =
+   return $ \create -> do
+      (param, final) <- create
+      finally (act param) final
+
+
+
+foreign import ccall safe "dynamic" callThe ::
+   Importer (LLVM.Ptr param -> Ptr a -> IO ())
+
+the ::
+   (Marshal.C p, Shape.Scalar z, Storable.C a) =>
+   Parametric p (Core.Array z a) -> IO (Rendered p a)
+the arr = do
+   func <-
+      Code.compile "the" $
+      Code.createFunction callThe "eval" $
+      \paramPtr resultPtr -> do
+         case arr $ Exp (Memory.load paramPtr) of
+            Core.Array z code ->
+               code (Shape.zeroIndex z) >>=
+               flip Storable.store resultPtr
+   withManagedParam $ \param ->
+      Marshal.with param $ \pptr ->
+      alloca $ \aptr -> func pptr aptr >> peek aptr
+
+foreign import ccall safe "dynamic" callTheMarshal ::
+   Importer (LLVM.Ptr param -> LLVM.Ptr a -> IO ())
+
+theMarshal ::
+   (Marshal.C p, Shape.Scalar z, Marshal.C a) =>
+   Parametric p (Core.Array z a) -> IO (Rendered p a)
+theMarshal arr = do
+   func <-
+      Code.compile "the-marshal" $
+      Code.createFunction callTheMarshal "eval" $
+      \paramPtr resultPtr -> do
+         case arr $ Exp (Memory.load paramPtr) of
+            Core.Array z code ->
+               code (Shape.zeroIndex z) >>=
+               flip Memory.store resultPtr
+   withManagedParam $ \param ->
+      Marshal.with param $ \pptr ->
+      Marshal.alloca $ \aptr ->
+         func pptr aptr >>
+         Marshal.peek aptr
+
+
+
+
+foreign import ccall safe "dynamic" callShaper ::
+   Importer (LLVM.Ptr param -> LLVM.Ptr shape -> IO Shape.Size)
+
+foreign import ccall safe "dynamic" callFill ::
+   Importer (LLVM.Ptr param -> LLVM.Ptr shape -> Ptr a -> IO ())
+
+
+{-
+Attention:
+The 'fill' function may alter the shape.
+An example is 'mapFilter'.
+-}
+materialize ::
+   (Shape.C sh, Marshal.C sh, Marshal.C p, Storable.C a) =>
+   String ->
+   (core -> Exp sh) ->
+   (core ->
+    LLVM.Value (MarshalPtr sh) -> LLVM.Value (Ptr a) ->
+    LLVM.CodeGenFunction () ()) ->
+   Parametric p core -> IO (Rendered p (Array sh a))
+materialize name shape fill core = do
+   (fsh, farr) <-
+      Code.compile name $
+      liftA2 (,)
+         (Code.createFunction callShaper "shape" $
+          \paramPtr resultPtr -> do
+            sh <- unExp $ shape $ core $ Exp (Memory.load paramPtr)
+            Memory.store sh resultPtr
+            Shape.size sh)
+         (Code.createFunction callFill "fill" $
+          \paramPtr shapePtr bufferPtr ->
+            fill (core $ Exp (Memory.load paramPtr)) shapePtr bufferPtr)
+
+   withManagedParam $ \param ->
+      Marshal.alloca $ \shptr ->
+      Marshal.with param $ \paramPtr -> do
+         fptr <- mallocArray =<< fsh paramPtr shptr
+         withForeignPtr fptr $ farr paramPtr shptr
+         sh <- Marshal.peek shptr
+         return (Array sh fptr)
+
+
+foreign import ccall safe "dynamic" callFillExpArray ::
+   Importer (LLVM.Ptr param -> Ptr final -> LLVM.Ptr shape -> Ptr a -> IO ())
+
+
+materializeExpArray ::
+   (Shape.C sh, Marshal.C sh, Marshal.C p, Storable.C a, Storable.C b) =>
+   String ->
+   (core -> Exp sh) ->
+   (core ->
+    LLVM.Value (Ptr b) ->
+    LLVM.Value (MarshalPtr sh) ->
+    LLVM.Value (Ptr a) ->
+    LLVM.CodeGenFunction () ()) ->
+   Parametric p core -> IO (Rendered p (b, Array sh a))
+materializeExpArray name shape fill core = do
+   (fsh, farr) <-
+      Code.compile name $
+      liftA2 (,)
+         (Code.createFunction callShaper "shape" $
+          \paramPtr resultPtr -> do
+            sh <- unExp $ shape $ core $ Exp (Memory.load paramPtr)
+            Memory.store sh resultPtr
+            Shape.size sh)
+         (Code.createFunction callFillExpArray "fill" $
+          \paramPtr finalPtr shapePtr bufferPtr ->
+            fill
+               (core $ Exp (Memory.load paramPtr))
+               finalPtr shapePtr bufferPtr)
+
+   withManagedParam $ \param ->
+      Marshal.alloca $ \shptr ->
+      alloca $ \finalPtr ->
+      Marshal.with param $ \paramPtr -> do
+         fptr <- mallocArray =<< fsh paramPtr shptr
+         withForeignPtr fptr $ farr paramPtr finalPtr shptr
+         sh <- Marshal.peek shptr
+         final <- peek finalPtr
+         return (final, Array sh fptr)
+
+
+foreign import ccall safe "dynamic" callShaper2 ::
+   Importer
+      (LLVM.Ptr param ->
+       LLVM.Ptr shapeA -> LLVM.Ptr shapeB -> Ptr Shape.Size -> IO ())
+
+foreign import ccall safe "dynamic" callFill2 ::
+   Importer
+      (LLVM.Ptr param ->
+       LLVM.Ptr shapeA -> Ptr a -> LLVM.Ptr shapeB -> Ptr b -> IO ())
+
+
+materialize2 ::
+   (Shape.C sha, Marshal.C sha,
+    Shape.C shb, Marshal.C shb,
+    Marshal.C p, Storable.C a, Storable.C b) =>
+   String ->
+   (core -> Exp (sha,shb)) ->
+   (core ->
+    (LLVM.Value (MarshalPtr sha), LLVM.Value (Ptr a)) ->
+    (LLVM.Value (MarshalPtr shb), LLVM.Value (Ptr b)) ->
+    LLVM.CodeGenFunction () ()) ->
+   Parametric p core -> IO (Rendered p (Array sha a, Array shb b))
+materialize2 name shape fill core = do
+   (fsh, farr) <-
+      Code.compile name $
+      liftA2 (,)
+         (Code.createFunction callShaper2 "shape" $
+          \paramPtr shapeAPtr shapeBPtr sizesPtr -> do
+            (sha,shb) <-
+               fmap MultiValue.unzip $ unExp $
+               shape $ core $ Exp (Memory.load paramPtr)
+            Memory.store sha shapeAPtr
+            Memory.store shb shapeBPtr
+            sizeAPtr <- LLVM.bitcast sizesPtr
+            flip LLVM.store sizeAPtr =<< Shape.size sha
+            sizeBPtr <- A.advanceArrayElementPtr sizeAPtr
+            flip LLVM.store sizeBPtr =<< Shape.size shb)
+         (Code.createFunction callFill2 "fill" $
+          \paramPtr shapeAPtr bufferAPtr shapeBPtr bufferBPtr ->
+            fill
+               (core $ Exp (Memory.load paramPtr))
+               (shapeAPtr, bufferAPtr) (shapeBPtr, bufferBPtr))
+
+   withManagedParam $ \param ->
+      Marshal.alloca $ \shaPtr ->
+      Marshal.alloca $ \shbPtr ->
+      allocaArray 2 $ \sizesPtr ->
+      Marshal.with param $ \paramPtr -> do
+         fsh paramPtr shaPtr shbPtr sizesPtr
+         afptr <- mallocArray =<< peekElemOff sizesPtr 0
+         bfptr <- mallocArray =<< peekElemOff sizesPtr 1
+         withForeignPtr afptr $ \aptr ->
+            withForeignPtr bfptr $ \bptr ->
+            farr paramPtr shaPtr aptr shbPtr bptr
+         sha <- Marshal.peek shaPtr
+         shb <- Marshal.peek shbPtr
+         return (Array sha afptr, Array shb bfptr)
+
+
+render ::
+   (Shape.C sh, Shape.Index sh ~ ix, Marshal.C sh,
+    Marshal.C p, Storable.C a) =>
+   Parametric p (Core.Array sh a) -> IO (Rendered p (Array sh a))
+render =
+   materialize "render" Core.shape
+      (\(Core.Array esh code) shapePtr bufferPtr -> do
+         let step ix p = flip Storable.storeNext p =<< code ix
+         sh <- Shape.load esh shapePtr
+         void $ Shape.loop step sh bufferPtr)
+
+
+data Scatter sh0 sh1 a =
+   Scatter {
+      scatterAccum :: Exp a -> Exp a -> Exp a,
+      scatterInit :: Core.Array sh1 a,
+      scatterMap :: Core.Array sh0 (Shape.Index sh1, a)
+   }
+
+scatter ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1, Marshal.C sh1,
+    Marshal.C p, Storable.C a) =>
+   Parametric p (Scatter sh0 sh1 a) -> IO (Rendered p (Array sh1 a))
+scatter =
+   materialize "scatter"
+      (Core.shape . scatterInit)
+      (\(Scatter accum arrInit arrMap) ->
+         Priv.scatter accum arrInit arrMap)
+
+
+
+data ScatterMaybe sh0 sh1 a =
+   ScatterMaybe {
+      scatterMaybeAccum :: Exp a -> Exp a -> Exp a,
+      scatterMaybeInit :: Core.Array sh1 a,
+      scatterMaybeMap :: Core.Array sh0 (Maybe (Shape.Index sh1, a))
+   }
+
+scatterMaybe ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1, Marshal.C sh1,
+    Marshal.C p, Storable.C a) =>
+   Parametric p (ScatterMaybe sh0 sh1 a) -> IO (Rendered p (Array sh1 a))
+scatterMaybe =
+   materialize "scatterMaybe"
+      (Core.shape . scatterMaybeInit)
+      (\(ScatterMaybe accum arrInit arrMap) ->
+         Priv.scatterMaybe accum arrInit arrMap)
+
+
+data MapAccumLSimple sh n acc a b =
+   MapAccumLSimple {
+      mapAccumLSimpleAccum :: Exp acc -> Exp a -> Exp (acc,b),
+      mapAccumLSimpleInit :: Core.Array sh acc,
+      mapAccumLSimpleArray :: Core.Array (sh, n) a
+   }
+
+mapAccumLSimple ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc, Marshal.C p, Storable.C a, Storable.C b) =>
+   Parametric p (MapAccumLSimple sh n acc a b) ->
+   IO (Rendered p (Array (sh,n) b))
+mapAccumLSimple =
+   materialize "mapAccumLSimple"
+      (Core.shape . mapAccumLSimpleArray)
+      (\(MapAccumLSimple f arrInit arrData) ->
+         Priv.mapAccumLSimple f arrInit arrData)
+
+
+data MapAccumLSequence n acc final a b =
+   MapAccumLSequence {
+      mapAccumLSequenceAccum :: Exp acc -> Exp a -> Exp (acc,b),
+      mapAccumLSequenceFinal :: Exp acc -> Exp final,
+      mapAccumLSequenceInit :: Exp acc,
+      mapAccumLSequenceArray :: Core.Array n a
+   }
+
+-- FIXME: check correct size of array of initial values
+mapAccumLSequence ::
+   (Shape.C n, Marshal.C n, MultiValue.C acc, Storable.C final,
+    Marshal.C p, Storable.C a, Storable.C b) =>
+   Parametric p (MapAccumLSequence n acc final a b) ->
+   IO (Rendered p (final, Array n b))
+mapAccumLSequence =
+   materializeExpArray "mapAccumLSequence"
+      (Core.shape . mapAccumLSequenceArray)
+      (\(MapAccumLSequence f final expInit arr) ->
+         Priv.mapAccumLSequence f final expInit arr)
+
+
+data MapAccumL sh n acc final a b =
+   MapAccumL {
+      mapAccumLAccum :: Exp acc -> Exp a -> Exp (acc,b),
+      mapAccumLFinal :: Exp acc -> Exp final,
+      mapAccumLInit :: Core.Array sh acc,
+      mapAccumLArray :: Core.Array (sh, n) a
+   }
+
+-- FIXME: check correct size of array of initial values
+mapAccumL ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc, Storable.C final,
+    Marshal.C p, Storable.C a, Storable.C b) =>
+   Parametric p (MapAccumL sh n acc final a b) ->
+   IO (Rendered p (Array sh final, Array (sh,n) b))
+mapAccumL =
+   materialize2 "mapAccumL"
+      (\core ->
+         Expr.zip
+            (Core.shape $ mapAccumLInit core)
+            (Core.shape $ mapAccumLArray core))
+      (\(MapAccumL f final arrInit arrData) ->
+         Priv.mapAccumL f final arrInit arrData)
+
+
+data FoldOuterL n sh a b =
+   FoldOuterL {
+      foldOuterLAccum :: Exp a -> Exp b -> Exp a,
+      foldOuterLInit :: Core.Array sh a,
+      foldOuterLArray :: Core.Array (n,sh) b
+   }
+
+-- FIXME: check correct size of array of initial values
+foldOuterL ::
+   (Shape.C n, Marshal.C n,
+    Shape.C sh, Marshal.C sh,
+    Marshal.C p, Storable.C a) =>
+   Parametric p (FoldOuterL n sh a b) -> IO (Rendered p (Array sh a))
+foldOuterL =
+   materialize "foldOuterL"
+      (Core.shape . foldOuterLInit)
+      (\(FoldOuterL f arrInit arrData) -> Priv.foldOuterL f arrInit arrData)
+
+
+data MapFilter n a b =
+   MapFilter {
+      mapFilterMap :: Exp a -> Exp b,
+      mapFilterPredicate :: Exp a -> Exp Bool,
+      mapFilterArray :: Core.Array n a
+   }
+
+mapFilter ::
+   (Shape.Sequence n, Marshal.C n, Marshal.C p, Storable.C b) =>
+   Parametric p (MapFilter n a b) -> IO (Rendered p (Array n b))
+mapFilter =
+   materialize "mapFilter"
+      (Core.shape . mapFilterArray)
+      (\(MapFilter f p arr) shapePtr bufferPtr ->
+         flip Memory.store shapePtr
+            =<< Priv.mapFilter f p arr shapePtr bufferPtr)
+
+
+data FilterOuter n sh a =
+   FilterOuter {
+      filterOuterPredicate :: Core.Array n Bool,
+      filterOuterArray :: Core.Array (n,sh) a
+   }
+
+-- FIXME: check correct size of row selection array
+filterOuter ::
+   (Shape.Sequence n, Marshal.C n,
+    Shape.C sh, Marshal.C sh,
+    Marshal.C p, Storable.C a) =>
+   Parametric p (FilterOuter n sh a) -> IO (Rendered p (Array (n,sh) a))
+filterOuter =
+   materialize "filterOuter"
+      (Core.shape . filterOuterArray)
+      (\(FilterOuter p arr) shapePtr bufferPtr ->
+         flip Memory.store shapePtr
+            =<< Priv.filterOuter p arr shapePtr bufferPtr)
+
+
+data AddDimension sh n a b =
+   AddDimension {
+      addDimensionSize :: Exp n,
+      addDimensionSelect :: Exp (Shape.Index n) -> Exp a -> Exp b,
+      addDimensionArray :: Core.Array sh a
+   }
+
+addDimension ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    Marshal.C p, Storable.C b) =>
+   Parametric p (AddDimension sh n a b) -> IO (Rendered p (Array (sh,n) b))
+addDimension =
+   materialize "addDimension"
+      (\r -> Expr.zip (Core.shape (addDimensionArray r)) (addDimensionSize r))
+      (\(AddDimension n select arr) -> Priv.addDimension n select arr)
diff --git a/src/Data/Array/Knead/Symbolic/PhysicalPrivate.hs b/src/Data/Array/Knead/Symbolic/PhysicalPrivate.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/PhysicalPrivate.hs
@@ -0,0 +1,259 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic.PhysicalPrivate where
+
+import qualified Data.Array.Knead.Symbolic.Private as Sym
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Code (getElementPtr)
+
+import LLVM.DSL.Expression (Exp, unExp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Control as C
+
+import qualified LLVM.Core as LLVM
+
+import Foreign.Ptr (Ptr, )
+
+import qualified Control.Applicative.HT as App
+import Control.Monad.HT (void, )
+import Control.Applicative ((<$>), )
+
+import Data.Tuple.HT (mapSnd, )
+
+import Prelude2010
+import Prelude ()
+
+
+
+type MarshalPtr a = LLVM.Ptr (Marshal.Struct a)
+
+writeArray ::
+   (Shape.C sh, Shape.Index sh ~ ix, Storable.C a) =>
+   MultiValue.T sh ->
+   (MultiValue.T ix -> LLVM.CodeGenFunction r (MultiValue.T a)) ->
+   LLVM.Value (Ptr a) ->
+   LLVM.CodeGenFunction r (LLVM.Value (Ptr a))
+writeArray sh code ptr = do
+   let clear ix p = flip Storable.storeNext p =<< code ix
+   Shape.loop clear sh ptr
+
+
+mapAccumLLoop ::
+   (MultiValue.C acc, Storable.C b,
+    Shape.C sh, Shape.Index sh ~ ix) =>
+   (MultiValue.T ix -> LLVM.CodeGenFunction r (MultiValue.T a)) ->
+   (Exp acc -> Exp a -> Exp (acc, b)) ->
+   MultiValue.T sh ->
+   LLVM.Value (Ptr b) -> MultiValue.T acc ->
+   LLVM.CodeGenFunction r (LLVM.Value (Ptr b), MultiValue.T acc)
+mapAccumLLoop code f n yPtr accInit = do
+   let step k0 (ptr0, acc0) = do
+         x <- code k0
+         (acc1,y) <- MultiValue.unzip <$> Expr.unliftM2 f acc0 x
+         ptr1 <- Storable.storeNext y ptr0
+         return (ptr1, acc1)
+   Shape.loop step n (yPtr, accInit)
+
+mapAccumLSimple ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc,
+    Storable.C x,
+    Storable.C y) =>
+   (Exp acc -> Exp x -> Exp (acc,y)) ->
+   Sym.Array sh acc -> Sym.Array (sh, n) x ->
+   LLVM.Value (MarshalPtr (sh,n)) ->
+   LLVM.Value (Ptr y) ->
+   LLVM.CodeGenFunction r ()
+mapAccumLSimple f (Sym.Array _ initCode) (Sym.Array esh code) sptr ptr = do
+   (sh, n) <- MultiValue.unzip <$> Shape.load esh sptr
+   let step ix ptrStart = do
+         accInit <- initCode ix
+         fst <$> mapAccumLLoop (code . MultiValue.zip ix) f n ptrStart accInit
+   void $ Shape.loop step sh ptr
+
+mapAccumLSequence ::
+   (Shape.C n, Marshal.C n,
+    MultiValue.C acc, Storable.C final,
+    Storable.C x,
+    Storable.C y) =>
+   (Exp acc -> Exp x -> Exp (acc,y)) ->
+   (Exp acc -> Exp final) ->
+   Exp acc -> Sym.Array n x ->
+   LLVM.Value (Ptr final) ->
+   LLVM.Value (MarshalPtr n) ->
+   LLVM.Value (Ptr y) ->
+   LLVM.CodeGenFunction r ()
+mapAccumLSequence f final initExp (Sym.Array esh code) accPtr sptr yPtr = do
+   n <- Shape.load esh sptr
+   accInit <- unExp initExp
+   accExit <- snd <$> mapAccumLLoop code f n yPtr accInit
+   flip Storable.store accPtr =<< Expr.unliftM1 final accExit
+
+mapAccumL ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc, Storable.C final,
+    Storable.C x,
+    Storable.C y) =>
+   (Exp acc -> Exp x -> Exp (acc,y)) ->
+   (Exp acc -> Exp final) ->
+   Sym.Array sh acc -> Sym.Array (sh, n) x ->
+   (LLVM.Value (MarshalPtr sh), LLVM.Value (Ptr final)) ->
+   (LLVM.Value (MarshalPtr (sh,n)), LLVM.Value (Ptr y)) ->
+   LLVM.CodeGenFunction r ()
+mapAccumL f final (Sym.Array _ initCode) (Sym.Array esh code)
+      (_, accPtr) (sptr, yPtr) = do
+   (sh, n) <- MultiValue.unzip <$> Shape.load esh sptr
+   let step ix (accPtr0, yPtrStart) = do
+         accInit <- initCode ix
+         (ptrStop, accExit) <-
+            mapAccumLLoop (code . MultiValue.zip ix) f n yPtrStart accInit
+         accPtr1 <-
+            flip Storable.storeNext accPtr0
+               =<< Expr.unliftM1 final accExit
+         return (accPtr1, ptrStop)
+   void $ Shape.loop step sh (accPtr,yPtr)
+
+foldOuterL ::
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    Storable.C a) =>
+   (Exp a -> Exp b -> Exp a) ->
+   Sym.Array sh a -> Sym.Array (n,sh) b ->
+   LLVM.Value (MarshalPtr sh) ->
+   LLVM.Value (Ptr a) ->
+   LLVM.CodeGenFunction r ()
+foldOuterL f (Sym.Array _ initCode) (Sym.Array esh code) sptr ptr = do
+   sh <- Shape.load (Expr.snd esh) sptr
+   n <- MultiValue.fst <$> unExp esh
+   void $ writeArray sh initCode ptr
+
+   let step k ix ptr0 = do
+         b <- code $ MultiValue.zip k ix
+         a0 <- Storable.load ptr0
+         a1 <- Expr.unliftM2 f a0 b
+         Storable.storeNext a1 ptr0
+   void $ Shape.loop (\k () -> void $ Shape.loop (step k) sh ptr) n ()
+
+{- |
+We need a scalar Shape type @n@.
+Scalar Shape types could be distinguished from other Shape types
+by the fact that you can convert any Index into a Shape.
+-}
+mapFilter ::
+   (Shape.Sequence n, Marshal.C n,
+    Storable.C b) =>
+   (Exp a -> Exp b) ->
+   (Exp a -> Exp Bool) ->
+   Sym.Array n a ->
+   LLVM.Value (MarshalPtr n) ->
+   LLVM.Value (Ptr b) ->
+   LLVM.CodeGenFunction r (MultiValue.T n)
+mapFilter f p (Sym.Array esh code) sptr ptr = do
+   n <- Shape.load esh sptr
+   let step ix (dstPtr,dstIx) = do
+         a <- code ix
+         MultiValue.Cons c <- Expr.unliftM1 p a
+         C.ifThen c (dstPtr,dstIx)
+            (App.lift2 (,)
+               (flip Storable.storeNext dstPtr =<< Expr.unliftM1 f a)
+               (MultiValue.inc dstIx))
+   Shape.sequenceShapeFromIndex . snd
+      =<< Shape.loop step n (ptr, MultiValue.zero)
+
+filterOuter ::
+   (Shape.Sequence n, Marshal.C n,
+    Shape.C sh, Marshal.C sh,
+    Storable.C a) =>
+   Sym.Array n Bool ->
+   Sym.Array (n,sh) a ->
+   LLVM.Value (MarshalPtr (n,sh)) ->
+   LLVM.Value (Ptr a) ->
+   LLVM.CodeGenFunction r (MultiValue.T (n,sh))
+filterOuter (Sym.Array _eish selectCode) (Sym.Array esh code) sptr ptr = do
+   (n,sh) <- MultiValue.unzip <$> Shape.load esh sptr
+   let step k (dstPtr0,dstK) = do
+         MultiValue.Cons c <- selectCode k
+         C.ifThen c (dstPtr0,dstK)
+            (do
+               dstPtr1 <- writeArray sh (code . MultiValue.zip k) dstPtr0
+               (,) dstPtr1 <$> MultiValue.inc dstK)
+   finalN <-
+      Shape.sequenceShapeFromIndex . snd
+         =<< Shape.loop step n (ptr, MultiValue.zero)
+   return $ MultiValue.zip finalN sh
+
+
+scatterMaybe ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    Marshal.C sh1,
+    Storable.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array sh1 a -> Sym.Array sh0 (Maybe (ix1, a)) ->
+   LLVM.Value (MarshalPtr sh1) ->
+   LLVM.Value (Ptr a) ->
+   LLVM.CodeGenFunction r ()
+scatterMaybe accum (Sym.Array esh codeInit) (Sym.Array eish codeMap)
+      sptr ptr = do
+
+   sh <- Shape.load esh sptr
+   void $ writeArray sh codeInit ptr
+
+   ish <- unExp eish
+   let fill ix () = do
+         (MultiValue.Cons c, (jx, a)) <-
+            mapSnd MultiValue.unzip . MultiValue.splitMaybe <$> codeMap ix
+         C.ifThen c () $ do
+            p <- getElementPtr sh ptr jx
+            flip Storable.store p
+               =<< Expr.unliftM2 (flip accum) a
+               =<< Storable.load p
+   Shape.loop fill ish ()
+
+scatter ::
+   (Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    Marshal.C sh1,
+    Storable.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Sym.Array sh1 a ->
+   Sym.Array sh0 (Shape.Index sh1, a) ->
+   LLVM.Value (MarshalPtr sh1) ->
+   LLVM.Value (Ptr a) ->
+   LLVM.CodeGenFunction r ()
+scatter accum (Sym.Array esh codeInit) (Sym.Array eish codeMap) sptr ptr = do
+   sh <- Shape.load esh sptr
+   void $ writeArray sh codeInit ptr
+
+   ish <- unExp eish
+   let fill ix () = do
+         (jx, a) <- MultiValue.unzip <$> codeMap ix
+         p <- getElementPtr sh ptr jx
+         flip Storable.store p
+            =<< Expr.unliftM2 (flip accum) a
+            =<< Storable.load p
+   Shape.loop fill ish ()
+
+addDimension ::
+   (Shape.C n, Marshal.C n, Shape.Index n ~ k,
+    Shape.C sh, Marshal.C sh,
+    Storable.C b) =>
+   Exp n ->
+   (Exp k -> Exp a -> Exp b) ->
+   Sym.Array sh a ->
+   LLVM.Value (MarshalPtr (sh,n)) ->
+   LLVM.Value (Ptr b) ->
+   LLVM.CodeGenFunction r ()
+addDimension en select (Sym.Array esh code) sptr ptr = do
+   (sh,n) <- MultiValue.unzip <$> Shape.load (Expr.zip esh en) sptr
+
+   let fill ix ptr0 = do
+         a <- code ix
+         writeArray n (\k -> Expr.unliftM2 select k a) ptr0
+   void $ Shape.loop fill sh ptr
diff --git a/src/Data/Array/Knead/Symbolic/Private.hs b/src/Data/Array/Knead/Symbolic/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Private.hs
@@ -0,0 +1,204 @@
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic.Private where
+
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+
+import LLVM.DSL.Expression (Exp(Exp))
+
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import qualified LLVM.Extra.Iterator as Iter
+import qualified LLVM.Extra.Maybe as Maybe
+import qualified LLVM.Core as LLVM
+
+import qualified Control.Category as Cat
+import qualified Control.Monad.HT as Monad
+import Control.Monad ((<=<), )
+
+import Prelude hiding (id, map, zipWith, replicate, )
+
+
+type Val = MultiValue.T
+type Code r a = LLVM.CodeGenFunction r (Val a)
+
+data Array sh a =
+   Array (Exp sh) (forall r. Val (Shape.Index sh) -> Code r a)
+
+shape :: Array sh a -> Exp sh
+shape (Array sh _) = sh
+
+(!) ::
+   (Shape.C sh,  Shape.Index sh  ~ ix) =>
+   Array sh a -> Exp ix -> Exp a
+(!) (Array _ code) (Exp ix) = Exp (code =<< ix)
+
+the :: (Shape.Scalar sh) => Array sh a -> Exp a
+the (Array z code) = Exp (code $ Shape.zeroIndex z)
+
+fromScalar :: (Shape.Scalar sh) => Exp a -> Array sh a
+fromScalar = fill Shape.scalar
+
+
+fill :: Exp sh -> Exp a -> Array sh a
+fill sh (Exp code) = Array sh (\_z -> code)
+
+
+{- |
+This class allows to implement functions without parameters
+for both simple and parameterized arrays.
+-}
+class C array where
+   lift0 :: Array sh a -> array sh a
+   lift1 :: (Array sha a -> Array shb b) -> array sha a -> array shb b
+   lift2 ::
+      (Array sha a -> Array shb b -> Array shc c) ->
+      array sha a -> array shb b -> array shc c
+
+instance C Array where
+   lift0 = Cat.id
+   lift1 = Cat.id
+   lift2 = Cat.id
+
+
+gather ::
+   (C array,
+    Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    MultiValue.C a) =>
+   array sh1 ix0 ->
+   array sh0 a ->
+   array sh1 a
+gather =
+   lift2 $ \(Array sh1 f) (Array _sh0 code) ->
+      Array sh1 (code <=< f)
+
+backpermute2 ::
+   (C array,
+    Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    Shape.C sh,  Shape.Index sh  ~ ix) =>
+   Exp sh ->
+   (Exp ix -> Exp ix0) ->
+   (Exp ix -> Exp ix1) ->
+   (Exp a -> Exp b -> Exp c) ->
+   array sh0 a -> array sh1 b -> array sh c
+backpermute2 sh projectIndex0 projectIndex1 f =
+   lift2 $ \(Array _sha codeA) (Array _shb codeB) ->
+      Array sh
+         (\ix ->
+            Monad.liftJoin2 (Expr.unliftM2 f)
+               (codeA =<< Expr.unliftM1 projectIndex0 ix)
+               (codeB =<< Expr.unliftM1 projectIndex1 ix))
+
+
+id ::
+   (C array, Shape.C sh, Shape.Index sh ~ ix) =>
+   Exp sh -> array sh ix
+id sh = lift0 $ Array sh return
+
+map ::
+   (C array, Shape.C sh) =>
+   (Exp a -> Exp b) ->
+   array sh a -> array sh b
+map f =
+   lift1 $ \(Array sh code) ->
+      Array sh (Expr.unliftM1 f <=< code)
+
+mapWithIndex ::
+   (C array, Shape.C sh, Shape.Index sh ~ ix) =>
+   (Exp ix -> Exp a -> Exp b) ->
+   array sh a -> array sh b
+mapWithIndex f =
+   lift1 $ \(Array sh code) ->
+      Array sh (\ix -> Expr.unliftM2 f ix =<< code ix)
+
+
+fold1Code ::
+   (Shape.C sh, Shape.Index sh ~ ix, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Exp sh ->
+   (Val ix -> Code r a) ->
+   Code r a
+fold1Code f (Exp nc) code = do
+   n <- nc
+   fmap Maybe.fromJust $
+      Shape.loop
+         (\i0 macc0 -> do
+            a <- code i0
+            acc1 <- Maybe.run macc0 (return a) (flip (Expr.unliftM2 f) a)
+            return $ Maybe.just acc1)
+         n Maybe.nothing
+
+fold1 ::
+   (C array, Shape.C sh0, Shape.C sh1, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   array (sh0, sh1) a -> array sh0 a
+fold1 f =
+   lift1 $ \(Array shs code) ->
+      case Expr.unzip shs of
+         (sh, s) -> Array sh $ fold1Code f s . MultiValue.curry code
+
+
+fold1All ::
+   (Shape.C sh, MultiValue.C a) =>
+   (Exp a -> Exp a -> Exp a) ->
+   Array sh a -> Exp a
+fold1All f (Array sh code) = Exp (fold1Code f sh code)
+
+
+findAllCode ::
+   (Shape.C sh, Shape.Index sh ~ ix, MultiValue.C a) =>
+   (Exp a -> Exp Bool) ->
+   Exp sh ->
+   (Val ix -> Code r a) ->
+   Code r (Maybe a)
+findAllCode p (Exp sh) code = do
+   n <- sh
+   finalFound <-
+      Iter.mapWhileState_
+         (\a _found -> do
+            MultiValue.Cons b <- Expr.unliftM1 p a
+            notb <- LLVM.inv b
+            return (notb, Maybe.fromBool b a))
+         (Iter.mapM code $ Shape.iterator n)
+         Maybe.nothing
+   Maybe.run finalFound
+      (return MultiValue.nothing)
+      (return . MultiValue.just)
+
+{- |
+In principle this can be implemented using fold1All
+but it has a short-cut semantics.
+@All@ means that it scans all dimensions
+but it does not mean that it finds all occurrences.
+If you want to get the index of the found element,
+please decorate the array elements with their indices before calling 'findAll'.
+-}
+findAll ::
+   (Shape.C sh, MultiValue.C a) =>
+   (Exp a -> Exp Bool) ->
+   Array sh a -> Exp (Maybe a)
+findAll p (Array sh code) = Exp (findAllCode p sh code)
+
+
+class Process proc where
+
+
+infixl 3 $:.
+
+{- |
+Use this for combining several dimension manipulators.
+E.g.
+
+> apply (passAny $:. pick 3 $:. pass $:. replicate 10) array
+
+The constraint @(Process proc0, Process proc1)@ is a bit weak.
+We like to enforce that the type constructor like @Slice.T@
+is the same in @proc0@ and @proc1@, and only the parameters differ.
+Currently this coherence is achieved,
+because we only provide functions of type @proc0 -> proc1@ with this condition.
+-}
+($:.) :: (Process proc0, Process proc1) => proc0 -> (proc0 -> proc1) -> proc1
+($:.) = flip ($)
diff --git a/src/Data/Array/Knead/Symbolic/Render.hs b/src/Data/Array/Knead/Symbolic/Render.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Render.hs
@@ -0,0 +1,177 @@
+{-# LANGUAGE TypeFamilies #-}
+{- |
+Apply operations on symbolic arrays to physical ones.
+-}
+module Data.Array.Knead.Symbolic.Render (
+   run,
+   MarshalExp(..),
+   MapFilter(..),
+   FilterOuter(..),
+   Scatter(..),
+   ScatterMaybe(..),
+   MapAccumLSimple(..),
+   MapAccumLSequence(..),
+   MapAccumL(..),
+   FoldOuterL(..),
+   AddDimension(..),
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.Render.Basic as Render
+import qualified Data.Array.Knead.Symbolic.Render.Argument as Arg
+import qualified Data.Array.Knead.Symbolic.PhysicalParametric as PhysP
+import qualified Data.Array.Knead.Symbolic.Physical as Phys
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import qualified Data.Array.Knead.Shape as Shape
+import Data.Array.Knead.Symbolic.PhysicalParametric
+         (MapFilter, FilterOuter,
+          MapAccumLSimple, MapAccumLSequence, MapAccumL, FoldOuterL,
+          Scatter, ScatterMaybe, AddDimension)
+
+import qualified LLVM.DSL.Render.Run as Run
+import LLVM.DSL.Expression (Exp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+
+import Prelude2010
+import Prelude ()
+
+
+
+class C f where
+   type Plain f
+   function :: (Marshal.C p) => Run.T IO p f (Plain f)
+
+instance
+   (Marshal.C sh, Shape.C sh, Storable.C a) =>
+      C (Core.Array sh a) where
+   type Plain (Core.Array sh a) = IO (Phys.Array sh a)
+   function = Run.Cons PhysP.render
+
+instance
+   (Shape.Sequence n, Marshal.C n,
+    Storable.C b, MultiValue.C b) =>
+      C (MapFilter n a b) where
+   type Plain (MapFilter n a b) = IO (Phys.Array n b)
+   function = Run.Cons PhysP.mapFilter
+
+instance
+   (Shape.Sequence n, Marshal.C n,
+    Shape.C sh, Marshal.C sh,
+    Storable.C a, MultiValue.C a) =>
+      C (FilterOuter n sh a) where
+   type Plain (FilterOuter n sh a) = IO (Phys.Array (n,sh) a)
+   function = Run.Cons PhysP.filterOuter
+
+instance
+   (Shape.C sh0, Marshal.C sh0,
+    Shape.C sh1, Marshal.C sh1,
+    Storable.C a, MultiValue.C a) =>
+      C (Scatter sh0 sh1 a) where
+   type Plain (Scatter sh0 sh1 a) = IO (Phys.Array sh1 a)
+   function = Run.Cons PhysP.scatter
+
+instance
+   (Shape.C sh0, Marshal.C sh0,
+    Shape.C sh1, Marshal.C sh1,
+    Storable.C a, MultiValue.C a) =>
+      C (ScatterMaybe sh0 sh1 a) where
+   type Plain (ScatterMaybe sh0 sh1 a) = IO (Phys.Array sh1 a)
+   function = Run.Cons PhysP.scatterMaybe
+
+instance
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc,
+    Storable.C a, MultiValue.C a,
+    Storable.C b, MultiValue.C b) =>
+      C (MapAccumLSimple sh n acc a b) where
+   type Plain (MapAccumLSimple sh n acc a b) = IO (Phys.Array (sh,n) b)
+   function = Run.Cons PhysP.mapAccumLSimple
+
+instance
+   (Shape.C n, Marshal.C n,
+    MultiValue.C acc,
+    Storable.C final, MultiValue.C final,
+    Storable.C a, MultiValue.C a,
+    Storable.C b, MultiValue.C b) =>
+      C (MapAccumLSequence n acc final a b) where
+   type Plain (MapAccumLSequence n acc final a b) = IO (final, Phys.Array n b)
+   function = Run.Cons PhysP.mapAccumLSequence
+
+instance
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    MultiValue.C acc,
+    Storable.C final, MultiValue.C final,
+    Storable.C a, MultiValue.C a,
+    Storable.C b, MultiValue.C b) =>
+      C (MapAccumL sh n acc final a b) where
+   type Plain (MapAccumL sh n acc final a b) =
+            IO (Phys.Array sh final, Phys.Array (sh,n) b)
+   function = Run.Cons PhysP.mapAccumL
+
+instance
+   (Shape.C n, Marshal.C n,
+    Shape.C sh, Marshal.C sh,
+    Storable.C a, MultiValue.C a,
+    Storable.C b, MultiValue.C b) =>
+      C (FoldOuterL n sh a b) where
+   type Plain (FoldOuterL n sh a b) = IO (Phys.Array sh a)
+   function = Run.Cons PhysP.foldOuterL
+
+instance
+   (Shape.C sh, Marshal.C sh,
+    Shape.C n, Marshal.C n,
+    Storable.C b, MultiValue.C b) =>
+      C (AddDimension sh n a b) where
+   type Plain (AddDimension sh n a b) = IO (Phys.Array (sh,n) b)
+   function = Run.Cons PhysP.addDimension
+
+
+instance (Storable.C a, MultiValue.C a) => C (Exp a) where
+   type Plain (Exp a) = IO a
+   function = Render.storable
+
+newtype MarshalExp a = MarshalExp {getMarshalExp :: Exp a}
+
+instance (Marshal.C a) => C (MarshalExp a) where
+   type Plain (MarshalExp a) = IO a
+   function = Run.premapDSL getMarshalExp Render.marshal
+
+instance (Argument arg, C func) => C (arg -> func) where
+   type Plain (arg -> func) = PlainArg arg -> Plain func
+   function = argument Render.*-> function
+
+
+class Argument a where
+   type PlainArg a
+   argument :: Arg.T (PlainArg a) a
+
+instance Argument () where
+   type PlainArg () = ()
+   argument = Arg.unit
+
+instance
+   (Shape.C sh, Marshal.C sh, Storable.C a) =>
+      Argument (Core.Array sh a) where
+   type PlainArg (Core.Array sh a) = Phys.Array sh a
+   argument = Arg.array
+
+instance (Marshal.C a) => Argument (Exp a) where
+   type PlainArg (Exp a) = a
+   argument = Arg.primitive
+
+instance (Argument a, Argument b) => Argument (a,b) where
+   type PlainArg (a,b) = (PlainArg a, PlainArg b)
+   argument = Arg.pair argument argument
+
+instance (Argument a, Argument b, Argument c) => Argument (a,b,c) where
+   type PlainArg (a,b,c) = (PlainArg a, PlainArg b, PlainArg c)
+   argument = Arg.triple argument argument argument
+
+
+
+run :: (C f) => f -> IO (Plain f)
+run = Render.run function
diff --git a/src/Data/Array/Knead/Symbolic/Render/Argument.hs b/src/Data/Array/Knead/Symbolic/Render/Argument.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Render/Argument.hs
@@ -0,0 +1,47 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ExistentialQuantification #-}
+module Data.Array.Knead.Symbolic.Render.Argument (
+   Arg.T(Arg.Cons),
+   Arg.unit,
+   Arg.primitive,
+   Arg.pair,
+   Arg.triple,
+   array,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.Physical as Phys
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Code (getElementPtr)
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+
+import qualified LLVM.DSL.Render.Argument as Arg
+import LLVM.DSL.Expression (unExp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+
+import Foreign.ForeignPtr (withForeignPtr, touchForeignPtr)
+
+import Prelude2010
+import Prelude ()
+
+
+
+array ::
+   (Shape.C sh, Marshal.C sh, Storable.C a) =>
+   Arg.T (Phys.Array sh a) (Core.Array sh a)
+array =
+   Arg.Cons
+      (Expr.uncurry $ \esh eptr ->
+         Core.Array esh
+            (\ix -> do
+               sh <- unExp esh
+               MultiValue.Cons ptr <- unExp eptr
+               Storable.load =<< getElementPtr sh ptr ix))
+      (\(Array.Array sh fptr) ->
+         withForeignPtr fptr $ \ptr ->
+         return ((sh, ptr), touchForeignPtr fptr))
diff --git a/src/Data/Array/Knead/Symbolic/Render/Basic.hs b/src/Data/Array/Knead/Symbolic/Render/Basic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Render/Basic.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE Rank2Types #-}
+{- |
+Apply operations on symbolic arrays to physical ones.
+
+This is an approach with no pre-defined direction of type dependencies.
+-}
+module Data.Array.Knead.Symbolic.Render.Basic (
+   run,
+   (*->),
+
+   storable,
+   marshal,
+   array,
+   scatter,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.Render.Argument as Arg
+import qualified Data.Array.Knead.Symbolic.PhysicalParametric as PhysP
+import qualified Data.Array.Knead.Symbolic.Physical as Phys
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import qualified Data.Array.Knead.Shape as Shape
+
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
+
+import qualified LLVM.DSL.Render.Run as Run
+import LLVM.DSL.Render.Run (run, (*->))
+import LLVM.DSL.Expression (Exp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+
+import Prelude2010
+import Prelude ()
+
+
+
+_example1raw ::
+   (Marshal.C sh, Shape.C sh, Marshal.C z, Marshal.C a, Storable.C b) =>
+   Run.T IO z (Exp a -> Core.Array sh b) (a -> IO (Phys.Array sh b))
+_example1raw = Arg.primitive *-> array
+
+_example2raw ::
+   (Marshal.C sh, Shape.C sh,
+    Marshal.C z, Marshal.C a, Marshal.C b, Storable.C c) =>
+   Run.T IO z
+      (Exp a -> Exp b -> Core.Array sh c)
+      (a -> b -> IO (Phys.Array sh c))
+_example2raw = Arg.primitive *-> Arg.primitive *-> array
+
+
+_example2 ::
+   (Marshal.C sh, Shape.C sh,
+    Marshal.C a, Marshal.C b, Storable.C c) =>
+   (Exp a -> Exp b -> Core.Array sh c) ->
+   IO (a -> b -> IO (Phys.Array sh c))
+_example2 = run (Arg.primitive *-> Arg.primitive *-> array)
+
+_example2exp ::
+   (Marshal.C a, Marshal.C b, Storable.C c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   IO (a -> b -> IO c)
+_example2exp = run (Arg.primitive *-> Arg.primitive *-> storable)
+
+_example2marshal ::
+   (Marshal.C a, Marshal.C b, Marshal.C c) =>
+   (Exp a -> Exp b -> Exp c) ->
+   IO (a -> b -> IO c)
+_example2marshal = run (Arg.primitive *-> Arg.primitive *-> marshal)
+
+_example2scatter ::
+   (Shape.C sh0, Shape.C sh1, Marshal.C sh1,
+    Marshal.C a, Marshal.C b, Storable.C c) =>
+   (Exp a -> Exp b -> PhysP.Scatter sh0 sh1 c) ->
+   IO (a -> b -> IO (Array.Array sh1 c))
+_example2scatter = run (Arg.primitive *-> Arg.primitive *-> scatter)
+
+
+
+
+singleton :: Exp a -> Core.Array () a
+singleton = Core.fromScalar
+
+storable :: (Marshal.C p, Storable.C a) => Run.T IO p (Exp a) (IO a)
+storable = Run.Cons $ PhysP.the . fmap singleton
+
+marshal :: (Marshal.C p, Marshal.C a) => Run.T IO p (Exp a) (IO a)
+marshal = Run.Cons $ PhysP.theMarshal . fmap singleton
+
+array ::
+   (Shape.C sh, Shape.Index sh ~ ix, Marshal.C sh,
+    Marshal.C p, Storable.C a) =>
+   Run.T IO p (Core.Array sh a) (IO (Phys.Array sh a))
+array = Run.Cons PhysP.render
+
+
+scatter ::
+   (Shape.C sh0, Shape.C sh1, Marshal.C sh1, Marshal.C p, Storable.C a) =>
+   Run.T IO p (PhysP.Scatter sh0 sh1 a) (IO (Array.Array sh1 a))
+scatter = Run.Cons PhysP.scatter
diff --git a/src/Data/Array/Knead/Symbolic/RenderAlt.hs b/src/Data/Array/Knead/Symbolic/RenderAlt.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/RenderAlt.hs
@@ -0,0 +1,117 @@
+{-# LANGUAGE TypeFamilies #-}
+{- |
+Apply operations on symbolic arrays to physical ones.
+
+In contrast to the "Data.Array.Knead.Symbolic.Render" module,
+here we map from Haskell types to LLVM ones.
+This is analogous to "Synthesizer.LLVM.Generator.Render".
+-}
+module Data.Array.Knead.Symbolic.RenderAlt (
+   run,
+   MarshalValue(..),
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.Render.Basic as Render
+import qualified Data.Array.Knead.Symbolic.Render.Argument as Arg
+import qualified Data.Array.Knead.Symbolic.PhysicalParametric as PhysP
+import qualified Data.Array.Knead.Symbolic.Physical as Phys
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import qualified Data.Array.Knead.Shape as Shape
+
+import qualified LLVM.DSL.Render.Run as Run
+import LLVM.DSL.Expression (Exp)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+
+import Data.Word (Word, Word32)
+
+import Prelude2010
+import Prelude ()
+
+
+
+class C f where
+   type DSL f
+   function :: (Marshal.C p) => Run.T IO p (DSL f) f
+
+instance (C_IO a) => C (IO a) where
+   type DSL (IO a) = DSL_IO a
+   function = buildIO
+
+
+class C_IO f where
+   type DSL_IO f
+   buildIO :: (Marshal.C p) => Run.T IO p (DSL_IO f) (IO f)
+
+instance
+   (Marshal.C sh, Shape.C sh, Storable.C a) =>
+      C_IO (Phys.Array sh a) where
+   type DSL_IO (Phys.Array sh a) = Core.Array sh a
+   buildIO = Run.Cons PhysP.render
+
+
+instance C_IO Float where
+   type DSL_IO Float = Exp Float
+   buildIO = Render.storable
+
+instance C_IO Word32 where
+   type DSL_IO Word32 = Exp Word32
+   buildIO = Render.storable
+
+newtype MarshalValue a = MarshalValue {getMarshalValue :: a}
+
+instance (Marshal.C a) => C_IO (MarshalValue a) where
+   type DSL_IO (MarshalValue a) = Exp a
+   buildIO = Run.postmapPlain (fmap MarshalValue) Render.marshal
+
+
+instance (Argument arg, C func) => C (arg -> func) where
+   type DSL (arg -> func) = DSLArg arg -> DSL func
+   function = argument Render.*-> function
+
+
+
+class Argument a where
+   type DSLArg a
+   argument :: Arg.T a (DSLArg a)
+
+instance Argument () where
+   type DSLArg () = ()
+   argument = Arg.unit
+
+instance
+   (Shape.C sh, Marshal.C sh, Storable.C a) =>
+      Argument (Phys.Array sh a) where
+   type DSLArg (Phys.Array sh a) = Core.Array sh a
+   argument = Arg.array
+
+
+instance Argument Float where
+   type DSLArg Float = Exp Float
+   argument = Arg.primitive
+
+instance Argument Int where
+   type DSLArg Int = Exp Int
+   argument = Arg.primitive
+
+instance Argument Word where
+   type DSLArg Word = Exp Word
+   argument = Arg.primitive
+
+instance Argument Word32 where
+   type DSLArg Word32 = Exp Word32
+   argument = Arg.primitive
+
+instance (Argument a, Argument b) => Argument (a,b) where
+   type DSLArg (a,b) = (DSLArg a, DSLArg b)
+   argument = Arg.pair argument argument
+
+instance (Argument a, Argument b, Argument c) => Argument (a,b,c) where
+   type DSLArg (a,b,c) = (DSLArg a, DSLArg b, DSLArg c)
+   argument = Arg.triple argument argument argument
+
+
+
+run :: (C f) => DSL f -> IO f
+run = Render.run function
diff --git a/src/Data/Array/Knead/Symbolic/ShapeDependent.hs b/src/Data/Array/Knead/Symbolic/ShapeDependent.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/ShapeDependent.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic.ShapeDependent where
+
+import qualified Data.Array.Knead.Symbolic.Private as Core
+import Data.Array.Knead.Symbolic.Private (Array(Array), )
+
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified Control.Monad.HT as Monad
+import Control.Monad ((<=<), )
+
+
+shape :: (Core.C array, Shape.C sh, Shape.Scalar z) => array sh a -> array z sh
+shape = Core.lift1 $ Core.fromScalar . Core.shape
+
+backpermute ::
+   (Core.C array,
+    Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1) =>
+   (Exp sh0 -> Exp sh1) ->
+   (Exp ix1 -> Exp ix0) ->
+   array sh0 a ->
+   array sh1 a
+backpermute createShape projectIndex =
+   Core.lift1 $ \(Array sh code) ->
+      Array (createShape sh)
+         (code <=< Expr.unliftM1 projectIndex)
+
+{- |
+This is between 'backpermute' and 'backpermute2'.
+You can access the shapes of two arrays,
+but only the content of one of them.
+This is necessary if the second array contributes only a virtual dimension.
+-}
+backpermuteExtra ::
+   (Core.C array,
+    Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    Shape.C sh,  Shape.Index sh  ~ ix) =>
+   (Exp sh0 -> Exp sh1 -> Exp sh) ->
+   (Exp ix -> Exp ix0) ->
+   array sh0 a -> array sh1 b -> array sh a
+backpermuteExtra newShape projectIndex =
+   Core.lift2 $ \(Array sh0 code) (Array sh1 _code) ->
+      Array (newShape sh0 sh1)
+         (\ix -> code =<< Expr.unliftM1 projectIndex ix)
+
+backpermute2 ::
+   (Core.C array,
+    Shape.C sh0, Shape.Index sh0 ~ ix0,
+    Shape.C sh1, Shape.Index sh1 ~ ix1,
+    Shape.C sh,  Shape.Index sh  ~ ix) =>
+   (Exp sh0 -> Exp sh1 -> Exp sh) ->
+   (Exp ix -> Exp ix0) ->
+   (Exp ix -> Exp ix1) ->
+   (Exp a -> Exp b -> Exp c) ->
+   array sh0 a -> array sh1 b -> array sh c
+backpermute2 combineShape projectIndex0 projectIndex1 f =
+   Core.lift2 $ \(Array sha codeA) (Array shb codeB) ->
+      Array (combineShape sha shb)
+         (\ix ->
+            Monad.liftJoin2 (Expr.unliftM2 f)
+               (codeA =<< Expr.unliftM1 projectIndex0 ix)
+               (codeB =<< Expr.unliftM1 projectIndex1 ix))
+
+fill ::
+   (Core.C array) =>
+   (Exp sh0 -> Exp sh1) -> Exp b ->
+   array sh0 a -> array sh1 b
+fill fsh a =
+   Core.lift1 $ \arr ->
+      Core.fill (fsh $ Core.shape arr) a
diff --git a/src/Data/Array/Knead/Symbolic/Slice.hs b/src/Data/Array/Knead/Symbolic/Slice.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Array/Knead/Symbolic/Slice.hs
@@ -0,0 +1,198 @@
+{- |
+Generate and apply index maps.
+This unifies the @replicate@ and @slice@ functions of the @accelerate@ package.
+However the structure of slicing and replicating cannot depend on parameters.
+If you need that, you must use 'ShapeDep.backpermute' and friends.
+-}
+{-
+Some notes on the design choice:
+
+Instead of the shallow embedding implemented by the 'T' type,
+we could maintain a symbolic representation of the Slice and Replicate pattern,
+like the accelerate package does.
+We actually used that representation in former versions.
+It has however some drawbacks:
+
+* We need additional type functions that map from the pattern
+  to the source and the target shape and we need a proof,
+  that the images of these type functions are actually shapes.
+  This worked already, but was rather cumbersome.
+
+* We need a way to store and pass this pattern through the Parameter handler.
+  This yields new problems:
+  We need a wrapper type for wrapping Index, Shape, Slice, Replicate, Fold patterns.
+  Then the question is whether we use one Wrap type with a phantom parameter
+  or whether we define a Wrap type for every pattern type.
+  That is, the options are to write either
+
+  > Wrap Shape (Z:.Int:.Int)
+
+  or
+
+  > Shape (Z:.Int:.Int)
+
+  The first one seems to save us many duplicate instances of
+  Storable, MultiValue etc.
+  and it allows us easily to reuse the (:.) for all kinds of patterns.
+  However, we need a way to restrict the element type of the (:.)-list elements.
+  We can define that using variable ConstraintKinds,
+  but e.g. we are not able to add a Storable superclass constraint
+  to the instance Storable (Wrap constr).
+  That is, we are left with the second option
+  and had to define a lot of similar Storable, MultiValue instances.
+-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE ExistentialQuantification #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Data.Array.Knead.Symbolic.Slice (
+   T,
+   Cubic,
+   apply,
+   passAny,
+   pass,
+   pick,
+   pickFst,
+   pickSnd,
+   extrude,
+   extrudeFst,
+   extrudeSnd,
+   transpose,
+   (Core.$:.),
+
+   id,
+   first,
+   second,
+   compose,
+   ) where
+
+import qualified Data.Array.Knead.Symbolic.ShapeDependent as ShapeDep
+import qualified Data.Array.Knead.Symbolic.Private as Core
+
+import qualified Data.Array.Knead.Shape.Cubic.Int as Index
+import qualified Data.Array.Knead.Shape.Cubic as Cubic
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Knead.Expression as Expr
+import Data.Array.Knead.Shape.Cubic ((#:.), (:.)((:.)), )
+import Data.Array.Knead.Expression (Exp, )
+
+import qualified LLVM.Extra.Multi.Value as MultiValue
+import LLVM.Extra.Multi.Value (atom, )
+
+import qualified Type.Data.Num.Unary as Unary
+
+import qualified Prelude as P
+import Prelude hiding (id, zipWith, zipWith3, zip, zip3, replicate, )
+
+
+
+{-
+This data type is almost identical to Core.Array.
+The only difference is,
+that the shape @sh1@ in T can depend on another shape @sh0@.
+-}
+data T sh0 sh1 =
+   forall ix0 ix1.
+   (Shape.Index sh0 ~ ix0, Shape.Index sh1 ~ ix1) =>
+   Cons
+      (Exp sh0 -> Exp sh1)
+      (Exp ix1 -> Exp ix0)
+
+{- |
+This is essentially a 'ShapeDep.backpermute'.
+-}
+apply ::
+   (Core.C array, Shape.C sh0, Shape.C sh1, MultiValue.C a) =>
+   T sh0 sh1 ->
+   array sh0 a ->
+   array sh1 a
+apply (Cons fsh fix) =
+   ShapeDep.backpermute fsh fix
+
+
+pickFst :: Exp (Shape.Index n) -> T (n,sh) sh
+pickFst i = Cons Expr.snd (Expr.zip i)
+
+pickSnd :: Exp (Shape.Index n) -> T (sh,n) sh
+pickSnd i = Cons Expr.fst (flip Expr.zip i)
+
+{- |
+Extrusion has the potential to do duplicate work.
+Only use it to add dimensions of size 1, e.g. numeric 1 or unit @()@
+or to duplicate slices of physical arrays.
+-}
+extrudeFst :: Exp n -> T sh (n,sh)
+extrudeFst n = Cons (Expr.zip n) Expr.snd
+
+extrudeSnd :: Exp n -> T sh (sh,n)
+extrudeSnd n = Cons (flip Expr.zip n) Expr.fst
+
+transpose :: T (sh0,sh1) (sh1,sh0)
+transpose = Cons Expr.swap Expr.swap
+
+
+-- Arrow combinators
+
+id :: T sh sh
+id = Cons P.id P.id
+
+first :: T sh0 sh1 -> T (sh0,sh) (sh1,sh)
+first (Cons fsh fix) = Cons (Expr.mapFst fsh) (Expr.mapFst fix)
+
+second :: T sh0 sh1 -> T (sh,sh0) (sh,sh1)
+second (Cons fsh fix) = Cons (Expr.mapSnd fsh) (Expr.mapSnd fix)
+
+infixr 1 `compose`
+
+compose :: T sh0 sh1 -> T sh1 sh2 -> T sh0 sh2
+compose (Cons fshA fixA) (Cons fshB fixB) = Cons (fshB . fshA) (fixA . fixB)
+
+
+type Cubic rank0 rank1 = T (Cubic.Shape rank0) (Cubic.Shape rank1)
+
+{- |
+Like @Any@ in @accelerate@.
+-}
+passAny :: Cubic rank rank
+passAny = Cons P.id P.id
+
+{- |
+Like @All@ in @accelerate@.
+-}
+pass ::
+   (Unary.Natural rank0, Unary.Natural rank1) =>
+   Cubic rank0 rank1 ->
+   Cubic (Unary.Succ rank0) (Unary.Succ rank1)
+pass (Cons fsh fix) =
+   Cons
+      (Expr.modify (atom:.atom) $ \(sh:.s) -> fsh sh :. s)
+      (Expr.modify (atom:.atom) $ \(ix:.i) -> fix ix :. i)
+
+{- |
+Like @Int@ in @accelerate/slice@.
+-}
+pick ::
+   (Unary.Natural rank0, Unary.Natural rank1) =>
+   Exp Index.Int ->
+   Cubic rank0 rank1 ->
+   Cubic (Unary.Succ rank0) rank1
+pick i (Cons fsh fix) =
+   Cons
+      (fsh . Cubic.tail)
+      (\ix -> fix ix #:. i)
+
+{- |
+Like @Int@ in @accelerate/replicate@.
+-}
+extrude ::
+   (Unary.Natural rank0, Unary.Natural rank1) =>
+   Exp Index.Int ->
+   Cubic rank0 rank1 ->
+   Cubic rank0 (Unary.Succ rank1)
+extrude n (Cons fsh fix) =
+   Cons
+      (\sh -> fsh sh #:. n)
+      (fix . Cubic.tail)
+
+
+instance Core.Process (T sh0 sh1) where
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,18 @@
+module Main where
+
+import qualified Test.Array as Array
+
+import qualified LLVM.Core as LLVM
+
+import Data.Tuple.HT (mapFst)
+
+import qualified Test.QuickCheck as QC
+
+
+main :: IO ()
+main = do
+   LLVM.initializeNativeTarget
+
+   mapM_ (\(msg,prop) -> putStr (msg++": ") >> prop >>= QC.quickCheck) $
+      map (mapFst ("Array."++)) Array.tests ++
+      []
diff --git a/test/Test/Array.hs b/test/Test/Array.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Array.hs
@@ -0,0 +1,101 @@
+module Test.Array where
+
+import qualified Data.Array.Knead.Symbolic.Render as Render
+import qualified Data.Array.Knead.Symbolic as Symb
+import qualified Data.Array.Knead.Symbolic.Slice as Slice
+import qualified Data.Array.Knead.Expression as Expr
+import qualified Data.Array.Knead.Shape as Shape
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as ComfortShape
+import Data.Array.Comfort.Storable (Array)
+
+import qualified LLVM.Extra.Multi.Value.Storable as Storable
+import qualified LLVM.Extra.Multi.Value.Marshal as Marshal
+import qualified LLVM.Extra.Multi.Value as MultiValue
+
+import qualified LLVM.Core as LLVM
+
+import qualified Type.Data.Num.Decimal as TypeNum
+
+import Foreign.Storable (Storable)
+
+import qualified Data.List.HT as ListHT
+import Data.Int (Int32, Int64)
+
+import Control.Applicative ((<$>))
+
+import qualified Test.QuickCheck.Monadic as QCMon
+import qualified Test.QuickCheck as QC
+
+
+type Dim = ComfortShape.ZeroBased Int64
+type Dim2 = (Dim, Dim)
+
+genArray :: (QC.Arbitrary a, Storable a) => QC.Gen (Array Dim2 a)
+genArray = do
+   m <- QC.choose (1,10)
+   n <- QC.choose (1,10)
+   let shape = (Shape.ZeroBased m, Shape.ZeroBased n)
+   Array.fromList shape <$> QC.vector (ComfortShape.size shape)
+
+
+rowSumSymb ::
+   (Shape.C sh0, Shape.C sh1, MultiValue.Additive a) =>
+   Symb.Array (sh0,sh1) a -> Symb.Array sh0 a
+rowSumSymb = Symb.fold1 Expr.add
+
+columnSumSymb ::
+   (Shape.C sh0, Shape.C sh1, MultiValue.Additive a) =>
+   Symb.Array (sh0,sh1) a -> Symb.Array sh1 a
+columnSumSymb = Symb.fold1 Expr.add . Slice.apply Slice.transpose
+
+
+getRows ::
+   (ComfortShape.C sh0, ComfortShape.C sh1, Storable a) =>
+   Array (sh0,sh1) a -> [[a]]
+getRows x =
+   ListHT.sliceVertical
+      (ComfortShape.size $ snd $ Array.shape x)
+      (Array.toList x)
+
+rowPred ::
+   (Num a, Eq a, Storable a,
+    ComfortShape.C sh0, ComfortShape.C sh1) =>
+   Array (sh0, sh1) a -> Array sh0 a -> Bool
+rowPred x y  =  Array.toList y == map sum (getRows x)
+
+columnPred ::
+   (Num a, Eq a, Storable a,
+    ComfortShape.C sh0, ComfortShape.C sh1) =>
+   Array (sh0, sh1) a -> Array sh1 a -> Bool
+columnPred x y  =  Array.toList y == foldl1 (zipWith (+)) (getRows x)
+
+run ::
+   (Shape.C sh0, Marshal.C sh0, Show sh0,
+    Shape.C sh1, Marshal.C sh1, Show sh1,
+    Show a, Num a, Eq a, Storable.C a) =>
+   QC.Gen (Array sh0 a) ->
+   (Symb.Array sh0 a -> Symb.Array sh1 a) ->
+   (Array sh0 a -> Array sh1 a -> Bool) ->
+   IO QC.Property
+run qcgen code predicate = do
+   act <- Render.run code
+   return $ QC.forAll qcgen $ \x ->
+      QCMon.monadicIO $ do
+         y <- QCMon.run $ act x
+         QCMon.assert $ predicate x y
+
+
+tests :: [(String, IO QC.Property)]
+tests =
+   ("rowSum",
+      run (genArray :: QC.Gen (Array Dim2 Int32)) rowSumSymb rowPred) :
+   ("columnSum",
+      run (genArray :: QC.Gen (Array Dim2 Int32)) columnSumSymb columnPred) :
+   ("rowSumV3",
+      run (genArray :: QC.Gen (Array Dim2 (LLVM.Vector TypeNum.D3 Int32)))
+         rowSumSymb rowPred) :
+   ("columnSumV3",
+      run (genArray :: QC.Gen (Array Dim2 (LLVM.Vector TypeNum.D3 Int32)))
+         columnSumSymb columnPred) :
+   []