diff --git a/Data/Yarr.hs b/Data/Yarr.hs
--- a/Data/Yarr.hs
+++ b/Data/Yarr.hs
@@ -71,6 +71,9 @@
     Currently there is only one option \"out of the box\" - to load image :)
     See "Data.Yarr.IO.Image" module in @yarr-image-io@ package.
 
+    Consider also "Data.Yarr.IO.List" module,
+    although it is very slow way to obtain manifest array in memory.
+
     /How to map and zip arrays:/
 
     See 'DefaultFusion' class and functions in "Data.Yarr.Flow" module.
@@ -181,7 +184,7 @@
     F, unsafeFromForeignPtr, toForeignPtr,
     
     -- ** Delayed
-    D, UArray(LinearDelayed, ShapeDelayed), delay,
+    D, UArray(LinearDelayed, ShapeDelayed), fromFunction, delay,
 
     -- ** Separate
     SE, fromSlices, unsafeMapSlices
@@ -190,7 +193,7 @@
 
 import Data.Yarr.Base hiding (Fusion(..))
 import Data.Yarr.Eval
-import Data.Yarr.Flow
+import Data.Yarr.Flow hiding (D, delay, SE)
 import Data.Yarr.Shape
 import Data.Yarr.Repr.Foreign
 import Data.Yarr.Repr.Delayed
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
@@ -13,6 +13,8 @@
 import Data.Yarr.Utils.Split
 
 
+instance Shape sh => PreferredWorkIndex CVL sh sh
+
 instance (BlockShape sh, UTarget tr tl sh a) =>
         Load CV CVL tr tl sh a where
     type LoadIndex CVL tl sh = sh
diff --git a/Data/Yarr/Eval.hs b/Data/Yarr/Eval.hs
--- a/Data/Yarr/Eval.hs
+++ b/Data/Yarr/Eval.hs
@@ -8,6 +8,8 @@
     -- * Load classes
     Load(..), RangeLoad(..),
     VecLoad(..), RangeVecLoad(..),
+
+    -- * Compute functions
     compute,
     dComputeP, dComputeS,
 
@@ -15,8 +17,11 @@
     L, SH,
 
     -- * Utility
-    entire
+    entire,
 
+    -- * Work index
+    WorkIndex(..), PreferredWorkIndex(..),
+
 ) where
 
 import GHC.Conc
