diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,14 @@
+# 0.5.4
+
+* Addition of `unsafeTransformStencil`
+* Add `zip4`, `unzip4`, `zipWith4`  and `izipWith4`
+* Make `Resize` a superclass of `Source`
+* Addition of `outerSlices`, `innerSlices`, `withinSlices` and `withinSlicesM`
+* Addition of `stackSlicesM`, `stackOuterSlicesM` and `stackInnerSlicesM`
+* Addition of `computeP`
+* Fix perfomrmance issue of folding functions applied to arrays with `Seq` computation
+  strategy.
+
 # 0.5.3
 
 * Fix `tanA` and `tanhA`. [#96](https://github.com/lehins/massiv/pull/96)
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.5.3.2
+version:             0.5.4.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
diff --git a/src/Data/Massiv/Array.hs b/src/Data/Massiv/Array.hs
--- a/src/Data/Massiv/Array.hs
+++ b/src/Data/Massiv/Array.hs
@@ -85,6 +85,7 @@
   , setComp
   , compute
   , computeS
+  , computeP
   , computeIO
   , computePrimM
   , computeAs
@@ -167,7 +168,6 @@
 import Prelude as P hiding (all, and, any, enumFromTo, foldl, foldr, mapM,
                             mapM_, maximum, minimum, or, product, replicate, splitAt,
                             sum, zip)
-
 
 {- $folding
 
diff --git a/src/Data/Massiv/Array/Delayed/Pull.hs b/src/Data/Massiv/Array/Delayed/Pull.hs
--- a/src/Data/Massiv/Array/Delayed/Pull.hs
+++ b/src/Data/Massiv/Array/Delayed/Pull.hs
@@ -110,6 +110,8 @@
 instance Functor (Array D ix) where
   fmap f (DArray c sz g) = DArray c sz (f . g)
   {-# INLINE fmap #-}
+  (<$) e (DArray c sz _) = DArray c sz (const e)
+  {-# INLINE (<$) #-}
 
 
 instance Index ix => Applicative (Array D ix) where
@@ -148,8 +150,7 @@
   {-# INLINE size #-}
   getComp = dComp
   {-# INLINE getComp #-}
-  loadArrayM !scheduler !arr =
-    splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
+  loadArrayM !scheduler !arr = splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
   {-# INLINE loadArrayM #-}
 
 instance Index ix => StrideLoad D ix e
diff --git a/src/Data/Massiv/Array/Delayed/Push.hs b/src/Data/Massiv/Array/Delayed/Push.hs
--- a/src/Data/Massiv/Array/Delayed/Push.hs
+++ b/src/Data/Massiv/Array/Delayed/Push.hs
@@ -56,9 +56,12 @@
 instance Index ix => Construct DL ix e where
   setComp c arr = arr {dlComp = c}
   {-# INLINE setComp #-}
-  makeArrayLinear comp sz f =
-    DLArray comp sz Nothing $ \scheduler startAt dlWrite ->
-      splitLinearlyWithStartAtM_ scheduler startAt (totalElem sz) (pure . f) dlWrite
+  makeArrayLinear comp sz f = DLArray comp sz Nothing load
+    where
+      load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+      load scheduler startAt dlWrite =
+        splitLinearlyWithStartAtM_ scheduler startAt (totalElem sz) (pure . f) dlWrite
+      {-# INLINE load #-}
   {-# INLINE makeArrayLinear #-}
 
 instance Index ix => Resize DL ix where
@@ -171,21 +174,25 @@
 --
 -- @since 0.3.1
 makeLoadArrayS ::
-     Index ix =>
-     Sz ix
+     forall ix e. Index ix
+  => Sz ix
   -- ^ Size of the resulting array
   -> e
   -- ^ Default value to use for all cells that might have been ommitted by the writing function
-  -> (forall m. Monad m => (ix -> e -> m Bool) -> m ())
+  -> (forall m. Monad m =>
+                  (ix -> e -> m Bool) -> m ())
   -- ^ Writing function that described which elements to write into the target array.
   -> Array DL ix e
-makeLoadArrayS sz defVal writer =
-  DLArray Seq sz (Just defVal) $ \_scheduler !startAt uWrite ->
-    let safeWrite !ix !e
-          | isSafeIndex sz ix = uWrite (startAt + toLinearIndex sz ix) e >> pure True
-          | otherwise = pure False
-        {-# INLINE safeWrite #-}
-     in writer safeWrite
+makeLoadArrayS sz defVal writer = DLArray Seq sz (Just defVal) load
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load _scheduler !startAt uWrite =
+      let safeWrite !ix !e
+            | isSafeIndex sz ix = uWrite (startAt + toLinearIndex sz ix) e >> pure True
+            | otherwise = pure False
+          {-# INLINE safeWrite #-}
+       in writer safeWrite
+    {-# INLINE load #-}
 {-# INLINE makeLoadArrayS #-}
 
 -- | Specify how an array should be loaded into memory. Unlike `makeLoadArrayS`, loading
@@ -194,7 +201,7 @@
 --
 -- @since 0.4.0
 makeLoadArray ::
-     Index ix
+     forall ix e. Index ix
   => Comp
   -- ^ Computation strategy to use. Directly affects the scheduler that gets created for
   -- the loading function.
@@ -208,13 +215,16 @@
   -- accepts a scheduler, that can be used for parallelization, as well as a safe element
   -- writing function.
   -> Array DL ix e
-makeLoadArray comp sz defVal writer =
-  DLArray comp sz (Just defVal) $ \scheduler !startAt uWrite ->
-    let safeWrite !ix !e
-          | isSafeIndex sz ix = uWrite (startAt + toLinearIndex sz ix) e >> pure True
-          | otherwise = pure False
-        {-# INLINE safeWrite #-}
-     in writer scheduler safeWrite
+makeLoadArray comp sz defVal writer = DLArray comp sz (Just defVal) load
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler !startAt uWrite =
+      let safeWrite !ix !e
+            | isSafeIndex sz ix = uWrite (startAt + toLinearIndex sz ix) e >> pure True
+            | otherwise = pure False
+          {-# INLINE safeWrite #-}
+       in writer scheduler safeWrite
+    {-# INLINE load #-}
 {-# INLINE makeLoadArray #-}
 
 -- | Specify how an array can be loaded/computed through creation of a `DL` array. Unlike
@@ -250,36 +260,51 @@
 --
 -- @since 0.5.2
 unsafeMakeLoadArrayAdjusted ::
+     forall ix e.
      Comp
   -> Sz ix
   -> Maybe e
-  -> (forall m. Monad m => Scheduler m () -> (Int -> e -> m ()) -> m ())
+  -> (forall m. Monad m =>
+                  Scheduler m () -> (Int -> e -> m ()) -> m ())
   -> Array DL ix e
-unsafeMakeLoadArrayAdjusted comp sz mDefVal writer =
-  DLArray comp sz mDefVal $ \scheduler !startAt uWrite ->
-    writer scheduler (\i -> uWrite (startAt + i))
+unsafeMakeLoadArrayAdjusted comp sz mDefVal writer = DLArray comp sz mDefVal load
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler !startAt uWrite = writer scheduler (\i -> uWrite (startAt + i))
+    {-# INLINE load #-}
 {-# INLINE unsafeMakeLoadArrayAdjusted #-}
 
 -- | Convert any `Load`able array into `DL` representation.
 --
 -- @since 0.3.0
-toLoadArray :: Load r ix e => Array r ix e -> Array DL ix e
-toLoadArray arr =
-  DLArray (getComp arr) (size arr) Nothing $ \scheduler startAt dlWrite ->
-    loadArrayM scheduler arr (dlWrite . (+ startAt))
+toLoadArray ::
+     forall r ix e. Load r ix e
+  => Array r ix e
+  -> Array DL ix e
+toLoadArray arr = DLArray (getComp arr) (size arr) (defaultElement arr) load
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler !startAt dlWrite = loadArrayM scheduler arr (dlWrite . (+ startAt))
+    {-# INLINE load #-}
 {-# INLINE[1] toLoadArray #-}
 {-# RULES "toLoadArray/id" toLoadArray = id #-}
 
 -- | Convert an array that can be loaded with stride into `DL` representation.
 --
 -- @since 0.3.0
-fromStrideLoad
-  :: StrideLoad r ix e => Stride ix -> Array r ix e -> Array DL ix e
+fromStrideLoad ::
+     forall r ix e. StrideLoad r ix e
+  => Stride ix
+  -> Array r ix e
+  -> Array DL ix e
 fromStrideLoad stride arr =
-  DLArray (getComp arr) newsz Nothing $ \scheduler startAt dlWrite ->
-    loadArrayWithStrideM scheduler stride newsz arr (\ !i -> dlWrite (i + startAt))
+  DLArray (getComp arr) newsz Nothing load
   where
     newsz = strideSize stride (size arr)
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler !startAt dlWrite =
+      loadArrayWithStrideM scheduler stride newsz arr (\ !i -> dlWrite (i + startAt))
+    {-# INLINE load #-}
 {-# INLINE fromStrideLoad #-}
 
 instance Index ix => Load DL ix e where
@@ -293,10 +318,12 @@
   {-# INLINE defaultElement #-}
 
 instance Functor (Array DL ix) where
-  fmap f arr =
-    arr
-      { dlLoad =
-          \scheduler startAt uWrite -> dlLoad arr scheduler startAt (\ !i e -> uWrite i (f e))
-      , dlDefault = f <$> dlDefault arr
-      }
+  fmap f arr = arr {dlLoad = loadFunctor arr f, dlDefault = f <$> dlDefault arr}
   {-# INLINE fmap #-}
+  (<$) e arr = arr {dlLoad = \_ _ _ -> pure (), dlDefault = Just e}
+  {-# INLINE (<$) #-}
+
+loadFunctor ::
+     Monad m => Array DL ix t1 -> (t1 -> t2) -> Scheduler m () -> Int -> (Int -> t2 -> m ()) -> m ()
+loadFunctor arr f scheduler startAt uWrite = dlLoad arr scheduler startAt (\ !i e -> uWrite i (f e))
+{-# INLINE loadFunctor #-}
diff --git a/src/Data/Massiv/Array/Manifest/Internal.hs b/src/Data/Massiv/Array/Manifest/Internal.hs
--- a/src/Data/Massiv/Array/Manifest/Internal.hs
+++ b/src/Data/Massiv/Array/Manifest/Internal.hs
@@ -23,6 +23,7 @@
   , toManifest
   , compute
   , computeS
+  , computeP
   , computeIO
   , computePrimM
   , computeAs
@@ -221,6 +222,19 @@
 computeS :: forall r ix e r' . (Mutable r ix e, Load r' ix e) => Array r' ix e -> Array r ix e
 computeS !arr = runST $ computePrimM arr
 {-# INLINE computeS #-}
+
+
+-- | Compute array in parallel using all cores disregarding predefined computation
+-- strategy. Computation stategy of the resulting array will match the source, despite
+-- that it is diregarded.
+--
+-- @since 0.5.4
+computeP ::
+     forall r ix e r'. (Mutable r ix e, Construct r' ix e, Load r' ix e)
+  => Array r' ix e
+  -> Array r ix e
+computeP arr = setComp (getComp arr) $ compute (setComp Par arr)
+{-# INLINE computeP #-}
 
 -- | Very similar to `compute`, but computes an array inside the `IO` monad. Despite being
 -- deterministic and referentially transparent, because this is an `IO` action it
diff --git a/src/Data/Massiv/Array/Mutable.hs b/src/Data/Massiv/Array/Mutable.hs
--- a/src/Data/Massiv/Array/Mutable.hs
+++ b/src/Data/Massiv/Array/Mutable.hs
@@ -418,7 +418,6 @@
 --   [ 10, 12 ]
 --
 -- @since 0.3.0
---
 createArrayS_ ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the newly created array
@@ -431,7 +430,6 @@
 -- | Just like `createArray_`, but together with `Array` it returns the result of the filling action.
 --
 -- @since 0.3.0
---
 createArrayS ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the newly created array
@@ -448,7 +446,6 @@
 -- | Just like `createArrayS_`, but restricted to `ST`.
 --
 -- @since 0.3.0
---
 createArrayST_ ::
      forall r ix e a. Mutable r ix e
   => Sz ix
@@ -461,7 +458,6 @@
 -- | Just like `createArrayS`, but restricted to `ST`.
 --
 -- @since 0.2.6
---
 createArrayST ::
      forall r ix e a. Mutable r ix e
   => Sz ix
@@ -477,8 +473,6 @@
 -- but is a stateful action, becuase it is restricted to `PrimMonad` thus allows for sharing
 -- the state between computation of each element.
 --
--- @since 0.2.6
---
 -- ====__Examples__
 --
 -- >>> import Data.Massiv.Array
@@ -496,6 +490,7 @@
 -- >>> readIORef ref
 -- 15
 --
+-- @since 0.2.6
 generateArrayS ::
      forall r ix e m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Resulting size of the array
@@ -603,7 +598,6 @@
 --   [ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34 ]
 --
 -- @since 0.3.0
---
 unfoldrPrimM_ ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the desired array
@@ -616,7 +610,6 @@
 -- | Same as `unfoldrPrimM_` but do the unfolding with index aware function.
 --
 -- @since 0.3.0
---
 iunfoldrPrimM_ ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the desired array
@@ -686,7 +679,6 @@
 --   [ 34, 21, 13, 8, 5, 3, 2, 1, 1, 0 ]
 --
 -- @since 0.3.0
---
 unfoldlPrimM_ ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the desired array
@@ -699,7 +691,6 @@
 -- | Same as `unfoldlPrimM_` but do the unfolding with index aware function.
 --
 -- @since 0.3.0
---
 iunfoldlPrimM_ ::
      forall r ix e a m. (Mutable r ix e, PrimMonad m)
   => Sz ix -- ^ Size of the desired array
diff --git a/src/Data/Massiv/Array/Ops/Construct.hs b/src/Data/Massiv/Array/Ops/Construct.hs
--- a/src/Data/Massiv/Array/Ops/Construct.hs
+++ b/src/Data/Massiv/Array/Ops/Construct.hs
@@ -225,21 +225,21 @@
 -- | Right unfold of a delayed load array with index aware function
 --
 -- @since 0.3.0
-iunfoldrS_
-  :: Construct DL ix e => Sz ix -> (a -> ix -> (e, a)) -> a -> Array DL ix e
-iunfoldrS_ sz f acc0 =
-  DLArray
-    { dlComp = Seq
-    , dlSize = sz
-    , dlDefault = Nothing
-    , dlLoad =
-        \_ startAt dlWrite ->
-          void $
-          loopM startAt (< (totalElem sz + startAt)) (+ 1) acc0 $ \ !i !acc -> do
-            let (e, acc') = f acc $ fromLinearIndex sz (i - startAt)
-            dlWrite i e
-            pure acc'
-    }
+iunfoldrS_ ::
+     forall ix e a. Construct DL ix e
+  => Sz ix
+  -> (a -> ix -> (e, a))
+  -> a
+  -> Array DL ix e
+iunfoldrS_ sz f acc0 = DLArray {dlComp = Seq, dlSize = sz, dlDefault = Nothing, dlLoad = load}
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load _ startAt dlWrite =
+      void $
+      loopM startAt (< (totalElem sz + startAt)) (+ 1) acc0 $ \ !i !acc ->
+        let (e, acc') = f acc $ fromLinearIndex sz (i - startAt)
+         in acc' <$ dlWrite i e
+    {-# INLINE load #-}
 {-# INLINE iunfoldrS_ #-}
 
 
@@ -255,20 +255,21 @@
 -- | Unfold sequentially from the right with an index aware function.
 --
 -- @since 0.3.0
-iunfoldlS_
-  :: Construct DL ix e => Sz ix -> (ix -> a -> (a, e)) -> a -> Array DL ix e
-iunfoldlS_ sz f acc0 =
-  DLArray
-    { dlComp = Seq
-    , dlSize = sz
-    , dlDefault = Nothing
-    , dlLoad =
-        \ _ startAt dlWrite ->
-          void $ loopDeepM startAt (< (totalElem sz + startAt)) (+ 1) acc0 $ \ !i !acc -> do
-            let (acc', e) = f (fromLinearIndex sz (i - startAt)) acc
-            dlWrite i e
-            pure acc'
-    }
+iunfoldlS_ ::
+     forall ix e a. Construct DL ix e
+  => Sz ix
+  -> (ix -> a -> (a, e))
+  -> a
+  -> Array DL ix e
+iunfoldlS_ sz f acc0 = DLArray {dlComp = Seq, dlSize = sz, dlDefault = Nothing, dlLoad = load}
+  where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load _ startAt dlWrite =
+      void $
+      loopDeepM startAt (< (totalElem sz + startAt)) (+ 1) acc0 $ \ !i !acc ->
+        let (acc', e) = f (fromLinearIndex sz (i - startAt)) acc
+         in acc' <$ dlWrite i e
+    {-# INLINE load #-}
 {-# INLINE iunfoldlS_ #-}
 
 
@@ -312,28 +313,28 @@
   -> Comp -- ^ Computation strategy.
   -> Sz ix -- ^ Resulting size of the array.
   -> Array DL ix e
-randomArray gen splitGen nextRandom comp sz =
-  unsafeMakeLoadArray comp sz Nothing $ \scheduler startAt writeAt ->
-    splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-      let slackStartAt = slackStart + startAt
-          writeRandom k genII = do
-            let (e, genII') = nextRandom genII
-            writeAt k e
-            pure genII'
-      genForSlack <-
-        loopM startAt (< slackStartAt) (+ chunkLength) gen $ \start genI -> do
-          let (genI0, genI1) =
-                if numWorkers scheduler == 1
-                  then (genI, genI)
-                  else splitGen genI
-          scheduleWork_ scheduler $
-            void $ loopM start (< (start + chunkLength)) (+ 1) genI0 writeRandom
-          pure genI1
-      when (slackStartAt < totalLength + startAt) $
-        scheduleWork_ scheduler $
-        void $ loopM slackStartAt (< totalLength + startAt) (+ 1) genForSlack writeRandom
+randomArray gen splitGen nextRandom comp sz = unsafeMakeLoadArray comp sz Nothing load
   where
     !totalLength = totalElem sz
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler startAt writeAt =
+      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+        let slackStartAt = slackStart + startAt
+            writeRandom k genII =
+              let (e, genII') = nextRandom genII
+               in genII' <$ writeAt k e
+        genForSlack <-
+          loopM startAt (< slackStartAt) (+ chunkLength) gen $ \start genI -> do
+            let (genI0, genI1) =
+                  if numWorkers scheduler == 1
+                    then (genI, genI)
+                    else splitGen genI
+            scheduleWork_ scheduler $
+              void $ loopM start (< (start + chunkLength)) (+ 1) genI0 writeRandom
+            pure genI1
+        when (slackStartAt < totalLength + startAt) $
+          scheduleWork_ scheduler $
+          void $ loopM slackStartAt (< totalLength + startAt) (+ 1) genForSlack writeRandom
 {-# INLINE randomArray #-}
 
 -- | Similar to `randomArray` but performs generation sequentially, which means it doesn't
@@ -421,7 +422,7 @@
 
 infix 4 ..., ..:
 
--- | Handy synonym for `rangeInclusive` `Seq`
+-- | Handy synonym for @`rangeInclusive` `Seq`@. Similar to @..@ for list.
 --
 -- >>> Ix1 4 ... 10
 -- Array D Seq (Sz1 7)
@@ -432,7 +433,7 @@
 (...) = rangeInclusive Seq
 {-# INLINE (...) #-}
 
--- | Handy synonym for `range` `Seq`
+-- | Handy synonym for @`range` `Seq`@
 --
 -- >>> Ix1 4 ..: 10
 -- Array D Seq (Sz1 6)
diff --git a/src/Data/Massiv/Array/Ops/Fold/Internal.hs b/src/Data/Massiv/Array/Ops/Fold/Internal.hs
--- a/src/Data/Massiv/Array/Ops/Fold/Internal.hs
+++ b/src/Data/Massiv/Array/Ops/Fold/Internal.hs
@@ -311,21 +311,23 @@
   -> b -- ^ Accumulator for chunks folding
   -> Array r ix e
   -> m b
-ifoldlIO f !initAcc g !tAcc !arr = do
-  let !sz = size arr
-      !totalLength = totalElem sz
-  results <-
-    withScheduler (getComp arr) $ \scheduler ->
-      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
-          scheduleWork scheduler $
-          iterLinearM sz start (start + chunkLength) 1 (<) initAcc $ \ !i ix !acc ->
-            f acc ix (unsafeLinearIndex arr i)
-        when (slackStart < totalLength) $
-          scheduleWork scheduler $
-          iterLinearM sz slackStart totalLength 1 (<) initAcc $ \ !i ix !acc ->
-            f acc ix (unsafeLinearIndex arr i)
-  F.foldlM g tAcc results
+ifoldlIO f !initAcc g !tAcc !arr
+  | getComp arr == Seq = ifoldlM f initAcc arr >>= g tAcc
+  | otherwise = do
+      let !sz = size arr
+          !totalLength = totalElem sz
+      results <-
+        withScheduler (getComp arr) $ \scheduler ->
+          splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+            loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+              scheduleWork scheduler $
+              iterLinearM sz start (start + chunkLength) 1 (<) initAcc $ \ !i ix !acc ->
+                f acc ix (unsafeLinearIndex arr i)
+            when (slackStart < totalLength) $
+              scheduleWork scheduler $
+              iterLinearM sz slackStart totalLength 1 (<) initAcc $ \ !i ix !acc ->
+                f acc ix (unsafeLinearIndex arr i)
+      F.foldlM g tAcc results
 {-# INLINE ifoldlIO #-}
 
 -- -- | Split an array into linear row-major vector chunks and apply an action to each of
@@ -361,19 +363,21 @@
 -- @since 0.1.0
 ifoldrIO :: (MonadUnliftIO m, Source r ix e) =>
            (ix -> e -> a -> m a) -> a -> (a -> b -> m b) -> b -> Array r ix e -> m b
-ifoldrIO f !initAcc g !tAcc !arr = do
-  let !sz = size arr
-      !totalLength = totalElem sz
-  results <-
-    withScheduler (getComp arr) $ \ scheduler ->
-      splitLinearly (numWorkers scheduler) totalLength $ \ chunkLength slackStart -> do
-        when (slackStart < totalLength) $
-          scheduleWork scheduler $
-          iterLinearM sz (totalLength - 1) slackStart (-1) (>=) initAcc $ \ !i ix !acc ->
-            f ix (unsafeLinearIndex arr i) acc
-        loopM_ slackStart (> 0) (subtract chunkLength) $ \ !start ->
-          scheduleWork scheduler $
-            iterLinearM sz (start - 1) (start - chunkLength) (-1) (>=) initAcc $ \ !i ix !acc ->
+ifoldrIO f !initAcc g !tAcc !arr
+  | getComp arr == Seq = ifoldrM f initAcc arr >>= (`g` tAcc)
+  | otherwise = do
+    let !sz = size arr
+        !totalLength = totalElem sz
+    results <-
+      withScheduler (getComp arr) $ \ scheduler ->
+        splitLinearly (numWorkers scheduler) totalLength $ \ chunkLength slackStart -> do
+          when (slackStart < totalLength) $
+            scheduleWork scheduler $
+            iterLinearM sz (totalLength - 1) slackStart (-1) (>=) initAcc $ \ !i ix !acc ->
               f ix (unsafeLinearIndex arr i) acc
-  F.foldlM (flip g) tAcc results
+          loopM_ slackStart (> 0) (subtract chunkLength) $ \ !start ->
+            scheduleWork scheduler $
+              iterLinearM sz (start - 1) (start - chunkLength) (-1) (>=) initAcc $ \ !i ix !acc ->
+                f ix (unsafeLinearIndex arr i) acc
+    F.foldlM (flip g) tAcc results
 {-# INLINE ifoldrIO #-}
diff --git a/src/Data/Massiv/Array/Ops/Map.hs b/src/Data/Massiv/Array/Ops/Map.hs
--- a/src/Data/Massiv/Array/Ops/Map.hs
+++ b/src/Data/Massiv/Array/Ops/Map.hs
@@ -52,12 +52,16 @@
   -- ** Zipping
   , zip
   , zip3
+  , zip4
   , unzip
   , unzip3
+  , unzip4
   , zipWith
   , zipWith3
+  , zipWith4
   , izipWith
   , izipWith3
+  , izipWith4
   , liftArray2
   -- *** Applicative
   , zipWithA
@@ -102,6 +106,19 @@
 zip3 = zipWith3 (,,)
 {-# INLINE zip3 #-}
 
+-- | Zip four arrays
+--
+-- @since 0.5.4
+zip4 ::
+     (Source r1 ix e1, Source r2 ix e2, Source r3 ix e3, Source r4 ix e4)
+  => Array r1 ix e1
+  -> Array r2 ix e2
+  -> Array r3 ix e3
+  -> Array r4 ix e4
+  -> Array D ix (e1, e2, e3, e4)
+zip4 = zipWith4 (,,,)
+{-# INLINE zip4 #-}
+
 -- | Unzip two arrays
 unzip :: Source r ix (e1, e2) => Array r ix (e1, e2) -> (Array D ix e1, Array D ix e2)
 unzip arr = (map fst arr, map snd arr)
@@ -113,6 +130,18 @@
 unzip3 arr = (map (\ (e, _, _) -> e) arr, map (\ (_, e, _) -> e) arr, map (\ (_, _, e) -> e) arr)
 {-# INLINE unzip3 #-}
 
+-- | Unzip four arrays
+--
+-- @since 0.5.4
+unzip4 :: Source r ix (e1, e2, e3, e4)
+       => Array r ix (e1, e2, e3, e4) -> (Array D ix e1, Array D ix e2, Array D ix e3, Array D ix e4)
+unzip4 arr =
+  ( map (\(e, _, _, _) -> e) arr
+  , map (\(_, e, _, _) -> e) arr
+  , map (\(_, _, e, _) -> e) arr
+  , map (\(_, _, _, e) -> e) arr)
+{-# INLINE unzip4 #-}
+
 --------------------------------------------------------------------------------
 -- zipWith ---------------------------------------------------------------------
 --------------------------------------------------------------------------------
@@ -161,6 +190,48 @@
           (coerce (size arr3)))) $ \ !ix ->
     f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix)
 {-# INLINE izipWith3 #-}
+
+
+
+-- | Just like `zipWith`, except zip four arrays with a function.
+--
+-- @since 0.5.4
+zipWith4 ::
+     (Source r1 ix e1, Source r2 ix e2, Source r3 ix e3, Source r4 ix e4)
+  => (e1 -> e2 -> e3 -> e4 -> e)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Array r3 ix e3
+  -> Array r4 ix e4
+  -> Array D ix e
+zipWith4 f = izipWith4 (\ _ e1 e2 e3 e4 -> f e1 e2 e3 e4)
+{-# INLINE zipWith4 #-}
+
+
+-- | Just like `zipWith4`, except with an index aware function.
+--
+-- @since 0.5.4
+izipWith4
+  :: (Source r1 ix e1, Source r2 ix e2, Source r3 ix e3, Source r4 ix e4)
+  => (ix -> e1 -> e2 -> e3 -> e4 -> e)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Array r3 ix e3
+  -> Array r4 ix e4
+  -> Array D ix e
+izipWith4 f arr1 arr2 arr3 arr4 =
+  DArray
+    (getComp arr1 <> getComp arr2 <> getComp arr3 <> getComp arr4)
+    (SafeSz
+       (liftIndex2
+          min
+          (liftIndex2
+             min
+             (liftIndex2 min (coerce (size arr1)) (coerce (size arr2)))
+             (coerce (size arr3)))
+          (coerce (size arr4)))) $ \ !ix ->
+    f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix) (unsafeIndex arr4 ix)
+{-# INLINE izipWith4 #-}
 
 
 -- | Similar to `zipWith`, except does it sequentially and using the `Applicative`. Note that
diff --git a/src/Data/Massiv/Array/Ops/Slice.hs b/src/Data/Massiv/Array/Ops/Slice.hs
--- a/src/Data/Massiv/Array/Ops/Slice.hs
+++ b/src/Data/Massiv/Array/Ops/Slice.hs
@@ -1,4 +1,6 @@
 {-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Slice
 -- Copyright   : (c) Alexey Kuleshevich 2018-2019
@@ -21,9 +23,15 @@
   , (<!>)
   , (<!?>)
   , (<??>)
+  -- ** Many slices
+  , outerSlices
+  , innerSlices
+  , withinSlices
+  , withinSlicesM
   ) where
 
 import Control.Monad (unless)
+import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Core.Common
 
 
@@ -147,12 +155,20 @@
 (<!?>) !arr (dim, i) = do
   (m, szl) <- pullOutSzM (size arr) dim
   unless (isSafeIndex m i) $ throwM $ IndexOutOfBoundsException m i
-  start <- setDimM zeroIndex dim i
   cutSz <- insertSzM szl dim oneSz
-  unsafeSlice arr start cutSz dim
+  internalInnerSlice dim cutSz arr i
 {-# INLINE (<!?>) #-}
 
 
+internalInnerSlice ::
+     (MonadThrow m, Slice r ix e) => Dim -> Sz ix -> Array r ix e -> Int -> m (Elt r ix e)
+internalInnerSlice dim cutSz arr i = do
+  start <- setDimM zeroIndex dim i
+  unsafeSlice arr start cutSz dim
+{-# INLINE internalInnerSlice #-}
+
+
+
 -- prop> arr !> i == arr <!> (dimensions (size arr), i)
 -- prop> arr <! i == arr <!> (1,i)
 --
@@ -174,3 +190,120 @@
 (<??>) :: (MonadThrow m, Slice r ix e) => m (Array r ix e) -> (Dim, Int) -> m (Elt r ix e)
 (<??>) !marr !ix = marr >>= (<!?> ix)
 {-# INLINE (<??>) #-}
+
+-- | Create a delayed array of outer slices.
+--
+-- ====__Examples__
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> A.mapM_ print $ outerSlices (0 ..: (3 :. 2))
+-- Array D Seq (Sz1 2)
+--   [ 0 :. 0, 0 :. 1 ]
+-- Array D Seq (Sz1 2)
+--   [ 1 :. 0, 1 :. 1 ]
+-- Array D Seq (Sz1 2)
+--   [ 2 :. 0, 2 :. 1 ]
+--
+-- @since 0.5.4
+outerSlices :: OuterSlice r ix e => Array r ix e -> Array D Ix1 (Elt r ix e)
+outerSlices arr = makeArray Seq k (unsafeOuterSlice arr)
+  where
+    (k, _) = unconsSz $ size arr
+-- TODO: move setComp to Load
+-- outerSlices arr = makeArray (getComp arr) k (unsafeOuterSlice arr')
+--   where
+--     arr' = setComp Seq arr
+--     (k, _) = unconsSz $ size arr
+{-# INLINE outerSlices #-}
+
+
+-- | Create a delayed array of inner slices.
+--
+-- ====__Examples__
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> A.mapM_ print $ innerSlices (0 ..: (3 :. 2))
+-- Array D Seq (Sz1 3)
+--   [ 0 :. 0, 1 :. 0, 2 :. 0 ]
+-- Array D Seq (Sz1 3)
+--   [ 0 :. 1, 1 :. 1, 2 :. 1 ]
+--
+-- @since 0.5.4
+innerSlices :: InnerSlice r ix e => Array r ix e -> Array D Ix1 (Elt r ix e)
+innerSlices arr = makeArray Seq k (unsafeInnerSlice arr sz)
+  where
+    sz@(_, k) = unsnocSz $ size arr
+-- TODO: move setComp to Load
+-- innerSlices arr = makeArray (getComp arr) k (unsafeInnerSlice arr' sz)
+--   where
+--     arr' = setComp Seq arr
+--     sz@(_, k) = unsnocSz $ size arr
+{-# INLINE innerSlices #-}
+
+-- | Create a delayed array of slices from within. Checks dimension at compile time.
+--
+-- ====__Examples__
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> arr = fromIx3 <$> (0 ..: (4 :> 3 :. 2))
+-- >>> print arr
+-- Array D Seq (Sz (4 :> 3 :. 2))
+--   [ [ [ (0,0,0), (0,0,1) ]
+--     , [ (0,1,0), (0,1,1) ]
+--     , [ (0,2,0), (0,2,1) ]
+--     ]
+--   , [ [ (1,0,0), (1,0,1) ]
+--     , [ (1,1,0), (1,1,1) ]
+--     , [ (1,2,0), (1,2,1) ]
+--     ]
+--   , [ [ (2,0,0), (2,0,1) ]
+--     , [ (2,1,0), (2,1,1) ]
+--     , [ (2,2,0), (2,2,1) ]
+--     ]
+--   , [ [ (3,0,0), (3,0,1) ]
+--     , [ (3,1,0), (3,1,1) ]
+--     , [ (3,2,0), (3,2,1) ]
+--     ]
+--   ]
+-- >>> A.mapM_ print $ withinSlices Dim2 arr
+-- Array D Seq (Sz (4 :. 2))
+--   [ [ (0,0,0), (0,0,1) ]
+--   , [ (1,0,0), (1,0,1) ]
+--   , [ (2,0,0), (2,0,1) ]
+--   , [ (3,0,0), (3,0,1) ]
+--   ]
+-- Array D Seq (Sz (4 :. 2))
+--   [ [ (0,1,0), (0,1,1) ]
+--   , [ (1,1,0), (1,1,1) ]
+--   , [ (2,1,0), (2,1,1) ]
+--   , [ (3,1,0), (3,1,1) ]
+--   ]
+-- Array D Seq (Sz (4 :. 2))
+--   [ [ (0,2,0), (0,2,1) ]
+--   , [ (1,2,0), (1,2,1) ]
+--   , [ (2,2,0), (2,2,1) ]
+--   , [ (3,2,0), (3,2,1) ]
+--   ]
+--
+-- @since 0.5.4
+withinSlices ::
+     (IsIndexDimension ix n, Slice r ix e)
+  => Dimension n
+  -> Array r ix e
+  -> Array D Ix1 (Elt r ix e)
+withinSlices dim = either throwImpossible id . withinSlicesM (fromDimension dim)
+{-# INLINE withinSlices #-}
+
+
+-- | Create a delayed array of slices from within. Same as `withinSlices`, but throws an
+-- error on invalid dimension.
+--
+-- /__Throws Exceptions__/: `IndexDimensionException`
+--
+-- @since 0.5.4
+withinSlicesM :: (MonadThrow m, Slice r ix e) => Dim -> Array r ix e -> m (Array D Ix1 (Elt r ix e))
+withinSlicesM dim arr = do
+  (k, szl) <- pullOutSzM (size arr) dim
+  cutSz <- insertSzM szl dim oneSz
+  pure $ makeArray Seq k (either throwImpossible id . internalInnerSlice dim cutSz arr)
+{-# INLINE withinSlicesM #-}
diff --git a/src/Data/Massiv/Array/Ops/Transform.hs b/src/Data/Massiv/Array/Ops/Transform.hs
--- a/src/Data/Massiv/Array/Ops/Transform.hs
+++ b/src/Data/Massiv/Array/Ops/Transform.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE ExplicitForAll #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# OPTIONS_GHC -fno-warn-redundant-constraints #-}
 -- |
@@ -43,6 +43,9 @@
   , concatOuterM
   , concatM
   , concat'
+  , stackSlicesM
+  , stackOuterSlicesM
+  , stackInnerSlicesM
   , splitAtM
   , splitAt'
   , splitExtractM
@@ -62,7 +65,7 @@
 import Control.Scheduler (traverse_)
 import Control.Monad as M (foldM_, unless, forM_)
 import Data.Bifunctor (bimap)
-import Data.Foldable as F (foldl', foldrM, toList)
+import Data.Foldable as F (foldl', foldrM, toList, length)
 import qualified Data.List as L (uncons)
 import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Array.Delayed.Push
@@ -393,7 +396,7 @@
 -- an allowed exception of the dimension they are being appended along, otherwise `Nothing` is
 -- returned.
 --
--- ===__Examples__
+-- ====__Examples__
 --
 -- Append two 2D arrays along both dimensions. Note that they do agree on inner dimensions.
 --
@@ -425,32 +428,34 @@
 -- *** Exception: SizeMismatchException: (Sz (2 :. 3)) vs (Sz (2 :. 4))
 --
 -- @since 0.3.0
-appendM :: (MonadThrow m, Source r1 ix e, Source r2 ix e) =>
-          Dim -> Array r1 ix e -> Array r2 ix e -> m (Array DL ix e)
+appendM ::
+     forall r1 r2 ix e m. (MonadThrow m, Source r1 ix e, Source r2 ix e)
+  => Dim
+  -> Array r1 ix e
+  -> Array r2 ix e
+  -> m (Array DL ix e)
 appendM n !arr1 !arr2 = do
   let !sz1 = size arr1
       !sz2 = size arr2
   (k1, szl1) <- pullOutSzM sz1 n
   (k2, szl2) <- pullOutSzM sz2 n
   unless (szl1 == szl2) $ throwM $ SizeMismatchException sz1 sz2
-  let k1' = unSz k1
+  let !k1' = unSz k1
   newSz <- insertSzM szl1 n (SafeSz (k1' + unSz k2))
+  let load :: Monad n => Scheduler n () -> Int -> (Int -> e -> n ()) -> n ()
+      load scheduler !startAt dlWrite = do
+        scheduleWork scheduler $
+          iterM_ zeroIndex (unSz sz1) (pureIndex 1) (<) $ \ix ->
+            dlWrite (startAt + toLinearIndex newSz ix) (unsafeIndex arr1 ix)
+        scheduleWork scheduler $
+          iterM_ zeroIndex (unSz sz2) (pureIndex 1) (<) $ \ix ->
+            let i = getDim' ix n
+                ix' = setDim' ix n (i + k1')
+             in dlWrite (startAt + toLinearIndex newSz ix') (unsafeIndex arr2 ix)
+      {-# INLINE load #-}
   return $
     DLArray
-      { dlComp = getComp arr1 <> getComp arr2
-      , dlSize = newSz
-      , dlDefault = Nothing
-      , dlLoad =
-          \scheduler startAt dlWrite -> do
-            scheduleWork scheduler $
-              iterM_ zeroIndex (unSz sz1) (pureIndex 1) (<) $ \ix ->
-                dlWrite (startAt + toLinearIndex newSz ix) (unsafeIndex arr1 ix)
-            scheduleWork scheduler $
-              iterM_ zeroIndex (unSz sz2) (pureIndex 1) (<) $ \ix ->
-                let i = getDim' ix n
-                    ix' = setDim' ix n (i + k1')
-                 in dlWrite (startAt + toLinearIndex newSz ix') (unsafeIndex arr2 ix)
-      }
+      {dlComp = getComp arr1 <> getComp arr2, dlSize = newSz, dlDefault = Nothing, dlLoad = load}
 {-# INLINE appendM #-}
 
 
@@ -476,13 +481,16 @@
 --
 -- @since 0.3.0
 concatM ::
-     (MonadThrow m, Foldable f, Source r ix e) => Dim -> f (Array r ix e) -> m (Array DL ix e)
+     forall r ix e f m. (MonadThrow m, Foldable f, Source r ix e)
+  => Dim
+  -> f (Array r ix e)
+  -> m (Array DL ix e)
 concatM n !arrsF =
   case L.uncons (F.toList arrsF) of
     Nothing -> pure empty
     Just (a, arrs) -> do
       let sz = unSz (size a)
-          szs = P.map (unSz . size) arrs
+          szs = unSz . size <$> arrs
       (k, szl) <- pullOutDimM sz n
       -- / remove the dimension out of all sizes along which concatenation will happen
       (ks, szls) <-
@@ -493,25 +501,192 @@
         (dropWhile ((== szl) . snd) $ P.zip szs szls)
       let kTotal = SafeSz $ F.foldl' (+) k ks
       newSz <- insertSzM (SafeSz szl) n kTotal
+      let load :: Monad n => Scheduler n () -> Int -> (Int -> e -> n ()) -> n ()
+          load scheduler startAt dlWrite =
+            let arrayLoader !kAcc (kCur, arr) = do
+                  scheduleWork scheduler $
+                    iforM_ arr $ \ix e ->
+                      let i = getDim' ix n
+                          ix' = setDim' ix n (i + kAcc)
+                       in dlWrite (startAt + toLinearIndex newSz ix') e
+                  pure (kAcc + kCur)
+             in M.foldM_ arrayLoader 0 $ (k, a) : P.zip ks arrs
+          {-# INLINE load #-}
       return $
-        DLArray
-          { dlComp = mconcat $ P.map getComp arrs
-          , dlSize = newSz
-          , dlDefault = Nothing
-          , dlLoad =
-              \scheduler startAt dlWrite ->
-                let arrayLoader !kAcc (kCur, arr) = do
-                      scheduleWork scheduler $
-                        iterM_ zeroIndex (unSz (size arr)) (pureIndex 1) (<) $ \ix ->
-                          let i = getDim' ix n
-                              ix' = setDim' ix n (i + kAcc)
-                           in dlWrite (startAt + toLinearIndex newSz ix') (unsafeIndex arr ix)
-                      pure (kAcc + kCur)
-                 in M.foldM_ arrayLoader 0 $ (k, a) : P.zip ks arrs
-          }
+        DLArray {dlComp = foldMap getComp arrsF, dlSize = newSz, dlDefault = Nothing, dlLoad = load}
 {-# INLINE concatM #-}
 
 
+-- | Stack slices on top of each other along the specified dimension.
+--
+-- /__Exceptions__/: `IndexDimensionException`, `SizeMismatchException`
+--
+-- ====__Examples__
+--
+-- Here are the three different ways to stack up two 2D Matrix pages into a 3D array.
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> x = compute (iterateN 3 succ 0) :: Matrix P Int
+-- >>> y = compute (iterateN 3 succ 9) :: Matrix P Int
+-- >>> x
+-- Array P Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+-- >>> y
+-- Array P Seq (Sz (3 :. 3))
+--   [ [ 10, 11, 12 ]
+--   , [ 13, 14, 15 ]
+--   , [ 16, 17, 18 ]
+--   ]
+-- >>> stackSlicesM 1 [x, y] :: IO (Array DL Ix3 Int)
+-- Array DL Seq (Sz (3 :> 3 :. 2))
+--   [ [ [ 1, 10 ]
+--     , [ 2, 11 ]
+--     , [ 3, 12 ]
+--     ]
+--   , [ [ 4, 13 ]
+--     , [ 5, 14 ]
+--     , [ 6, 15 ]
+--     ]
+--   , [ [ 7, 16 ]
+--     , [ 8, 17 ]
+--     , [ 9, 18 ]
+--     ]
+--   ]
+-- >>> stackSlicesM 2 [x, y] :: IO (Array DL Ix3 Int)
+-- Array DL Seq (Sz (3 :> 2 :. 3))
+--   [ [ [ 1, 2, 3 ]
+--     , [ 10, 11, 12 ]
+--     ]
+--   , [ [ 4, 5, 6 ]
+--     , [ 13, 14, 15 ]
+--     ]
+--   , [ [ 7, 8, 9 ]
+--     , [ 16, 17, 18 ]
+--     ]
+--   ]
+-- >>> stackSlicesM 3 [x, y] :: IO (Array DL Ix3 Int)
+-- Array DL Seq (Sz (2 :> 3 :. 3))
+--   [ [ [ 1, 2, 3 ]
+--     , [ 4, 5, 6 ]
+--     , [ 7, 8, 9 ]
+--     ]
+--   , [ [ 10, 11, 12 ]
+--     , [ 13, 14, 15 ]
+--     , [ 16, 17, 18 ]
+--     ]
+--   ]
+--
+-- @since 0.5.4
+stackSlicesM ::
+     forall r ix e f m. (Foldable f, MonadThrow m, Source r (Lower ix) e, Index ix)
+  => Dim
+  -> f (Array r (Lower ix) e)
+  -> m (Array DL ix e)
+stackSlicesM dim !arrsF = do
+  case L.uncons (F.toList arrsF) of
+    Nothing -> pure empty
+    Just (a, arrs) -> do
+      let sz = size a
+          len = SafeSz (F.length arrsF)
+      -- / make sure all arrays have the same size
+      M.forM_ arrsF $ \arr ->
+        let sz' = size arr
+         in unless (sz == sz') $ throwM (SizeMismatchException sz sz')
+      newSz <- insertSzM sz dim len
+      let load :: Monad n => Scheduler n () -> Int -> (Int -> e -> n ()) -> n ()
+          load scheduler startAt dlWrite =
+            let loadIndex k ix = dlWrite (toLinearIndex newSz (insertDim' ix dim k) + startAt)
+                arrayLoader !k arr = (k + 1) <$ scheduleWork scheduler (imapM_ (loadIndex k) arr)
+             in M.foldM_ arrayLoader 0 arrsF
+          {-# INLINE load #-}
+      return $
+        DLArray {dlComp = foldMap getComp arrs, dlSize = newSz, dlDefault = Nothing, dlLoad = load}
+{-# INLINE stackSlicesM #-}
+
+-- | Specialized `stackOuterM` to handling stacking from the outside. It is the inverse of
+-- `Data.Massiv.Array.outerSlices`.
+--
+-- /__Exceptions__/: `SizeMismatchException`
+--
+-- ====__Examples__
+--
+-- In this example we stack vectors as row of a matrix from top to bottom:
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> x = compute (iterateN 3 succ 0) :: Matrix P Int
+-- >>> x
+-- Array P Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+-- >>> rows = outerSlices x
+-- >>> A.mapM_ print rows
+-- Array M Seq (Sz1 3)
+--   [ 1, 2, 3 ]
+-- Array M Seq (Sz1 3)
+--   [ 4, 5, 6 ]
+-- Array M Seq (Sz1 3)
+--   [ 7, 8, 9 ]
+-- >>> stackOuterSlicesM rows :: IO (Matrix DL Int)
+-- Array DL Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+--
+-- @since 0.5.4
+stackOuterSlicesM ::
+     forall r ix e f m. (Foldable f, MonadThrow m, Source r (Lower ix) e, Index ix)
+  => f (Array r (Lower ix) e)
+  -> m (Array DL ix e)
+stackOuterSlicesM = stackSlicesM (dimensions (Proxy :: Proxy ix))
+{-# INLINE stackOuterSlicesM #-}
+
+-- | Specialized `stackOuterM` to handling stacking from the inside. It is the inverse of
+-- `Data.Massiv.Array.outerSlices`.
+--
+-- /__Exceptions__/: `SizeMismatchException`
+--
+-- ====__Examples__
+--
+-- In this example we stack vectors as columns of a matrix from left to right:
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> x = compute (iterateN 3 succ 0) :: Matrix P Int
+-- >>> x
+-- Array P Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+-- >>> columns = innerSlices x
+-- >>> A.mapM_ print columns
+-- Array M Seq (Sz1 3)
+--   [ 1, 4, 7 ]
+-- Array M Seq (Sz1 3)
+--   [ 2, 5, 8 ]
+-- Array M Seq (Sz1 3)
+--   [ 3, 6, 9 ]
+-- >>> stackInnerSlicesM columns :: IO (Matrix DL Int)
+-- Array DL Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+--
+-- @since 0.5.4
+stackInnerSlicesM ::
+     forall r ix e f m. (Foldable f, MonadThrow m, Source r (Lower ix) e, Index ix)
+  => f (Array r (Lower ix) e)
+  -> m (Array DL ix e)
+stackInnerSlicesM = stackSlicesM 1
+{-# INLINE stackInnerSlicesM #-}
+
+
 -- | /O(1)/ - Split an array into two at an index along a specified dimension.
 --
 -- /Related/: `splitAt'`, `splitExtractM`, `Data.Massiv.Vector.sliceAt'`, `Data.Massiv.Vector.sliceAtM`
@@ -669,51 +844,84 @@
 -- | Discard elements from the source array according to the stride.
 --
 -- @since 0.3.0
-downsample :: Source r ix e => Stride ix -> Array r ix e -> Array DL ix e
+downsample ::
+     forall r ix e. Source r ix e
+  => Stride ix
+  -> Array r ix e
+  -> Array DL ix e
 downsample stride arr =
-  DLArray
-    { dlComp = getComp arr
-    , dlSize = resultSize
-    , dlDefault = defaultElement arr
-    , dlLoad =
-        \scheduler startAt dlWrite ->
-          splitLinearlyWithStartAtM_
-            scheduler
-            startAt
-            (totalElem resultSize)
-            (pure . unsafeLinearWriteWithStride)
-            dlWrite
-    }
+  DLArray {dlComp = getComp arr, dlSize = resultSize, dlDefault = defaultElement arr, dlLoad = load}
   where
     resultSize = strideSize stride (size arr)
     strideIx = unStride stride
     unsafeLinearWriteWithStride =
       unsafeIndex arr . liftIndex2 (*) strideIx . fromLinearIndex resultSize
     {-# INLINE unsafeLinearWriteWithStride #-}
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler startAt dlWrite =
+      splitLinearlyWithStartAtM_
+        scheduler
+        startAt
+        (totalElem resultSize)
+        (pure . unsafeLinearWriteWithStride)
+        dlWrite
+    {-# INLINE load #-}
 {-# INLINE downsample #-}
 
 
--- | Insert the same element into a `Load`able array according to the stride.
+-- | Insert the same element into a `Load`able array according to the supplied stride.
 --
+-- ====__Examples__
+--
+-- >>> import Data.Massiv.Array as A
+-- >>> arr = iterateN (Sz2 3 2) succ (0 :: Int)
+-- >>> arr
+-- Array DL Seq (Sz (3 :. 2))
+--   [ [ 1, 2 ]
+--   , [ 3, 4 ]
+--   , [ 5, 6 ]
+--   ]
+-- >>> upsample 0 (Stride (2 :. 3)) arr
+-- Array DL Seq (Sz (6 :. 6))
+--   [ [ 1, 0, 0, 2, 0, 0 ]
+--   , [ 0, 0, 0, 0, 0, 0 ]
+--   , [ 3, 0, 0, 4, 0, 0 ]
+--   , [ 0, 0, 0, 0, 0, 0 ]
+--   , [ 5, 0, 0, 6, 0, 0 ]
+--   , [ 0, 0, 0, 0, 0, 0 ]
+--   ]
+-- >>> upsample 9 (Stride (1 :. 2)) arr
+-- Array DL Seq (Sz (3 :. 4))
+--   [ [ 1, 9, 2, 9 ]
+--   , [ 3, 9, 4, 9 ]
+--   , [ 5, 9, 6, 9 ]
+--   ]
+--
 -- @since 0.3.0
-upsample
-  :: Load r ix e => e -> Stride ix -> Array r ix e -> Array DL ix e
+upsample ::
+     forall r ix e. Load r ix e
+  => e -- ^ Element to use for filling the newly added cells
+  -> Stride ix -- ^ Fill cells according to this stride
+  -> Array r ix e -- ^ Array that will have cells added to
+  -> Array DL ix e
 upsample !fillWith safeStride arr =
   DLArray
     { dlComp = getComp arr
     , dlSize = newsz
     , dlDefault = Just fillWith
-    , dlLoad =
-        \scheduler startAt dlWrite -> do
-          M.forM_ (defaultElement arr) $ \prevFillWith ->
-            loopM_
-              startAt
-              (< totalElem sz)
-              (+ 1)
-              (\i -> dlWrite (adjustLinearStride (i + startAt)) prevFillWith)
-          loadArrayM scheduler arr (\i -> dlWrite (adjustLinearStride (i + startAt)))
+    , dlLoad = load
     }
   where
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler !startAt dlWrite = do
+      M.forM_ (defaultElement arr) $ \prevFillWith ->
+        loopM_
+          startAt
+          (< totalElem sz)
+          (+ 1)
+          (\i -> dlWrite (adjustLinearStride (i + startAt)) prevFillWith)
+      loadArrayM scheduler arr (\i -> dlWrite (adjustLinearStride (i + startAt)))
+    {-# INLINE load #-}
     adjustLinearStride = toLinearIndex newsz . timesStride . fromLinearIndex sz
     {-# INLINE adjustLinearStride #-}
     timesStride !ix = liftIndex2 (*) stride ix
@@ -819,21 +1027,23 @@
 --
 -- @since 0.3.1
 zoomWithGrid ::
-     Source r ix e
+     forall r ix e. Source r ix e
   => e -- ^ Value to use for the grid
   -> Stride ix -- ^ Scaling factor
   -> Array r ix e -- ^ Source array
   -> Array DL ix e
-zoomWithGrid gridVal (Stride zoomFactor) arr =
-  unsafeMakeLoadArray Seq newSz (Just gridVal) $ \scheduler _ writeElement ->
-    iforSchedulerM_ scheduler arr $ \ !ix !e -> do
-      let !kix = liftIndex2 (*) ix kx
-      mapM_ (\ !ix' -> writeElement (toLinearIndex newSz ix') e) $
-        range Seq (liftIndex (+1) kix) (liftIndex2 (+) kix kx)
+zoomWithGrid gridVal (Stride zoomFactor) arr = unsafeMakeLoadArray Seq newSz (Just gridVal) load
   where
-    !kx = liftIndex (+1) zoomFactor
+    !kx = liftIndex (+ 1) zoomFactor
     !lastNewIx = liftIndex2 (*) kx $ unSz (size arr)
-    !newSz = Sz (liftIndex (+1) lastNewIx)
+    !newSz = Sz (liftIndex (+ 1) lastNewIx)
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler _ writeElement =
+      iforSchedulerM_ scheduler arr $ \ !ix !e ->
+        let !kix = liftIndex2 (*) ix kx
+         in mapM_ (\ !ix' -> writeElement (toLinearIndex newSz ix') e) $
+            range Seq (liftIndex (+ 1) kix) (liftIndex2 (+) kix kx)
+    {-# INLINE load #-}
 {-# INLINE zoomWithGrid #-}
 
 -- | Increaze the size of the array accoridng to the stride multiplier while replicating
@@ -871,17 +1081,19 @@
 --
 -- @since 0.4.4
 zoom ::
-     Source r ix e
+     forall r ix e. Source r ix e
   => Stride ix -- ^ Scaling factor
   -> Array r ix e -- ^ Source array
   -> Array DL ix e
-zoom (Stride zoomFactor) arr =
-  unsafeMakeLoadArray Seq newSz Nothing $ \scheduler _ writeElement ->
-    iforSchedulerM_ scheduler arr $ \ !ix !e -> do
-      let !kix = liftIndex2 (*) ix zoomFactor
-      mapM_ (\ !ix' -> writeElement (toLinearIndex newSz ix') e) $
-        range Seq kix (liftIndex2 (+) kix zoomFactor)
+zoom (Stride zoomFactor) arr = unsafeMakeLoadArray Seq newSz Nothing load
   where
     !lastNewIx = liftIndex2 (*) zoomFactor $ unSz (size arr)
     !newSz = Sz lastNewIx
+    load :: Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ()
+    load scheduler _ writeElement =
+      iforSchedulerM_ scheduler arr $ \ !ix !e ->
+        let !kix = liftIndex2 (*) ix zoomFactor
+         in mapM_ (\ !ix' -> writeElement (toLinearIndex newSz ix') e) $
+            range Seq kix (liftIndex2 (+) kix zoomFactor)
+    {-# INLINE load #-}
 {-# INLINE zoom #-}
diff --git a/src/Data/Massiv/Array/Stencil.hs b/src/Data/Massiv/Array/Stencil.hs
--- a/src/Data/Massiv/Array/Stencil.hs
+++ b/src/Data/Massiv/Array/Stencil.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE TypeFamilies #-}
@@ -274,6 +275,7 @@
 idStencil = makeUnsafeStencil oneSz zeroIndex $ \ _ get -> get zeroIndex
 {-# INLINE idStencil #-}
 
+
 -- | Stencil that does a left fold in a row-major order. Regardless of the supplied size
 -- resulting stencil will be centered at zero, although by using `Padding` it is possible
 -- to overcome this limitation.
@@ -334,6 +336,9 @@
 {-# INLINE foldrStencil #-}
 
 
+-- | Create a stencil that will fold all elements in the region monoidally.
+--
+-- @since 0.4.3
 foldStencil :: (Monoid e, Index ix) => Sz ix -> Stencil ix e e
 foldStencil = foldlStencil mappend mempty
 {-# INLINE foldStencil #-}
diff --git a/src/Data/Massiv/Array/Stencil/Internal.hs b/src/Data/Massiv/Array/Stencil/Internal.hs
--- a/src/Data/Massiv/Array/Stencil/Internal.hs
+++ b/src/Data/Massiv/Array/Stencil/Internal.hs
@@ -35,6 +35,7 @@
   , stencilFunc   :: (ix -> Value e) -> ix -> Value a
   }
 
+
 instance Index ix => NFData (Stencil ix e a) where
   rnf (Stencil sz ix f) = sz `deepseq` ix `deepseq` f `seq` ()
 
diff --git a/src/Data/Massiv/Array/Stencil/Unsafe.hs b/src/Data/Massiv/Array/Stencil/Unsafe.hs
--- a/src/Data/Massiv/Array/Stencil/Unsafe.hs
+++ b/src/Data/Massiv/Array/Stencil/Unsafe.hs
@@ -2,6 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RecordWildCards #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Unsafe
 -- Copyright   : (c) Alexey Kuleshevich 2018-2019
@@ -13,6 +14,7 @@
 module Data.Massiv.Array.Stencil.Unsafe
   ( -- * Stencil
     makeUnsafeStencil
+  , unsafeTransformStencil
   , unsafeMapStencil
   -- ** Deprecated
   , mapStencilUnsafe
@@ -123,3 +125,80 @@
       Value $ inline $ relStencil ix (unValue . getVal . liftIndex2 (+) ix)
     {-# INLINE stencil #-}
 {-# INLINE makeUnsafeStencil #-}
+
+
+-- | Perform an arbitrary transformation of a stencil. This stencil modifier can be used for
+-- example to turn a vector stencil into a matrix stencil implement, or transpose a matrix
+-- stencil. It is really easy to get this wrong, so be extremely careful.
+--
+-- ====__Examples__
+--
+-- Convert a 1D stencil into a row or column 2D stencil:
+--
+-- >>> import Data.Massiv.Array
+-- >>> import Data.Massiv.Array.Unsafe
+-- >>> let arr = compute $ iterateN 3 succ 0 :: Array P Ix2 Int
+-- >>> arr
+-- Array P Seq (Sz (3 :. 3))
+--   [ [ 1, 2, 3 ]
+--   , [ 4, 5, 6 ]
+--   , [ 7, 8, 9 ]
+--   ]
+-- >>> let rowStencil = unsafeTransformStencil (\(Sz n) -> Sz (1 :. n)) (0 :.) $ \ f getVal (i :. j) -> f (getVal . (i :.)) j
+-- >>> applyStencil noPadding (rowStencil (sumStencil (Sz1 3))) arr
+-- Array DW Seq (Sz (3 :. 1))
+--   [ [ 6 ]
+--   , [ 15 ]
+--   , [ 24 ]
+--   ]
+-- >>> let columnStencil = unsafeTransformStencil (\(Sz n) -> Sz (n :. 1)) (:. 0) $ \ f getVal (i :. j) -> f (getVal . (:. j)) i
+-- >>> applyStencil noPadding (columnStencil (sumStencil (Sz1 3))) arr
+-- Array DW Seq (Sz (1 :. 3))
+--   [ [ 12, 15, 18 ]
+--   ]
+--
+-- @since 0.5.4
+unsafeTransformStencil ::
+     (Sz ix' -> Sz ix)
+  -- ^ Forward modifier for the size
+  -> (ix' -> ix)
+  -- ^ Forward index modifier
+  -> (((ix' -> Value e) -> ix' -> Value a) -> (ix -> Value e) -> ix -> Value a)
+  -- ^ Inverse stencil function modifier
+  -> Stencil ix' e a
+  -- ^ Original stencil.
+  -> Stencil ix e a
+unsafeTransformStencil transformSize transformIndex transformFunc Stencil {..} =
+  Stencil
+    { stencilSize = transformSize stencilSize
+    , stencilCenter = transformIndex stencilCenter
+    , stencilFunc = transformFunc stencilFunc
+    }
+{-# INLINE unsafeTransformStencil #-}
+
+
+
+{-
+
+Invalid stencil transformer function.
+
+TODO: figure out if there is a safe way to do stencil index trnasformation.
+
+
+transformStencil ::
+     (Default e, Index ix)
+  => (Sz ix' -> Sz ix)
+  -- ^ Forward modifier for the size
+  -> (ix' -> ix)
+  -- ^ Forward index modifier
+  -> (ix -> ix')
+  -- ^ Inverse index modifier
+  -> Stencil ix' e a
+  -- ^ Original stencil.
+  -> Stencil ix e a
+transformStencil transformSize transformIndex transformIndex' stencil =
+  validateStencil def $! unsafeTransformStencil transformSize transformIndex transformIndex' stencil
+{-# INLINE transformStencil #-}
+
+
+-}
diff --git a/src/Data/Massiv/Core/Common.hs b/src/Data/Massiv/Core/Common.hs
--- a/src/Data/Massiv/Core/Common.hs
+++ b/src/Data/Massiv/Core/Common.hs
@@ -235,7 +235,7 @@
 
 
 -- | Arrays that can be used as source to practically any manipulation function.
-class Load r ix e => Source r ix e where
+class (Resize r ix, Load r ix e) => Source r ix e where
   {-# MINIMAL (unsafeIndex|unsafeLinearIndex), unsafeLinearSlice #-}
 
   -- | Lookup element in the array. No bounds check is performed and access of
@@ -392,7 +392,7 @@
 -- | Manifest arrays are backed by actual memory and values are looked up versus
 -- computed as it is with delayed arrays. Because of this fact indexing functions
 -- @(`!`)@, @(`!?`)@, etc. are constrained to manifest arrays only.
-class (Load r ix e, Source r ix e) => Manifest r ix e where
+class Source r ix e => Manifest r ix e where
 
   unsafeLinearIndexM :: Array r ix e -> Int -> e
 
diff --git a/src/Data/Massiv/Vector.hs b/src/Data/Massiv/Vector.hs
--- a/src/Data/Massiv/Vector.hs
+++ b/src/Data/Massiv/Vector.hs
@@ -959,12 +959,9 @@
 cons :: Load r Ix1 e => e -> Vector r e -> Vector DL e
 cons e v =
   let dv = toLoadArray v
-   in dv
-        { dlSize = SafeSz (1 + unSz (dlSize dv))
-        , dlLoad =
-            \scheduler startAt uWrite ->
-              uWrite startAt e >> dlLoad dv scheduler (startAt + 1) uWrite
-        }
+      load scheduler startAt uWrite = uWrite startAt e >> dlLoad dv scheduler (startAt + 1) uWrite
+      {-# INLINE load #-}
+   in dv {dlSize = SafeSz (1 + unSz (dlSize dv)), dlLoad = load}
 {-# INLINE cons #-}
 
 -- | /O(1)/ - Add an element to the vector from the right side
@@ -974,12 +971,9 @@
 snoc v e =
   let dv = toLoadArray v
       !k = unSz (size dv)
-   in dv
-        { dlSize = SafeSz (1 + k)
-        , dlLoad =
-            \scheduler startAt uWrite ->
-              dlLoad dv scheduler startAt uWrite >> uWrite (k + startAt) e
-        }
+      load scheduler startAt uWrite = dlLoad dv scheduler startAt uWrite >> uWrite (k + startAt) e
+      {-# INLINE load #-}
+   in dv {dlSize = SafeSz (1 + k), dlLoad = load}
 {-# INLINE snoc #-}
 
 
@@ -1433,7 +1427,7 @@
 
 -- | Similar to `filterM`, but map with an index aware function.
 --
--- Corresponds to: @filterM (uncurry f) . imap (,)@
+-- Corresponds to: @`filterM` (uncurry f) . `simap` (,)@
 --
 -- ==== __Examples__
 --