diff --git a/Data/Yarr/Base.hs b/Data/Yarr/Base.hs
--- a/Data/Yarr/Base.hs
+++ b/Data/Yarr/Base.hs
@@ -210,26 +210,29 @@
 -- | Internal implementation class. Generalizes @linear-@ and simple
 -- indexing and writing function in 'USource' and 'UTarget' classes.
 class (Shape sh, Shape i) => WorkIndex sh i where
+    toWork :: sh -> i
     gindex :: USource r l sh a => UArray r l sh a -> i -> IO a
     gwrite :: UTarget tr tl sh a => UArray tr tl sh a -> i -> a -> IO ()
     gsize :: USource r l sh a => UArray r l sh a -> i
+    gsize = toWork. extent
+    {-# INLINE gsize #-}
 
 instance Shape sh => WorkIndex sh sh where
+    toWork = id
     gindex = index
     gwrite = write
-    gsize = extent
     {-# INLINE gindex #-}
     {-# INLINE gwrite #-}
-    {-# INLINE gsize #-}
+    
 
 #define WI_INT_INST(sh)           \
 instance WorkIndex sh Int where { \
+    toWork = size;                \
     gindex = linearIndex;         \
     gwrite = linearWrite;         \
-    gsize = size . extent;        \
+    {-# INLINE toWork #-};        \
     {-# INLINE gindex #-};        \
     {-# INLINE gwrite #-};        \
-    {-# INLINE gsize #-};         \
 }
 
 WI_INT_INST(Dim2)
diff --git a/Data/Yarr/Convolution/Eval.hs b/Data/Yarr/Convolution/Eval.hs
--- a/Data/Yarr/Convolution/Eval.hs
+++ b/Data/Yarr/Convolution/Eval.hs
@@ -7,7 +7,7 @@
 import Data.Yarr.Convolution.Repr
 import Data.Yarr.Repr.Separate
 
-import Data.Yarr.Utils.FixedVector as V
+import Data.Yarr.Utils.FixedVector as V hiding (zero)
 import Data.Yarr.Utils.Fork
 import Data.Yarr.Utils.Parallel
 import Data.Yarr.Utils.Split
diff --git a/Data/Yarr/Eval.hs b/Data/Yarr/Eval.hs
--- a/Data/Yarr/Eval.hs
+++ b/Data/Yarr/Eval.hs
@@ -25,7 +25,7 @@
 import Data.Yarr.Base as B
 import Data.Yarr.Shape as S
 
-import Data.Yarr.Utils.FixedVector as V
+import Data.Yarr.Utils.FixedVector as V hiding (zero)
 import Data.Yarr.Utils.Fork
 import Data.Yarr.Utils.Parallel
 import Data.Yarr.Utils.Primitive as P
diff --git a/Data/Yarr/Flow.hs b/Data/Yarr/Flow.hs
--- a/Data/Yarr/Flow.hs
+++ b/Data/Yarr/Flow.hs
@@ -13,10 +13,15 @@
     dzipElems, dzipElemsM,
 
     -- * High level shortcuts
-    traverse, zipElems, mapElems, mapElemsM
+    traverse, zipElems, mapElems, mapElemsM,
+
+    -- ** Cartesian products
+    cartProduct2, icartProduct2, icartProduct2M,
+    cartProduct3, icartProduct3, icartProduct3M,
 ) where
 
 import Data.Yarr.Base
+import Data.Yarr.Shape
 import Data.Yarr.Fusion
 import Data.Yarr.Repr.Delayed
 import Data.Yarr.Repr.Separate
@@ -127,3 +132,84 @@
     -> UArray D l sh d    -- ^ Result array
 {-# INLINE dzipWith3 #-}
 dzipWith3 f arr1 arr2 arr3 = dzip3 f (delay arr1) (delay arr2) (delay arr3)
+
+
+cartProduct2
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b)
+    => (a -> b -> c)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray D SH Dim2 c
+{-# INLINE cartProduct2 #-}
+cartProduct2 f = icartProduct2M (\_ x y -> return (f x y))
+
+icartProduct2
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b)
+    => (Dim2 -> a -> b -> c)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray D SH Dim2 c
+{-# INLINE icartProduct2 #-}
+icartProduct2 f = icartProduct2M (\ix x y -> return (f ix x y))
+
+icartProduct2M
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b)
+    => (Dim2 -> a -> b -> IO c)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray D SH Dim2 c
+{-# INLINE icartProduct2M #-}
+icartProduct2M f arr1 arr2 =
+    ShapeDelayed
+        (extent arr1, extent arr2)
+        (touchArray arr1 >> touchArray arr2)
+        (force arr1 >> force arr2)
+        indexF
+  where
+    {-# INLINE indexF #-}
+    indexF ix@(i, j) = do
+        a <- arr1 `index` i
+        b <- arr2 `index` j
+        f ix a b
+
+cartProduct3
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b, USource r3 l3 Dim1 c)
+    => (a -> b -> c -> d)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray r3 l3 Dim1 c
+    -> UArray D SH Dim3 d
+{-# INLINE cartProduct3 #-}
+cartProduct3 f = icartProduct3M (\_ x y z -> return (f x y z))
+
+icartProduct3
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b, USource r3 l3 Dim1 c)
+    => (Dim3 -> a -> b -> c -> d)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray r3 l3 Dim1 c
+    -> UArray D SH Dim3 d
+{-# INLINE icartProduct3 #-}
+icartProduct3 f = icartProduct3M (\ix x y z -> return (f ix x y z))
+
+icartProduct3M
+    :: (USource r1 l1 Dim1 a, USource r2 l2 Dim1 b, USource r3 l3 Dim1 c)
+    => (Dim3 -> a -> b -> c -> IO d)
+    -> UArray r1 l1 Dim1 a
+    -> UArray r2 l2 Dim1 b
+    -> UArray r3 l3 Dim1 c
+    -> UArray D SH Dim3 d
+{-# INLINE icartProduct3M #-}
+icartProduct3M f arr1 arr2 arr3 =
+    ShapeDelayed
+        (extent arr1, extent arr2, extent arr3)
+        (touchArray arr1 >> touchArray arr2 >> touchArray arr3)
+        (force arr1 >> force arr2 >> force arr3)
+        indexF
+  where
+    {-# INLINE indexF #-}
+    indexF ix@(i, j, k) = do
+        a <- arr1 `index` i
+        b <- arr2 `index` j
+        c <- arr3 `index` k
+        f ix a b c
diff --git a/Data/Yarr/Repr/Delayed.hs b/Data/Yarr/Repr/Delayed.hs
--- a/Data/Yarr/Repr/Delayed.hs
+++ b/Data/Yarr/Repr/Delayed.hs
@@ -13,7 +13,11 @@
     UArray(..),
 
     -- * Misc
-    L, SH, fromFunction, delay, delayShaped,
+    L, SH, fromFunction, fromLinearFunction,
+    delay, delayShaped, delayLinear,
+
+    -- ** Const arrays
+    linearConst, shapedConst,
 ) where
 
 import Prelude as P
@@ -234,7 +238,7 @@
 -- that has specialized implementation in the library (mapping, zipping, etc).
 --
 -- Suitable to obtain arrays of constant element,
--- of indices (@fromFunction sh 'id'@), and so on.
+-- of indices (@fromFunction sh 'return'@), and so on.
 fromFunction
     :: Shape sh
     => sh               -- ^ Extent of array
@@ -243,6 +247,24 @@
 {-# INLINE fromFunction #-}
 fromFunction sh f = ShapeDelayed sh (return ()) (return ()) f
 
+fromLinearFunction
+    :: Shape sh
+    => sh               -- ^ Extent of array
+    -> (Int -> IO a)     -- ^ Linear ndexing function
+    -> UArray D L sh a -- ^ Result array
+{-# INLINE fromLinearFunction #-}
+fromLinearFunction sh f = LinearDelayed sh (return ()) (return ()) f
+
+linearConst
+    :: Shape sh => sh -> a -> UArray D L sh a
+{-# INLINE linearConst #-}
+linearConst sh x = fromLinearFunction sh (const (return x))
+
+shapedConst
+    :: Shape sh => sh -> a -> UArray D SH sh a
+{-# INLINE shapedConst #-}
+shapedConst sh x = fromFunction sh (const (return x))
+
 -- | Wraps @('index' arr)@ into Delayed representation. Normally you shouldn't need
 -- to use this function. It may be dangerous for performance, because
 -- preferred 'Data.Yarr.Eval.Load'ing type of source array is ignored.
@@ -250,6 +272,11 @@
 {-# INLINE delayShaped #-}
 delayShaped arr =
     ShapeDelayed (extent arr) (touchArray arr) (force arr) (index arr)
+
+delayLinear :: USource r l sh a => UArray r l sh a -> UArray D L sh a
+{-# INLINE delayLinear #-}
+delayLinear arr =
+    LinearDelayed (extent arr) (touchArray arr) (force arr) (linearIndex arr)
 
 -- | In opposite to 'D'elayed (source) Delayed Target holds abstract /writing/
 -- function: @(sh -> a -> IO ())@. It may be used to perform arbitrarily tricky
diff --git a/Data/Yarr/Shape.hs b/Data/Yarr/Shape.hs
--- a/Data/Yarr/Shape.hs
+++ b/Data/Yarr/Shape.hs
@@ -5,7 +5,7 @@
     module Data.Yarr.WorkTypes,
 
     -- * Shape and BlockShape
-    Block, Shape(..), BlockShape(..),
+    Block, Shape(..), BlockShape(..), MultiShape(..),
 
     -- * Shape instances
     Dim1, Dim2, Dim3,
@@ -40,6 +40,7 @@
     zero :: sh
     -- | 'Dim1' @size@ is 'id', @size (3, 5) == 15@
     size :: sh -> Int
+    inc :: sh -> sh
     -- | @(1, 2, 3) \`plus\` (0, 0, 1) == (1, 2, 4)@
     plus :: sh -> sh -> sh
     -- | @(1, 2) \`minus\` (1, 0) == (0, 2)@ 
@@ -141,12 +142,18 @@
         -> VecList (BorderCount sh) (Block sh) -- ^ Shavings
 
 
+class (Shape sh, Shape lsh) => MultiShape sh lsh | sh -> lsh, lsh -> sh where
+    lower :: sh -> lsh
+    inner :: sh -> Int
+    combine :: lsh -> Int -> sh
 
+
 type Dim1 = Int
 
 instance Shape Dim1 where
     zero = 0
     size = id
+    inc = succ
     plus = (+)
     offset off i = i - off
     fromLinear _ i = i
@@ -183,6 +190,7 @@
 
     {-# INLINE zero #-}
     {-# INLINE size #-}
+    {-# INLINE inc #-}
     {-# INLINE plus #-}
     {-# INLINE offset #-}
     {-# INLINE fromLinear #-}
@@ -218,6 +226,7 @@
 instance Shape Dim2 where
     zero = (0, 0)
     size (h, w) = h * w
+    inc (h, w) = (h, w + 1)
     plus (y1, x1) (y2, x2) = (y1 + y2, x1 + x2)
     offset (offY, offX) (y, x) = (y - offY, x - offX)
     fromLinear (_, w) i = i `quotRem` w
@@ -347,6 +356,7 @@
 
     {-# INLINE zero #-}
     {-# INLINE size #-}
+    {-# INLINE inc #-}
     {-# INLINE plus #-}
     {-# INLINE offset #-}
     {-# INLINE fromLinear #-}
@@ -378,6 +388,14 @@
 
     {-# INLINE clipBlock #-}
 
+instance MultiShape Dim2 Dim1 where
+    lower = fst
+    inner = snd
+    combine = (,)
+    {-# INLINE lower #-}
+    {-# INLINE inner #-}
+    {-# INLINE combine #-}
+
 -- | 2D-unrolled filling to maximize profit from
 -- \"Global value numbering\" LLVM optimization.
 --
@@ -436,6 +454,7 @@
 instance Shape Dim3 where
     zero = (0, 0, 0)
     size (d, h, w) = d * h * w
+    inc (d, h, w) = (d, h, w + 1)
     plus (z1, y1, x1) (z2, y2, x2) = (z1 + z2, y1 + y2, x1 + x2)
     offset (offZ, offY, offX) (z, y, x) = (z - offZ, y - offY, x - offX)
     fromLinear (_, h, w) i =
@@ -562,6 +581,7 @@
 
     {-# INLINE zero #-}
     {-# INLINE size #-}
+    {-# INLINE inc #-}
     {-# INLINE plus #-}
     {-# INLINE offset #-}
     {-# INLINE fromLinear #-}
@@ -579,3 +599,11 @@
     {-# INLINE unrolledFoldl #-}
     {-# INLINE foldr #-}
     {-# INLINE unrolledFoldr #-}
+
+instance MultiShape Dim3 Dim2 where
+    lower (z, y, _) = (z, y)
+    inner (_, _, x) = x
+    combine (z, y) x = (z, y, x)
+    {-# INLINE lower #-}
+    {-# INLINE inner #-}
+    {-# INLINE combine #-}
diff --git a/Data/Yarr/Utils/FixedVector.hs b/Data/Yarr/Utils/FixedVector.hs
--- a/Data/Yarr/Utils/FixedVector.hs
+++ b/Data/Yarr/Utils/FixedVector.hs
@@ -5,7 +5,7 @@
     Fn, arity,
     
     -- * Missed utility
-    zipWith3, zipWithM_, apply, all, any,
+    zipWith3, zipWithM_, apply, all, any, zero,
     iifoldl, iifoldM,
 
     -- * Aliases and shortcuts
@@ -27,7 +27,9 @@
 
 ) where
 
-import Prelude hiding (foldl, zipWith, zipWith3, all, any, sequence_)
+import Prelude hiding (
+    foldl, zipWith, zipWith3,
+    all, any, sequence_, replicate)
 
 import Control.DeepSeq
 
@@ -92,6 +94,10 @@
 any :: Vector v a => (a -> Bool) -> v a -> Bool
 {-# INLINE any #-}
 any p = foldl (\a x -> a || (p x)) False
+
+zero :: (Vector v a, Num a) => v a
+{-# INLINE zero #-}
+zero = replicate 0
 
 
 iifoldl
diff --git a/Data/Yarr/Walk.hs b/Data/Yarr/Walk.hs
--- a/Data/Yarr/Walk.hs
+++ b/Data/Yarr/Walk.hs
@@ -1,12 +1,12 @@
 
 module Data.Yarr.Walk (
     -- * Fold combinators
-    -- | See source of these 4 functions
-    -- to construct more similar ones,
-    -- if you need.
-    reduceL, reduceLeftM,
-    reduceR, reduceRightM,
+    reduceL, ireduceL, reduceLeftM,
+    reduceR, ireduceR, reduceRightM,
 
+    -- * Inner dim reducers
+    reduceInner, ireduceInner,
+
     -- * Combinators to walk with mutable state
     -- | Added specially to improve performance
     -- of tasks like histogram filling.
@@ -29,11 +29,12 @@
 import Data.Yarr.Base
 import Data.Yarr.Shape as S
 import Data.Yarr.Eval
+import Data.Yarr.Repr.Delayed
 
 import Data.Yarr.Walk.Internal
 
 
--- | /O(0)/
+-- | /O(1)/
 reduceLeftM
     :: Foldl i a b        -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
     -> (b -> a -> IO b)   -- ^ Monadic left reduce
@@ -42,7 +43,7 @@
 {-# INLINE reduceLeftM #-}
 reduceLeftM foldl rf = foldl (\b _ a -> rf b a)
 
--- | /O(0)/
+-- | /O(1)/
 reduceL
     :: Foldl i a b        -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
     -> (b -> a -> b)      -- ^ Pure left reduce
@@ -51,7 +52,16 @@
 {-# INLINE reduceL #-}
 reduceL foldl rf = foldl (\b _ a -> return $ rf b a)
 
--- | /O(0)/
+-- | /O(1)/
+ireduceL
+    :: Foldl i a b        -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
+    -> (b -> i -> a -> b) -- ^ Pure indexed left reduce
+    -> StatefulWalk i a b -- ^ Result stateful walk to be passed
+                          -- to walk runners
+{-# INLINE ireduceL #-}
+ireduceL foldl rf = foldl (\b i a -> return $ rf b i a)
+
+-- | /O(1)/
 reduceRightM
     :: Foldr i a b         -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
     -> (a -> b -> IO b)    -- ^ Monadic right reduce
@@ -60,7 +70,7 @@
 {-# INLINE reduceRightM #-}
 reduceRightM foldr rf = foldr (\_ a b -> rf a b)
 
--- | /O(0)/
+-- | /O(1)/
 reduceR
     :: Foldr i a b        -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
     -> (a -> b -> b)      -- ^ Pure right reduce
@@ -69,8 +79,17 @@
 {-# INLINE reduceR #-}
 reduceR foldr rf = foldr (\_ a b -> return $ rf a b)
 
+-- | /O(1)/
+ireduceR
+    :: Foldr i a b        -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
+    -> (i -> a -> b -> b) -- ^ Pure indexed right reduce
+    -> StatefulWalk i a b -- ^ Result stateful walk to be passed
+                          -- to walk runners
+{-# INLINE ireduceR #-}
+ireduceR foldr rf = foldr (\i a b -> return $ rf i a b)
 
--- | /O(0)/
+
+-- | /O(1)/
 mutate
     :: Fill i a           -- ^ 'S.fill' or curried 'S.unrolledFill'.
                           -- If mutating is associative,
@@ -82,7 +101,7 @@
 {-# INLINE mutate #-}
 mutate fill mf = imutate fill (\s i -> mf s)
 
--- | /O(0)/ Version of 'mutate', accepts mutating function
+-- | /O(1)/ Version of 'mutate', accepts mutating function
 -- which additionaly accepts array index.
 imutate
     :: Fill i a               -- ^ 'S.fill' or curried 'S.unrolledFill'.
@@ -97,6 +116,25 @@
     fill index (imf s) start end
     return s
 
+-- | /O(1)/
+reduceInner
+    :: (USource r l sh a, MultiShape sh lsh, PreferredWorkIndex l sh i)
+    => StatefulWalk i a b
+    -> (lsh -> IO b)
+    -> UArray r l sh a
+    -> UArray D SH lsh b
+{-# INLINE reduceInner #-}
+reduceInner = anyReduceInner
+
+-- | /O(1)/
+ireduceInner
+    :: (USource r l sh a, MultiShape sh lsh)
+    => StatefulWalk sh a b
+    -> (lsh -> IO b)
+    -> UArray r l sh a
+    -> UArray D SH lsh b
+{-# INLINE ireduceInner #-}
+ireduceInner = anyReduceInner
 
 
 -- | /O(n)/ Walk with state,
diff --git a/Data/Yarr/Walk/Internal.hs b/Data/Yarr/Walk/Internal.hs
--- a/Data/Yarr/Walk/Internal.hs
+++ b/Data/Yarr/Walk/Internal.hs
@@ -9,11 +9,27 @@
 import Data.Yarr.Base
 import Data.Yarr.Shape as S
 import Data.Yarr.Eval
+import Data.Yarr.Repr.Delayed
 
-import Data.Yarr.Utils.FixedVector as V hiding (toList)
+import Data.Yarr.Utils.FixedVector as V hiding (toList, zero)
 import Data.Yarr.Utils.Fork
 import Data.Yarr.Utils.Parallel
 
+anyReduceInner
+    :: (USource r l sh a, MultiShape sh lsh, WorkIndex sh i)
+    => StatefulWalk i a b
+    -> (lsh -> IO b)
+    -> UArray r l sh a
+    -> UArray D SH lsh b
+{-# INLINE anyReduceInner #-}
+anyReduceInner fold getZ arr =
+    ShapeDelayed (lower sh) (touchArray arr) (force arr) ix
+  where
+    sh = extent arr
+    ix lsh =
+        fold (getZ lsh) (gindex arr)
+             (toWork (combine lsh 0))
+             (toWork (combine (inc lsh) (inner sh)))
 
 anyWalk
     :: (USource r l sh a, WorkIndex sh i)
diff --git a/Debug/Yarr.hs b/Debug/Yarr.hs
--- a/Debug/Yarr.hs
+++ b/Debug/Yarr.hs
@@ -18,7 +18,7 @@
 
 import Data.Yarr.Base hiding (fmap)
 import Data.Yarr.Shape
-import Data.Yarr.Utils.FixedVector as V
+import Data.Yarr.Utils.FixedVector as V hiding (zero)
 
 -- | Yarr something to stderr.
 yarr :: String -> IO ()
diff --git a/yarr.cabal b/yarr.cabal
--- a/yarr.cabal
+++ b/yarr.cabal
@@ -1,5 +1,5 @@
 Name:                yarr
-Version:             1.3.1
+Version:             1.3.2
 Synopsis:            Yet another array library
 Description:
     Yarr is a new blazing fast dataflow framework (array library),
@@ -56,7 +56,7 @@
         ghc-prim         == 0.3.*,
         deepseq          == 1.3.*,
         fixed-vector     == 0.1.2.1,
-        primitive        >= 0.2,
+        primitive        >= 0.2 && < 0.6,
         template-haskell == 2.8.*,
         missing-foreign  == 0.1.1