@@ -47,6 +52,37 @@
 {-# INLINE threads #-}
 threads = return
 
+-- | Internal implementation class.
+class (Shape sh, Shape i) => WorkIndex sh i where
+    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
+
+instance Shape sh => WorkIndex sh sh where
+    gindex = index
+    gwrite = write
+    gsize = extent
+    {-# INLINE gindex #-}
+    {-# INLINE gwrite #-}
+    {-# INLINE gsize #-}
+
+#define WI_INT_INST(sh)           \
+instance WorkIndex sh Int where { \
+    gindex = linearIndex;         \
+    gwrite = linearWrite;         \
+    gsize = size . extent;        \
+    {-# INLINE gindex #-};        \
+    {-# INLINE gwrite #-};        \
+    {-# INLINE gsize #-};         \
+}
+
+WI_INT_INST(Dim2)
+WI_INT_INST(Dim3)
+
+-- | Internal implementation class.
+class WorkIndex sh i => PreferredWorkIndex l sh i | l sh -> i where
+
+
 -- | This class abstracts pair of array types,
 -- which could be loaded one to another.
 --
@@ -73,7 +109,8 @@
 -- it should have only 3 parameters: @Load l tl sh@.
 -- But Convoluted ('Data.Yarr.Convolution.Repr.CV') representation is
 -- tightly connected with it's load type.
-class (USource r l sh a, UTarget tr tl sh a, Shape (LoadIndex l tl sh)) =>
+class (USource r l sh a, UTarget tr tl sh a,
+       WorkIndex sh (LoadIndex l tl sh)) =>
         Load r l tr tl sh a where
     -- | Used in @fill@ parameter function.
     -- There are two options for this type to be: @sh@ itself or @Int@.
@@ -100,6 +137,14 @@
           -> UArray r l sh a            -- ^ Source array
           -> UArray tr tl sh a          -- ^ Target array
           -> IO ()
+    loadP fill threads arr tarr = do
+        force arr
+        force tarr
+        !ts <- threads
+        parallel_ ts $
+            makeFork ts zero (gsize arr) (fill (gindex arr) (gwrite tarr))
+        touchArray arr
+        touchArray tarr
 
     -- | /O(n)/ Sequential analog of 'loadP' function.
     -- Loads source to target 'entire'ly.
@@ -111,7 +156,17 @@
           -> UArray r l sh a            -- ^ Source array
           -> UArray tr tl sh a          -- ^ Target array
           -> IO ()
+    loadS fill arr tarr = do
+        force arr
+        force tarr
+        fill (gindex arr) (gwrite tarr) zero (gsize arr)
+        touchArray arr
+        touchArray tarr
 
+    {-# INLINE loadP #-}
+    {-# INLINE loadS #-}
+
+
 -- | Class abstracts pair of arrays which could be loaded in
 -- just specified range of indices.
 --
@@ -120,8 +175,7 @@
 -- cube for 'Dim3'. Thus, it is specified by pair of indices:
 -- \"top-left\" (minimum is 'zero') and \"bottom-right\" (maximum is
 -- @('entire' arr tarr)@) corners.
-class (Load r l tr tl sh a, LoadIndex l tl sh ~ sh) =>
-        RangeLoad r l tr tl sh a where
+class Load r l tr tl sh a => RangeLoad r l tr tl sh a where
 
     -- | /O(n)/ Loads elements from source to target in specified range
     -- in parallel.
@@ -141,6 +195,14 @@
         -> sh                -- ^ Top-left 
         -> sh                -- ^ and bottom-right corners of range to load
         -> IO ()
+    rangeLoadP fill threads arr tarr start end = do
+        force arr
+        force tarr
+        !ts <- threads
+        parallel_ ts $
+            makeFork ts start end (fill (index arr) (write tarr))
+        touchArray arr
+        touchArray tarr
 
     -- | /O(n)/ Sequentially loads elements from source to target in specified range.
     rangeLoadS
@@ -150,8 +212,19 @@
         -> sh                -- ^ Top-left
         -> sh                -- ^ and bottom-right corners of range to load
         -> IO ()
+    rangeLoadS fill arr tarr start end = do
+        force arr
+        force tarr
+        fill (index arr) (write tarr) start end
+        touchArray arr
+        touchArray tarr
 
+    {-# INLINE rangeLoadP #-}
+    {-# INLINE rangeLoadS #-}
 
+
+
+
 -- | Class abstracts /separated in time and space/ loading 'slices' of one array type
 -- to another. Result of running functions with @-Slices-@ infix
 -- /is always identical/ to result of running corresponding function from
@@ -183,8 +256,7 @@
 --  * @e@ - vector element type, common for source and target arrays
 --
 class (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e,
-       Load slr l tslr tl sh e, Shape (LoadIndex l tl sh),
-       Dim v ~ Dim v2) =>
+       Load slr l tslr tl sh e, Dim v ~ Dim v2) =>
         VecLoad r slr l tr tslr tl sh v v2 e where
 
     -- | /O(n)/ Entirely, slice-wise loads vectors from source to target 
@@ -202,6 +274,18 @@
         -> UArray r l sh (v e)        -- ^ Source array of vectors
         -> UArray tr tl sh (v2 e)     -- ^ Target array of vectors
         -> IO ()
+    loadSlicesP fill threads arr tarr = do
+        force arr
+        force tarr
+        !ts <- threads
+        parallel_ ts $
+            makeForkSlicesOnce
+                ts (V.replicate (zero, gsize arr))
+                (V.zipWith
+                    (\sl tsl -> fill (gindex sl) (gwrite tsl))
+                    (slices arr) (slices tarr))
+        touchArray arr
+        touchArray tarr
 
     -- | /O(n)/ Sequentially loads vectors from source to target, slice by slice.
     loadSlicesS
@@ -209,11 +293,21 @@
         -> UArray r l sh (v e)        -- ^ Source array of vectors
         -> UArray tr tl sh (v2 e)     -- ^ Target array of vectors
         -> IO ()
+    loadSlicesS fill arr tarr = do
+        force arr
+        force tarr
+        V.zipWithM_ (loadS fill) (slices arr) (slices tarr)
+        touchArray arr
+        touchArray tarr
 
+    {-# INLINE loadSlicesP #-}
+    {-# INLINE loadSlicesS #-}
+
+
 -- | This class extends 'VecLoad' just like 'RangeLoad' extends 'Load'.
 -- It abstracts slice-wise loading from one array type to
 -- another in specified range.
-class (VecLoad r slr l tr tslr tl sh v v2 e, LoadIndex l tl sh ~ sh) =>
+class (VecLoad r slr l tr tslr tl sh v v2 e, RangeLoad slr l tslr tl sh e) =>
         RangeVecLoad r slr l tr tslr tl sh v v2 e where
 
     -- | /O(n)/ Loads vectors from source to target in specified range, slice-wise,
@@ -226,6 +320,18 @@
         -> sh                     -- ^ Top-left
         -> sh                     -- ^ and bottom-right corners of range to load
         -> IO ()
+    rangeLoadSlicesP fill threads arr tarr start end = do
+        force arr
+        force tarr
+        !ts <- threads
+        parallel_ ts $
+            makeForkSlicesOnce
+                ts (V.replicate (start, end))
+                (V.zipWith
+                    (\sl tsl -> fill (index sl) (write tsl))
+                    (slices arr) (slices tarr))
+        touchArray arr
+        touchArray tarr
 
     -- | /O(n)/ Sequentially loads vector elements from source to target
     -- in specified range, slice by slice.
@@ -236,7 +342,19 @@
         -> sh                     -- ^ Top-left
         -> sh                     -- ^ and bottom-right corners of range to load
         -> IO ()
+    rangeLoadSlicesS fill arr tarr start end = do
+        force arr
+        force tarr
+        V.zipWithM_
+            (\sl tsl -> rangeLoadS fill sl tsl start end)
+            (slices arr) (slices tarr)
+        touchArray arr
+        touchArray tarr
 
+    {-# INLINE rangeLoadSlicesP #-}
+    {-# INLINE rangeLoadSlicesS #-}
+
+
 -- | /O(n)/ This function simplifies the most common way of loading
 -- arrays.
 --
@@ -265,6 +383,9 @@
     load arr marr
     freeze marr
 
+-- | Most common parallel use case of 'compute'.
+--
+-- @dComputeP = 'compute' ('loadP' 'S.fill' 'caps')@
 dComputeP
     :: (USource r l sh a, Manifest tr mtr tl sh a,
         Load r l mtr tl sh a)
@@ -273,6 +394,10 @@
 {-# INLINE dComputeP #-}
 dComputeP = compute (loadP fill caps)
 
+
+-- | Most common sequential use case of 'compute'.
+--
+-- @dComputeS = 'compute' ('loadS' 'S.fill')@
 dComputeS
     :: (USource r l sh a, Manifest tr mtr tl sh a,
         Load r l mtr tl sh a)
@@ -293,55 +418,21 @@
 -- functions defined by default.
 data L
 
-instance (USource r L sh a, UTarget tr L sh a) => Load r L tr L sh a where
+instance WorkIndex sh Int => PreferredWorkIndex L sh Int
 
+instance (USource r L sh a, UTarget tr L sh a, WorkIndex sh Int) =>
+        Load r L tr L sh a where
     type LoadIndex L L sh = Int
-    
-    loadP lfill threads arr tarr = do
-        force arr
-        force tarr
-        !ts <- threads
-        parallel_ ts $
-            makeFork ts 0 (size (extent arr))
-                     (lfill (linearIndex arr) (linearWrite tarr))
-        touchArray arr
-        touchArray tarr
 
-    loadS lfill arr tarr = do
-        force arr
-        force tarr
-        lfill (linearIndex arr) (linearWrite tarr) 0 (size (extent arr))
-        touchArray arr
-        touchArray tarr
-
-    {-# INLINE loadP #-}
-    {-# INLINE loadS #-}
+instance Load r L tr L sh a => RangeLoad r L tr L sh a
 
 instance (UVecSource r slr L sh v e, UVecTarget tr tslr L sh v2 e,
           Load slr L tslr L sh e, Dim v ~ Dim v2) =>
-        VecLoad r slr L tr tslr L sh v v2 e where
-    loadSlicesP lfill threads arr tarr = do
-        force arr
-        force tarr
-        !ts <- threads
-        parallel_ ts $
-            makeForkSlicesOnce
-                ts (V.replicate (0, size (extent arr)))
-                (V.zipWith
-                    (\sl tsl -> lfill (linearIndex sl) (linearWrite tsl))
-                    (slices arr) (slices tarr))
-        touchArray arr
-        touchArray tarr
+        VecLoad r slr L tr tslr L sh v v2 e
 
-    loadSlicesS lfill arr tarr = do
-        force arr
-        force tarr
-        V.zipWithM_ (loadS lfill) (slices arr) (slices tarr)
-        touchArray arr
-        touchArray tarr
+instance (VecLoad r slr L tr tslr L sh v v2 e, RangeLoad slr L tslr L sh e) =>
+        RangeVecLoad r slr L tr tslr L sh v v2 e
 
-    {-# INLINE loadSlicesP #-}
-    {-# INLINE loadSlicesS #-}
 
 -- | General shape load type index. 'UArray's with 'SH' load type index
 -- specialize 'index' and 'write' and leave 'linearIndex' and 'linearWrite'
@@ -354,119 +445,22 @@
 -- Integral division is very expensive operation even on modern CPUs.
 data SH
 
-#define SH_LOAD_INST(l,tl)                                               \
-instance (USource r l sh a, UTarget tr tl sh a) =>                       \
-        Load r l tr tl sh a where {                                      \
-    type LoadIndex l tl sh = sh;                                         \
-    loadP fill threads arr tarr =                                        \
-        shRangeLoadP fill threads arr tarr zero (entire arr tarr);       \
-    loadS fill arr tarr =                                                \
-        shRangeLoadS fill arr tarr zero (entire arr tarr);               \
-    {-# INLINE loadP #-};                                                \
-    {-# INLINE loadS #-};                                                \
-};                                                                       \
-instance (USource r l sh a, UTarget tr tl sh a) =>                       \
-        RangeLoad r l tr tl sh a where {                                 \
-    rangeLoadP = shRangeLoadP;                                           \
-    rangeLoadS = shRangeLoadS;                                           \
-    {-# INLINE rangeLoadP #-};                                           \
-    {-# INLINE rangeLoadS #-};                                           \
-};                                                                       \
-instance (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e,      \
-          Load slr l tslr tl sh e, Dim v ~ Dim v2) =>                    \
-        VecLoad r slr l tr tslr tl sh v v2 e where {                     \
-    loadSlicesP fill threads arr tarr =                                  \
-        shRangeLoadSlicesP fill threads arr tarr zero (entire arr tarr); \
-    loadSlicesS fill arr tarr =                                          \
-        shRangeLoadSlicesS fill arr tarr zero (entire arr tarr);         \
-    {-# INLINE loadSlicesP #-};                                          \
-    {-# INLINE loadSlicesS #-};                                          \
-};                                                                       \
-instance (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e,      \
-          Load slr l tslr tl sh e, Dim v ~ Dim v2) =>                    \
-        RangeVecLoad r slr l tr tslr tl sh v v2 e where {                \
-    rangeLoadSlicesP = shRangeLoadSlicesP;                               \
-    rangeLoadSlicesS = shRangeLoadSlicesS;                               \
-    {-# INLINE rangeLoadSlicesP #-};                                     \
-    {-# INLINE rangeLoadSlicesS #-};                                     \
-}
+instance Shape sh => PreferredWorkIndex SH sh sh
 
+#define SH_LOAD_INST(l,tl)                                          \
+instance (USource r l sh a, UTarget tr tl sh a) =>                  \
+        Load r l tr tl sh a where {                                 \
+    type LoadIndex l tl sh = sh;                                    \
+};                                                                  \
+instance Load r l tr tl sh a => RangeLoad r l tr tl sh a;           \
+instance (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e, \
+          Load slr l tslr tl sh e, Dim v ~ Dim v2) =>               \
+        VecLoad r slr l tr tslr tl sh v v2 e;                       \
+instance (VecLoad r slr l tr tslr tl sh v v2 e,                     \
+          RangeLoad slr l tslr tl sh e) =>                          \
+        RangeVecLoad r slr l tr tslr tl sh v v2 e;                  \
+
 SH_LOAD_INST(SH,L)
 SH_LOAD_INST(L,SH)
 SH_LOAD_INST(SH,SH)
 
-
-shRangeLoadP
-    :: (USource r l sh a, UTarget tr tl sh a)
-    => Fill sh a
-    -> Threads
-    -> UArray r l sh a
-    -> UArray tr tl sh a
-    -> sh -> sh
-    -> IO ()
-{-# INLINE shRangeLoadP #-}
-shRangeLoadP fill threads arr tarr start end = do
-    force arr
-    force tarr
-    !ts <- threads
-    parallel_ ts $
-        makeFork ts start end (fill (index arr) (write tarr))
-    touchArray arr
-    touchArray tarr
-
-shRangeLoadS
-    :: (USource r l sh a, UTarget tr tl sh a)
-    => Fill sh a
-    -> UArray r l sh a
-    -> UArray tr tl sh a
-    -> sh -> sh
-    -> IO ()
-{-# INLINE shRangeLoadS #-}
-shRangeLoadS fill arr tarr start end = do
-    force arr
-    force tarr
-    fill (index arr) (write tarr) start end
-    touchArray arr
-    touchArray tarr
-
-
-shRangeLoadSlicesP
-    :: (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e,
-        Dim v ~ Dim v2)
-    => Fill sh e
-    -> Threads
-    -> UArray r l sh (v e)
-    -> UArray tr tl sh (v2 e)
-    -> sh -> sh
-    -> IO ()
-{-# INLINE shRangeLoadSlicesP #-}
-shRangeLoadSlicesP fill threads arr tarr start end = do
-    force arr
-    force tarr
-    !ts <- threads
-    parallel_ ts $
-        makeForkSlicesOnce
-            ts (V.replicate (start, end))
-            (V.zipWith
-                (\sl tsl -> fill (index sl) (write tsl))
-                (slices arr) (slices tarr))
-    touchArray arr
-    touchArray tarr
-
-shRangeLoadSlicesS
-    :: (UVecSource r slr l sh v e, UVecTarget tr tslr tl sh v2 e,
-        Dim v ~ Dim v2)
-    => Fill sh e
-    -> UArray r l sh (v e)
-    -> UArray tr tl sh (v2 e)
-    -> sh -> sh
-    -> IO ()
-{-# INLINE shRangeLoadSlicesS #-}
-shRangeLoadSlicesS fill arr tarr start end = do
-    force arr
-    force tarr
-    V.zipWithM_
-        (\sl tsl -> shRangeLoadS fill sl tsl start end)
-        (slices arr) (slices tarr)
-    touchArray arr
-    touchArray tarr
diff --git a/Data/Yarr/Fold.hs b/Data/Yarr/Fold.hs
deleted file mode 100644
--- a/Data/Yarr/Fold.hs
+++ /dev/null
@@ -1,190 +0,0 @@
-
-module Data.Yarr.Fold (
-    -- * Fold support
-    Fold,
-    reduceL, reduceLeftM,
-    reduceR, reduceRightM,
-
-    -- * Fold runners
-    runFold, runFoldP,
-    runFoldSlicesSeparate, runFoldSlicesSeparateP,
-
-    -- * Shortcuts
-    toList,
-) where
-
-import Prelude as P
-import Control.Monad as M
-import Data.List (groupBy)
-import Data.Function (on)
-
-import Data.Yarr.Base
-import Data.Yarr.Shape as S
-import Data.Yarr.Eval
-import Data.Yarr.Convolution
-
-import Data.Yarr.Utils.FixedVector as V hiding (toList)
-import Data.Yarr.Utils.Fork
-import Data.Yarr.Utils.Parallel
-
-
-class Shape fsh => Reduce l ash fsh | l ash -> fsh where
-    getIndex :: USource r l ash a => UArray r l ash a -> (fsh -> IO a)
-    getSize :: USource r l ash a => UArray r l ash a -> fsh
-
-#define SH_GET_INDEX(l,sh)      \
-instance Reduce l sh sh where { \
-    getIndex = index;           \
-    getSize = extent;           \
-    {-# INLINE getIndex #-};    \
-    {-# INLINE getSize #-};     \
-}
-
-SH_GET_INDEX(SH, Dim1)
-SH_GET_INDEX(SH, Dim2)
-SH_GET_INDEX(SH, Dim3)
-
-SH_GET_INDEX(CVL, Dim1)
-SH_GET_INDEX(CVL, Dim2)
-SH_GET_INDEX(CVL, Dim3)
-
-
-#define LINEAR_GET_INDEX(l,sh)   \
-instance Reduce l sh Int where { \
-    getIndex = linearIndex;      \
-    getSize = size . extent;     \
-    {-# INLINE getIndex #-};     \
-    {-# INLINE getSize #-};      \
-}
-
-LINEAR_GET_INDEX(L, Dim1)
-LINEAR_GET_INDEX(L, Dim2)
-LINEAR_GET_INDEX(L, Dim3)
-
-
--- | Curried 'Foldl' or 'Foldr'.
--- Generalizes both partially applied left and right folds.
---
--- See source of following 4 functions to construct more similar ones,
--- if you need.
-type Fold sh a b = 
-       IO b         -- ^ Zero
-    -> (sh -> IO a) -- ^ Get
-    -> sh           -- ^ Start
-    -> sh           -- ^ End
-    -> IO b         -- ^ Result
-
--- | /O(0)/
-reduceLeftM
-    :: Foldl sh a b     -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
-    -> (b -> a -> IO b) -- ^ Monaric left reduce
-    -> Fold sh a b      -- ^ Curried fold to be passed to 'runFold' functions.
-{-# INLINE reduceLeftM #-}
-reduceLeftM foldl rf = foldl (\b _ a -> rf b a)
-
--- | /O(0)/
-reduceL
-    :: Foldl sh a b  -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
-    -> (b -> a -> b) -- ^ Pure left reduce
-    -> Fold sh a b   -- ^ Curried fold to be passed to 'runFold' functions.
-{-# INLINE reduceL #-}
-reduceL foldl rf = foldl (\b _ a -> return $ rf b a)
-
--- | /O(0)/
-reduceRightM
-    :: Foldr sh a b     -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
-    -> (a -> b -> IO b) -- ^ Monaric right reduce
-    -> Fold sh a b      -- ^ Curried fold to be passed to 'runFold' functions.
-{-# INLINE reduceRightM #-}
-reduceRightM foldr rf = foldr (\_ a b -> rf a b)
-
--- | /O(0)/
-reduceR
-    :: Foldr sh a b  -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
-    -> (a -> b -> b) -- ^ Pure right reduce
-    -> Fold sh a b   -- ^ Curried fold to be passed to 'runFold' functions.
-{-# INLINE reduceR #-}
-reduceR foldr rf = foldr (\_ a b -> return $ rf a b)
-
-
--- | /O(n)/
---
--- Example:
---
--- @'toList' = runFold ('reduceR' 'S.foldr' (:)) (return [])@
-runFold
-    :: (USource r l sh a, Reduce l sh fsh)
-    => Fold fsh a b     -- ^ Curried folding worker function
-    -> IO b             -- ^ Monadic fold zero. Wrap pure zero in 'return'.
-    -> UArray r l sh a  -- ^ Source array 
-    -> IO b             -- ^ Fold result
-{-# INLINE runFold #-}
-runFold fold mz arr = do
-    force arr
-    res <- fold mz (getIndex arr) zero (getSize arr)
-    touchArray arr
-    return res
-
--- | /O(n)/ Run associative fold in parallel.
---
--- Example -- associative image histogram filling in the test:
--- <https://github.com/leventov/yarr/blob/master/tests/lum-equalization.hs>
-runFoldP
-    :: (USource r l sh a, Reduce l sh fsh)
-    => Threads          -- ^ Number of threads to parallelize folding on
-    -> Fold fsh a b     -- ^ Curried folding worker function
-    -> IO b             -- ^ Monadic fold zero. Wrap pure zero in 'return'.
-    -> (b -> b -> IO b) -- ^ Associative monadic result joining function
-    -> UArray r l sh a  -- ^ Source array
-    -> IO b             -- ^ Fold result
-{-# INLINE runFoldP #-}
-runFoldP threads fold mz join arr = do
-    force arr
-    ts <- threads
-    (r:rs) <- parallel ts $
-                makeFork ts zero (getSize arr) (fold mz (getIndex arr))
-    touchArray arr
-
-    M.foldM join r rs
-
--- | /O(n)/ 
-runFoldSlicesSeparate
-    :: (UVecSource r slr l sh v e, Reduce l sh fsh)
-    => Fold fsh e b           -- ^ Curried folding function to work on slices
-    -> IO b                   -- ^ Monadic fold zero. Wrap pure zero in 'return'.
-    -> UArray r l sh (v e)    -- ^ Source array of vectors
-    -> IO (VecList (Dim v) b) -- ^ Vector of fold results
-{-# INLINE runFoldSlicesSeparate #-}
-runFoldSlicesSeparate fold mz arr =
-    V.mapM (\sl -> runFold fold mz sl) (slices arr)
-
--- | /O(n)/ Run associative fold over slices of array of vectors in parallel.
-runFoldSlicesSeparateP
-    :: (UVecSource r slr l sh v e, Reduce l sh fsh)
-    => Threads                -- ^ Number of threads to parallelize folding on
-    -> Fold fsh e b           -- ^ Curried folding function to work on slices
-    -> IO b                   -- ^ Monadic fold zero. Wrap pure zero in 'return'.
-    -> (b -> b -> IO b)       -- ^ Associative monadic result joining function
-    -> UArray r l sh (v e)    -- ^ Source array of vectors
-    -> IO (VecList (Dim v) b) -- ^ Vector of fold results
-{-# INLINE runFoldSlicesSeparateP #-}
-runFoldSlicesSeparateP threads fold mz join arr = do
-    force arr
-    ts <- threads
-    trs <- parallel ts $
-            makeForkSlicesOnce
-                ts
-                (V.replicate (zero, getSize arr))
-                (V.map (\sl -> fold mz (getIndex sl)) (slices arr))
-    touchArray arr
-
-    let rsBySlices = P.map (P.map snd) $ groupBy ((==) `on` fst) $ concat trs
-    rs <- M.mapM (\(r:rs) -> M.foldM join r rs) rsBySlices
-    return (VecList rs)
-
--- | /O(n)/ Covert array to list.
-toList
-    :: (USource r l sh a, Reduce l sh fsh)
-    => UArray r l sh a
-    -> IO [a]
-toList = runFold (reduceR S.foldr (:)) (return [])
diff --git a/Data/Yarr/IO/List.hs b/Data/Yarr/IO/List.hs
new file mode 100644
--- /dev/null
+++ b/Data/Yarr/IO/List.hs
@@ -0,0 +1,40 @@
+
+module Data.Yarr.IO.List where
+
+import Control.Monad
+
+import Data.Yarr.Base
+import Data.Yarr.Shape as S
+import Data.Yarr.Eval
+import Data.Yarr.Work
+
+import Debug.Yarr
+
+-- | /O(n)/ Covert array to flat list.
+-- Multidimentional arrays are flatten in column-major order:
+-- 
+-- \[(elem at (0, .., 0, 1)), (elem at (0, .., 0, 2)), ...\]
+toList
+    :: (USource r l sh a, PreferredWorkIndex l sh i)
+    => UArray r l sh a -> IO [a]
+{-# INLINE toList #-}
+toList = work (reduceR S.foldr (:)) (return [])
+
+-- | /O(n)/ Loads manifest array into memory, with elements
+-- from flatten list.
+--
+-- Use this function in the last resort, there are plenty of
+-- methods to 'Load' array, from 'Data.Yarr.D'elayed array for example.
+fromList
+    :: Manifest r mr l sh a
+    => sh                   -- ^ Extent of array
+    -> [a]                  -- ^ Flatten elements
+    -> IO (UArray r l sh a) -- ^ Result manifest array
+{-# INLINE fromList #-}
+fromList sh xs =
+    if (length xs) /= (size sh)
+        then yerr "fromList: list length doesn't correspond size of array shape"
+        else do
+            arr <- new sh
+            zipWithM_ (linearWrite arr) [0..] xs
+            freeze arr
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,7 @@
     UArray(..),
 
     -- * Misc
-    L, SH, delay, delayShaped,
+    L, SH, fromFunction, delay, delayShaped,
 ) where
 
 import Prelude as P
@@ -223,6 +223,21 @@
       => UArray r l sh a -> UArray D l sh a
 {-# INLINE delay #-}
 delay = B.fmap id
+
+-- | Wrap indexing function into delayed representation.
+-- 
+-- Use this function carefully, don't implement through it something
+-- 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.
+fromFunction
+    :: Shape sh
+    => sh               -- ^ Extent of array
+    -> (sh -> IO a)     -- ^ Indexing function
+    -> UArray D SH sh a -- ^ Result array
+{-# INLINE fromFunction #-}
+fromFunction sh f = ShapeDelayed sh (return ()) (return ()) f
 
 -- | Wraps @('index' arr)@ into Delayed representation. Normally you shouldn't need
 -- to use this function. It may be dangerous for performance, because
diff --git a/Data/Yarr/Shape.hs b/Data/Yarr/Shape.hs
--- a/Data/Yarr/Shape.hs
+++ b/Data/Yarr/Shape.hs
@@ -1,44 +1,32 @@
 {-# LANGUAGE InstanceSigs #-}
 
-module Data.Yarr.Shape where
+module Data.Yarr.Shape (
+    -- * Flow types hierarchy
+    module Data.Yarr.WorkTypes,
 
+    -- * Shape and BlockShape
+    Block, Shape(..), BlockShape(..),
+
+    -- * Shape instances
+    Dim1, Dim2, Dim3,
+
+    -- * Specialized flow
+    dim2BlockFill,
+) where
+
 import Prelude as P hiding (foldl, foldr)
 import GHC.Exts
 
 import Control.DeepSeq
 
+import Data.Yarr.WorkTypes
+
 import Data.Yarr.Utils.FixedVector as V hiding (foldl, foldr)
 import Data.Yarr.Utils.LowLevelFlow
 import Data.Yarr.Utils.Primitive
 import Data.Yarr.Utils.Split
 
--- | Alias to frequently used get-write-from-to arguments combo.
---
--- Passed as 1st parameter of all 'Data.Yarr.Eval.Load'ing functions
--- from "Data.Yarr.Eval" module.
-type Fill sh a =
-       (sh -> IO a)       -- ^ Get
-    -> (sh -> a -> IO ()) -- ^ Write
-    -> sh                 -- ^ Start
-    -> sh                 -- ^ End
-    -> IO ()
 
-type Foldl sh a b =
-       (b -> sh -> a -> IO b) -- ^ Generalized left reduce
-    -> IO b                   -- ^ Zero
-    -> (sh -> IO a)           -- ^ Get
-    -> sh                     -- ^ Start
-    -> sh                     -- ^ End
-    -> IO b                   -- ^ Result
-
-type Foldr sh a b =
-       (sh -> a -> b -> IO b) -- ^ Generalized right reduce
-    -> IO b                   -- ^ Zero
-    -> (sh -> IO a)           -- ^ Get
-    -> sh                     -- ^ Start
-    -> sh                     -- ^ End
-    -> IO b                   -- ^ Result
-
 -- | Mainly for internal use.
 -- Abstracts top-left -- bottom-right pair of indices.
 type Block sh = (sh, sh)
@@ -96,27 +84,32 @@
     blockSize :: Block sh -> Int
     insideBlock :: Block sh -> sh -> Bool
 
-    -- | Following 6 functions shouldn't be called directly,
-    -- they are intented to be passed as first argument
-    -- to 'Data.Yarr.Eval.Load' and functions from
-    -- "Data.Yarr.Fold" module.
     makeChunkRange :: Int -> sh -> sh -> (Int -> Block sh)
 
+    -- | Standard left fold wothout unrolling.
+    --
+    -- This one and 5 following functions shouldn't be called directly,
+    -- they are intented to be passed as first argument
+    -- to 'Data.Yarr.Eval.Load' and functions from
+    -- "Data.Yarr.Work" module.
     foldl :: Foldl sh a b
 
     unrolledFoldl
         :: forall a b uf. Arity uf
-        => uf                     -- ^ Unroll factor
-        -> (a -> IO ())           -- ^ 'touch' or 'noTouch'
-        -> Foldl sh a b
+        => uf           -- ^ Unroll factor
+        -> (a -> IO ()) -- ^ 'touch' or 'noTouch'
+        -> Foldl sh a b -- ^ Result curried function
+                        -- to be passed to working functions
 
+    -- | Standard right folding function without unrolling.
     foldr :: Foldr sh a b
 
     unrolledFoldr
         :: forall a b uf. Arity uf
-        => uf                     -- ^ Unroll factor
-        -> (a -> IO ())           -- ^ 'touch' or 'noTouch'
-        -> Foldr sh a b
+        => uf           -- ^ Unroll factor
+        -> (a -> IO ()) -- ^ 'touch' or 'noTouch'
+        -> Foldr sh a b -- ^ Result curried function
+                        -- to be passed to working functions
 
     -- | Standard fill without unrolling.
     -- To avoid premature optimization just type @fill@
@@ -126,10 +119,10 @@
 
     unrolledFill
         :: forall a uf. Arity uf
-        => uf                 -- ^ Unroll factor
-        -> (a -> IO ())       -- ^ 'touch' or 'noTouch'
-        -> Fill sh a          -- ^ Result curried function
-                              --   to pass to loading functions.
+        => uf           -- ^ Unroll factor
+        -> (a -> IO ()) -- ^ 'touch' or 'noTouch'
+        -> Fill sh a    -- ^ Result curried function
+                        -- to by passed to loading functions
 
     {-# INLINE minus #-}
     {-# INLINE intersectBlocks #-}
@@ -382,19 +375,19 @@
 
     {-# INLINE clipBlock #-}
 
--- | 2D-unrolling to maximize profit from
+-- | 2D-unrolled filling to maximize profit from
 -- \"Global value numbering\" LLVM optimization.
 --
 -- Example:
 --
--- @blurred <- 'Data.Yarr.Eval.compute' ('Data.Yarr.Eval.loadP' (dim2BlockFill 'n1' 'n4' 'touch')) delayedBlurred@
+-- @blurred <- 'Data.Yarr.Eval.compute' ('Data.Yarr.Eval.loadS' (dim2BlockFill 'n1' 'n4' 'touch')) delayedBlurred@
 dim2BlockFill
     :: forall a bsx bsy. (Arity bsx, Arity bsy)
     => bsx                  -- ^ Block size by x. Use 'n1' - 'n8' values.
     -> bsy                  -- ^ Block size by y
     -> (a -> IO ())         -- ^ 'touch' or 'noTouch'
     -> Fill Dim2 a          -- ^ Result curried function
-                            --   to pass to loading functions.
+                            --   to be passed to loading functions.
 {-# INLINE dim2BlockFill #-}
 dim2BlockFill blockSizeX blockSizeY tch =
     \get write ->
@@ -583,5 +576,3 @@
     {-# INLINE unrolledFoldl #-}
     {-# INLINE foldr #-}
     {-# INLINE unrolledFoldr #-}
-
-
diff --git a/Data/Yarr/Utils/Fork.hs b/Data/Yarr/Utils/Fork.hs
--- a/Data/Yarr/Utils/Fork.hs
+++ b/Data/Yarr/Utils/Fork.hs
@@ -10,10 +10,11 @@
 
 makeForkEachSlice
     :: (Shape sh, Arity n, v ~ VecList n)
-    => Int
-    -> sh -> sh
-    -> v (sh -> sh -> IO a)
-    -> (Int -> IO (v a))
+    => Int               -- ^ Number of threads to fork work on
+    -> sh                -- ^ Start
+    -> sh                -- ^ End
+    -> v (Work sh a)     -- ^ Slice works
+    -> (Int -> IO (v a)) -- ^ Thread work, returns piece of result for each slice
 {-# INLINE makeForkEachSlice #-}
 makeForkEachSlice threads start end rangeWorks =
     let {-# INLINE etWork #-}
@@ -23,10 +24,11 @@
 
 makeForkSlicesOnce
     :: (Shape sh, Arity n)
-    => Int
-    -> VecList n (sh, sh)
-    -> VecList n (sh -> sh -> IO a)
-    -> (Int -> IO [(Int, a)])
+    => Int                    -- ^ Number of threads to fork work on
+    -> VecList n (sh, sh)     -- ^ (start, end) for each slice
+    -> VecList n (Work sh a)  -- ^ Slice works
+    -> (Int -> IO [(Int, a)]) -- ^ Thread work, returns pieces of results:
+                              -- [(slice number, result)]
 {-# INLINE makeForkSlicesOnce #-}
 makeForkSlicesOnce !threads ranges rangeWorks =
     let !slices = V.length rangeWorks
@@ -72,9 +74,11 @@
 
 makeFork
     :: Shape sh
-    => Int
-    -> sh -> sh
-    -> ((sh -> sh -> IO a) -> (Int -> IO a))
+    => Int            -- ^ Number of threads to fork work on
+    -> sh             -- ^ Start
+    -> sh             -- ^ End
+    -> (Work sh a)    -- ^ Work
+    -> (Int -> IO a)  -- ^ Thread work
 {-# INLINE makeFork #-}
 makeFork chunks start end =
     let {-# INLINE chunkRange #-}
diff --git a/Data/Yarr/Utils/LowLevelFlow.hs b/Data/Yarr/Utils/LowLevelFlow.hs
--- a/Data/Yarr/Utils/LowLevelFlow.hs
+++ b/Data/Yarr/Utils/LowLevelFlow.hs
@@ -56,14 +56,15 @@
 {-# INLINE foldl# #-}
 foldl# reduce mz get start# end# =
     let {-# INLINE go# #-}
-        go# i# b
+        go# !b i#
             | i# >=# end# = return b
             | otherwise   = do
                 let i = (I# i#)
                 a <- get i
                 b' <- reduce b i a
-                go# (i# +# 1#) b'
-    in mz >>= go# start#
+                go# b' (i# +# 1#)
+    in do z <- mz
+          go# z start#
 
 unrolledFoldl#
     :: forall a b uf. Arity uf
@@ -79,8 +80,8 @@
     let !(I# uf#) = arity unrollFactor
         lim# = end# -# uf#
         {-# INLINE go# #-}
-        go# i# b
-            | i# ># lim# = rest# i# b
+        go# !b i#
+            | i# ># lim# = rest# b i#
             | otherwise  = do
                 let is :: VecList uf Int
                     is = V.generate (+ (I# i#))
@@ -89,19 +90,20 @@
                 b' <- V.foldM
                         (\b (i, a) -> reduce b i a) b
                         (V.zipWith (,) is as)
-                go# (i# +# uf#) b'
+                go# b' (i# +# uf#)
 
         {-# INLINE rest# #-}
-        rest# i# b
+        rest# !b i#
             | i# >=# end# = return b
             | otherwise   = do
                 let i = (I# i#)
                 a <- get i
                 tch a
                 b' <- reduce b i a
-                rest# (i# +# 1#) b'
+                rest# b' (i# +# 1#)
 
-    in mz >>= go# start#
+    in do z <- mz
+          go# z start#
 
 
 foldr#
@@ -113,14 +115,15 @@
 {-# INLINE foldr# #-}
 foldr# reduce mz get start# end# =
     let {-# INLINE go# #-}
-        go# i# b
+        go# !b i#
             | i# <# start# = return b
             | otherwise    = do
                 let i = (I# i#)
                 a <- get i
                 b' <- reduce i a b
-                go# (i# -# 1#) b'
-    in mz >>= go# (end# -# 1#)
+                go# b' (i# -# 1#)
+    in do z <- mz
+          go# z (end# -# 1#)
 
 unrolledFoldr#
     :: forall a b uf. Arity uf
@@ -136,8 +139,8 @@
     let !(I# uf#) = arity unrollFactor
         lim# = start# +# uf# -# 1#
         {-# INLINE go# #-}
-        go# i# b
-            | i# <# lim# = rest# i# b
+        go# !b i#
+            | i# <# lim# = rest# b i#
             | otherwise  = do
                 let is :: VecList uf Int
                     is = V.generate ((I# i#) -)
@@ -146,17 +149,18 @@
                 b' <- V.foldM
                         (\b (i, a) -> reduce i a b) b
                         (V.zipWith (,) is as)
-                go# (i# -# uf#) b'
+                go# b' (i# -# uf#)
 
         {-# INLINE rest# #-}
-        rest# i# b
+        rest# !b i#
             | i# <# start# = return b
             | otherwise    = do
                 let i = (I# i#)
                 a <- get i
                 tch a
                 b' <- reduce i a b
-                rest# (i# -# 1#) b'
+                rest# b' (i# -# 1#)
 
-    in mz >>= go# (end# -# 1#)
+    in do z <- mz
+          go# z (end# -# 1#)
 
diff --git a/Data/Yarr/Utils/Split.hs b/Data/Yarr/Utils/Split.hs
--- a/Data/Yarr/Utils/Split.hs
+++ b/Data/Yarr/Utils/Split.hs
@@ -2,9 +2,10 @@
 module Data.Yarr.Utils.Split where
 
 makeSplitIndex
-    :: Int
-    -> Int -> Int
-    -> (Int -> Int)
+    :: Int          -- ^ Number of chunks to split range on (@n@)
+    -> Int          -- ^ Start of range
+    -> Int          -- ^ End of range
+    -> (Int -> Int) -- ^ Split index function
 {-# INLINE makeSplitIndex #-}
 makeSplitIndex chunks start end =
     let !len = end - start
@@ -14,8 +15,11 @@
                  in \c -> if c < chunkLeftover
                         then start + c * (chunkLen + 1)
                         else start + c * chunkLen + chunkLeftover
-
-evenChunks :: [a] -> Int -> [[a]]
+-- | Well-known missed in "Data.List.Split" function.
+evenChunks
+    :: [a]    -- ^ List to split
+    -> Int    -- ^ Number of chuncks (@n@)
+    -> [[a]]  -- ^ Exactly @n@ even chunks of the initial list
 {-# INLINE evenChunks #-}
 evenChunks xs n =
     let len = length xs
diff --git a/Data/Yarr/Work.hs b/Data/Yarr/Work.hs
new file mode 100644
--- /dev/null
+++ b/Data/Yarr/Work.hs
@@ -0,0 +1,270 @@
+
+module Data.Yarr.Work (
+    -- * Fold combinators
+    -- | See source of these 4 functions
+    -- to construct more similar ones,
+    -- if you need.
+    reduceL, reduceLeftM,
+    reduceR, reduceRightM,
+
+    -- * Combinators to work with mutable state
+    -- | Added specially to improve performance
+    -- of tasks like histogram filling.
+    --
+    -- Unfortunately, GHC doesn't figure that folding state
+    -- isn't changed as ADT in such cases and doesn't lift
+    -- it's evaluation higher from folding routine.
+    mutate, imutate,
+
+    -- * Work runners
+    work, iwork, rangeWork,
+    workP, iworkP, rangeWorkP,
+    workOnSlicesSeparate, iworkOnSlicesSeparate, rangeWorkOnSlicesSeparate,
+    workOnSlicesSeparateP, iworkOnSlicesSeparateP, rangeWorkOnSlicesSeparateP,
+
+    -- * Aliases for work types
+    StatefulWork, Foldl, Foldr,
+) where
+
+import Data.Yarr.Base
+import Data.Yarr.Shape as S
+import Data.Yarr.Eval
+
+import Data.Yarr.Work.Internal
+
+
+-- | /O(0)/
+reduceLeftM
+    :: Foldl i a b        -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
+    -> (b -> a -> IO b)   -- ^ Monadic left reduce
+    -> StatefulWork i a b -- ^ Result stateful work to be passed
+                          -- to work runners
+{-# INLINE reduceLeftM #-}
+reduceLeftM foldl rf = foldl (\b _ a -> rf b a)
+
+-- | /O(0)/
+reduceL
+    :: Foldl i a b        -- ^ 'S.foldl' or curried 'S.unrolledFoldl'
+    -> (b -> a -> b)      -- ^ Pure left reduce
+    -> StatefulWork i a b -- ^ Result stateful work to be passed
+                          -- to work runners
+{-# INLINE reduceL #-}
+reduceL foldl rf = foldl (\b _ a -> return $ rf b a)
+
+-- | /O(0)/
+reduceRightM
+    :: Foldr i a b         -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
+    -> (a -> b -> IO b)    -- ^ Monadic right reduce
+    -> StatefulWork i a b  -- ^ Result stateful work to be passed
+                           -- to work runners
+{-# INLINE reduceRightM #-}
+reduceRightM foldr rf = foldr (\_ a b -> rf a b)
+
+-- | /O(0)/
+reduceR
+    :: Foldr i a b        -- ^ 'S.foldr' or curried 'S.unrolledFoldr'
+    -> (a -> b -> b)      -- ^ Pure right reduce
+    -> StatefulWork i a b -- ^ Result stateful work to be passed
+                          -- to work runners
+{-# INLINE reduceR #-}
+reduceR foldr rf = foldr (\_ a b -> return $ rf a b)
+
+
+-- | /O(0)/
+mutate
+    :: Fill i a           -- ^ 'S.fill' or curried 'S.unrolledFill'.
+                          -- If mutating is associative,
+                          -- 'S.dim2BlockFill' is also acceptable.
+    -> (s -> a -> IO ())  -- ^ (state -> array element -> (state has changed))
+                          -- -- State mutating function
+    -> StatefulWork i a s -- ^ Result stateful work to be passed
+                          -- to work runners
+{-# INLINE mutate #-}
+mutate fill mf = imutate fill (\s i -> mf s)
+
+-- | /O(0)/ Version of 'mutate', accepts mutating function
+-- which additionaly accepts array index.
+imutate
+    :: Fill i a               -- ^ 'S.fill' or curried 'S.unrolledFill'.
+                              -- If mutating is associative,
+                              -- 'S.dim2BlockFill' is also acceptable.
+    -> (s -> i -> a -> IO ()) -- ^ Indexed state mutating function
+    -> StatefulWork i a s     -- ^ Result stateful work to be passed
+                              -- to work runners
+{-# INLINE imutate #-}
+imutate fill imf ms index start end = do
+    s <- ms
+    fill index (imf s) start end
+    return s
+
+
+
+-- | /O(n)/ Run non-indexed stateful work.
+--
+-- Example:
+--
+-- @'Data.Yarr.IO.List.toList' = work ('reduceR' 'S.foldr' (:)) (return [])@
+work
+    :: (USource r l sh a, PreferredWorkIndex l sh i)
+    => StatefulWork i a s -- ^ Stateful working function
+    -> IO s               -- ^ Monadic initial state (fold zero).
+                          -- Wrap pure state in 'return'.
+    -> UArray r l sh a    -- ^ Source array
+    -> IO s               -- ^ Final state (fold result)
+{-# INLINE work #-}
+work = anyWork
+
+-- | /O(n)/ Run indexed stateful work.
+--
+-- Example:
+--
+-- @res \<- iwork ('S.foldl' (\\s i a -> ...)) foldZero sourceArray@
+iwork
+    :: USource r l sh a
+    => StatefulWork sh a s -- ^ Stateful working function
+    -> IO s                -- ^ Monadic initial state (fold zero).
+                           -- Wrap pure state in 'return'.
+    -> UArray r l sh a     -- ^ Source array
+    -> IO s                -- ^ Final state (fold result)
+{-# INLINE iwork #-}
+iwork = anyWork
+
+-- | /O(n)/ Run stateful work in specified range of indices.
+rangeWork
+    :: USource r l sh a
+    => StatefulWork sh a s -- ^ Stateful working function
+    -> IO s                -- ^ Monadic initial state (fold zero).
+                           -- Wrap pure state in 'return'.
+    -> UArray r l sh a     -- ^ Source array
+    -> sh                  -- ^ Top-left
+    -> sh                  -- ^ and bottom-right corners of range to work in
+    -> IO s                -- ^ Final state (fold result)
+{-# INLINE rangeWork #-}
+rangeWork = anyRangeWork
+
+
+-- | /O(n)/ Run associative non-indexed stateful work in parallel.
+--
+-- Example -- associative image histogram filling in the test:
+-- <https://github.com/leventov/yarr/blob/master/tests/lum-equalization.hs>
+workP
+    :: (USource r l sh a, PreferredWorkIndex l sh i)
+    => Threads            -- ^ Number of threads to parallelize work on
+    -> StatefulWork i a s -- ^ Associative stateful working function
+    -> IO s               -- ^ Monadic zero state.
+                          -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)   -- ^ Associative monadic state joining function
+    -> UArray r l sh a    -- ^ Source array
+    -> IO s               -- ^ Gathered state (fold result)
+{-# INLINE workP #-}
+workP = anyWorkP
+
+-- | /O(n)/ Run associative indexed stateful work in parallel.
+iworkP
+    :: USource r l sh a
+    => Threads             -- ^ Number of threads to parallelize work on
+    -> StatefulWork sh a s -- ^ Associative stateful working function
+    -> IO s                -- ^ Monadic zero state.
+                           -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)    -- ^ Associative monadic state joining function
+    -> UArray r l sh a     -- ^ Source array
+    -> IO s                -- ^ Gathered state (fold result)
+{-# INLINE iworkP #-}
+iworkP = anyWorkP
+
+-- | /O(n)/ Run associative stateful work in specified range in parallel.
+rangeWorkP
+    :: USource r l sh a
+    => Threads             -- ^ Number of threads to parallelize work on
+    -> StatefulWork sh a s -- ^ Associative stateful working function
+    -> IO s                -- ^ Monadic zero state.
+                           -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)    -- ^ Associative monadic state joining function
+    -> UArray r l sh a     -- ^ Source array
+    -> sh                  -- ^ Top-left
+    -> sh                  -- ^ and bottom-right corners of range to work in
+    -> IO s                -- ^ Gathered state (fold result)
+{-# INLINE rangeWorkP #-}
+rangeWorkP = anyRangeWorkP
+
+
+-- | /O(n)/ Run non-indexed stateful work over each slice of array of vectors.
+workOnSlicesSeparate
+    :: (UVecSource r slr l sh v e, PreferredWorkIndex l sh i)
+    => StatefulWork i e s     -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic initial state (fold zero).
+                              -- Wrap pure state in 'return'.
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> IO (VecList (Dim v) s) -- ^ Vector of final states (fold results)
+{-# INLINE workOnSlicesSeparate #-}
+workOnSlicesSeparate = anyWorkOnSlicesSeparate
+
+-- | /O(n)/ Run indexed stateful work over each slice of array of vectors.
+iworkOnSlicesSeparate
+    :: UVecSource r slr l sh v e
+    => StatefulWork sh e s    -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic initial state (fold zero).
+                              -- Wrap pure state in 'return'.
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> IO (VecList (Dim v) s) -- ^ Vector of final states (fold results)
+{-# INLINE iworkOnSlicesSeparate #-}
+iworkOnSlicesSeparate = anyWorkOnSlicesSeparate
+
+-- | /O(n)/ Run stateful work in specified range
+-- over each slice of array of vectors.
+rangeWorkOnSlicesSeparate
+    :: UVecSource r slr l sh v e
+    => StatefulWork sh e s    -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic initial state (fold zero).
+                              -- Wrap pure state in 'return'.
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> sh                     -- ^ Top-left
+    -> sh                     -- ^ and bottom-right corners of range to work in
+    -> IO (VecList (Dim v) s) -- ^ Vector of final states (fold results)
+{-# INLINE rangeWorkOnSlicesSeparate #-}
+rangeWorkOnSlicesSeparate = anyRangeWorkOnSlicesSeparate
+
+
+-- | /O(n)/ Run associative non-indexed stateful work
+-- over slices of array of vectors in parallel.
+workOnSlicesSeparateP
+    :: (UVecSource r slr l sh v e, PreferredWorkIndex l sh i)
+    => Threads                -- ^ Number of threads to parallelize work on
+    -> StatefulWork i e s     -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic zero state.
+                              -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)       -- ^ Associative monadic state joining function
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> IO (VecList (Dim v) s) -- ^ Vector of gathered per slice results
+{-# INLINE workOnSlicesSeparateP #-}
+workOnSlicesSeparateP = anyWorkOnSlicesSeparateP
+
+-- | /O(n)/ Run associative indexed stateful work
+-- over slices of array of vectors in parallel.
+iworkOnSlicesSeparateP
+    :: UVecSource r slr l sh v e
+    => Threads                -- ^ Number of threads to parallelize work on
+    -> StatefulWork sh e s    -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic zero state.
+                              -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)       -- ^ Associative monadic state joining function
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> IO (VecList (Dim v) s) -- ^ Vector of gathered per slice results
+{-# INLINE iworkOnSlicesSeparateP #-}
+iworkOnSlicesSeparateP = anyWorkOnSlicesSeparateP
+
+-- | /O(n)/ Run associative stateful work in specified range
+-- over slices of array of vectors in parallel.
+rangeWorkOnSlicesSeparateP
+    :: UVecSource r slr l sh v e
+    => Threads                -- ^ Number of threads to parallelize work on
+    -> StatefulWork sh e s    -- ^ Stateful slice-wise working function
+    -> IO s                   -- ^ Monadic zero state.
+                              -- Wrap pure state in 'return'.
+    -> (s -> s -> IO s)       -- ^ Associative monadic state joining function
+    -> UArray r l sh (v e)    -- ^ Source array of vectors
+    -> sh                     -- ^ Top-left
+    -> sh                     -- ^ and bottom-right corners of range to work in
+    -> IO (VecList (Dim v) s) -- ^ Vector of gathered per slice results
+{-# INLINE rangeWorkOnSlicesSeparateP #-}
+rangeWorkOnSlicesSeparateP = anyRangeWorkOnSlicesSeparateP
diff --git a/Data/Yarr/Work/Internal.hs b/Data/Yarr/Work/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Data/Yarr/Work/Internal.hs
@@ -0,0 +1,135 @@
+
+module Data.Yarr.Work.Internal where
+
+import Prelude as P
+import Control.Monad as M
+import Data.List (groupBy)
+import Data.Function (on)
+
+import Data.Yarr.Base
+import Data.Yarr.Shape as S
+import Data.Yarr.Eval
+
+import Data.Yarr.Utils.FixedVector as V hiding (toList)
+import Data.Yarr.Utils.Fork
+import Data.Yarr.Utils.Parallel
+
+
+anyWork
+    :: (USource r l sh a, WorkIndex sh i)
+    => StatefulWork i a s
+    -> IO s
+    -> UArray r l sh a
+    -> IO s
+{-# INLINE anyWork #-}
+anyWork fold mz arr = anyRangeWork fold mz arr zero (gsize arr)
+
+anyRangeWork
+    :: (USource r l sh a, WorkIndex sh i)
+    => StatefulWork i a s
+    -> IO s
+    -> UArray r l sh a
+    -> i -> i
+    -> IO s
+{-# INLINE anyRangeWork #-}
+anyRangeWork fold mz arr start end = do
+    force arr
+    res <- fold mz (gindex arr) start end
+    touchArray arr
+    return res
+
+
+anyWorkP
+    :: (USource r l sh a, WorkIndex sh i)
+    => Threads
+    -> StatefulWork i a s
+    -> IO s
+    -> (s -> s -> IO s)
+    -> UArray r l sh a
+    -> IO s
+{-# INLINE anyWorkP #-}
+anyWorkP threads fold mz join arr =
+    anyRangeWorkP threads fold mz join arr zero (gsize arr)
+
+anyRangeWorkP
+    :: (USource r l sh a, WorkIndex sh i)
+    => Threads
+    -> StatefulWork i a s
+    -> IO s
+    -> (s -> s -> IO s)
+    -> UArray r l sh a
+    -> i -> i
+    -> IO s
+{-# INLINE anyRangeWorkP #-}
+anyRangeWorkP threads fold mz join arr start end = do
+    force arr
+    ts <- threads
+    (r:rs) <- parallel ts $
+                makeFork ts start end (fold mz (gindex arr))
+    touchArray arr
+
+    M.foldM join r rs
+
+
+anyWorkOnSlicesSeparate
+    :: (UVecSource r slr l sh v e, WorkIndex sh i)
+    => StatefulWork i e s
+    -> IO s
+    -> UArray r l sh (v e)
+    -> IO (VecList (Dim v) s)
+{-# INLINE anyWorkOnSlicesSeparate #-}
+anyWorkOnSlicesSeparate fold mz arr =
+    anyRangeWorkOnSlicesSeparate fold mz arr zero (gsize arr)
+
+anyRangeWorkOnSlicesSeparate
+    :: (UVecSource r slr l sh v e, WorkIndex sh i)
+    => StatefulWork i e s
+    -> IO s
+    -> UArray r l sh (v e)
+    -> i -> i
+    -> IO (VecList (Dim v) s)
+{-# INLINE anyRangeWorkOnSlicesSeparate #-}
+anyRangeWorkOnSlicesSeparate fold mz arr start end = do
+    force arr
+    rs <- V.mapM (\sl -> anyRangeWork fold mz sl start end) (slices arr)
+    touchArray arr
+    return rs
+
+anyWorkOnSlicesSeparateP
+    :: (UVecSource r slr l sh v e, WorkIndex sh i)
+    => Threads
+    -> StatefulWork i e s
+    -> IO s
+    -> (s -> s -> IO s)
+    -> UArray r l sh (v e)
+    -> IO (VecList (Dim v) s)
+{-# INLINE anyWorkOnSlicesSeparateP #-}
+anyWorkOnSlicesSeparateP threads fold mz join arr =
+    anyRangeWorkOnSlicesSeparateP threads fold mz join arr zero (gsize arr)
+
+anyRangeWorkOnSlicesSeparateP
+    :: (UVecSource r slr l sh v e, WorkIndex sh i)
+    => Threads
+    -> StatefulWork i e s
+    -> IO s
+    -> (s -> s -> IO s)
+    -> UArray r l sh (v e)
+    -> i -> i
+    -> IO (VecList (Dim v) s)
+{-# INLINE anyRangeWorkOnSlicesSeparateP #-}
+anyRangeWorkOnSlicesSeparateP threads fold mz join arr start end = do
+    force arr
+    let sls = slices arr
+    V.mapM force sls
+
+    ts <- threads
+    trs <- parallel ts $
+            makeForkSlicesOnce
+                ts
+                (V.replicate (start, end))
+                (V.map (\sl -> fold mz (gindex sl)) sls)
+    touchArray arr
+
+    let rsBySlices = P.map (P.map snd) $ groupBy ((==) `on` fst) $ concat trs
+    rs <- M.mapM (\(r:rs) -> M.foldM join r rs) rsBySlices
+    return (VecList rs)
diff --git a/Data/Yarr/WorkTypes.hs b/Data/Yarr/WorkTypes.hs
new file mode 100644
--- /dev/null
+++ b/Data/Yarr/WorkTypes.hs
@@ -0,0 +1,45 @@
+
+module Data.Yarr.WorkTypes where
+
+-- | Generalizes interval works: 'Fill's, 'StatefulWork's.
+--
+-- To be passed to functions from "Data.Yarr.Utils.Fork" module
+-- and called directly.
+type Work sh a =
+       sh   -- ^ Start (lower index)
+    -> sh   -- ^ End (higher index)
+    -> IO a -- ^ Result
+
+-- | Alias to frequently used get-write-from-to arguments combo.
+--
+-- To be passed as 1st parameter of all 'Data.Yarr.Eval.Load'ing functions
+-- from "Data.Yarr.Eval" module.
+type Fill sh a =
+       (sh -> IO a)       -- ^ Indexing function
+    -> (sh -> a -> IO ()) -- ^ Writing function
+    -> Work sh ()         -- ^ Curried result function -- worker
+
+
+-- | Generalizes both partially applied left and right folds,
+-- as well as works on mutable state.
+--
+-- To be passed to fold runners from "Data.Yarr.Work" module.
+type StatefulWork sh a s = 
+       IO s         -- ^ Initial state
+    -> (sh -> IO a) -- ^ Indexing function
+    -> Work sh s    -- ^ Curried result function -- worker,
+                    -- emits final state
+
+-- | Generalizes left to right folds.
+--
+-- To be passed to fold combinators from "Data.Yarr.Work" module.
+type Foldl sh a b =
+       (b -> sh -> a -> IO b) -- ^ Generalized left reduce
+    -> StatefulWork sh a b    -- ^ Curried result stateful work
+
+-- | Generalizes right to left folds.
+--
+-- To be passed to fold combinators from "Data.Yarr.Work" module.
+type Foldr sh a b =
+       (sh -> a -> b -> IO b) -- ^ Generalized right reduce
+    -> StatefulWork sh a b    -- ^ Curried result stateful work
diff --git a/yarr.cabal b/yarr.cabal
--- a/yarr.cabal
+++ b/yarr.cabal
@@ -1,5 +1,5 @@
 Name:                yarr
-Version:             0.9.2
+Version:             1.2.3
 Synopsis:            Yet another array library
 Description:
     Yarr is a new blazing fast dataflow framework (array library),
@@ -15,14 +15,14 @@
     .
     > let greyImage = zipElems (\r g b -> 0.21 * r + 0.71 * g + 0.07 * b) image
     .
-    The library is considerably faster than @repa@.
+    In some cases the library is considerably faster than @repa@.
     See benchmark results: <https://github.com/leventov/yarr/blob/master/tests/bench-results.md>
     .
     Shortcoming by design: lack of pure indexing interface.
     .
     /Changes in version 0.9.2:/
     .
-        * Safe folds -- see "Data.Yarr.Fold"
+        * Safe folds -- see "Data.Yarr.Work"
     .
         * Issue with slice-wise loading with unrolled filling function solved
     .
@@ -69,13 +69,14 @@
         Data.Yarr.Base
         Data.Yarr.Eval
         Data.Yarr.Flow
-        Data.Yarr.Fold
+        Data.Yarr.Work
         Data.Yarr.Shape
         Data.Yarr.Repr.Foreign
         Data.Yarr.Repr.Boxed
         Data.Yarr.Repr.Delayed
         Data.Yarr.Repr.Separate
         Data.Yarr.Convolution
+        Data.Yarr.IO.List
         Data.Yarr.Utils.FixedVector
         Data.Yarr.Utils.Fork
         Data.Yarr.Utils.Parallel
@@ -85,8 +86,6 @@
         Debug.Yarr
 
     other-modules:
-        Data.Yarr.Utils.Storable
-
         -- re-exported in Utils.FixedVector
         Data.Yarr.Utils.FixedVector.Arity
         Data.Yarr.Utils.FixedVector.VecTuple
@@ -98,3 +97,9 @@
         Data.Yarr.Convolution.Repr
         Data.Yarr.Convolution.Eval
         Data.Yarr.Convolution.StaticStencils
+
+        -- re-exported in Data.Yarr.Shape
+        Data.Yarr.WorkTypes
+
+        Data.Yarr.Work.Internal
+        Data.Yarr.Utils.Storable