diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,45 @@
+# 1.0.2
+
+* Addition of `Iterator` type class and related fucntions:
+  * Addition of `RowMajor`, `RowMajorLinear` and `RowMajorUnbalanced` iterators.
+  * Switch parallel left fold to new iterator
+* Improvements to functions that do the looping:
+  * Addition of `loopNextA_` and `loopNextM`
+  * Deprecate `loopM_` in favor of `loopA_`
+  * Addition of `loopA` and `loopF` for applicative iterators
+  * Addition of `iloopA_` and `iloopM`
+  * Addition of `nextMaybeF`
+  * Addition of `iterLinearST_`, `iterLinearAccST_` and `iterLinearAccST`
+  * Addition of an optimized `scheduleMassivWork` for internal use
+  * Addition of a new workhorse: `splitWorkWithFactorST`
+  * Addition of a new workhorse: `splitWorkWithFactorST`
+* Changes to `Index` class:
+  * Deprecate `iterM_` in favor of `iterA_`
+  * Adititon of sequential iterators: `iterTargetRowMajorA_`,
+    `iterTargetRowMajorAccM` needed for `RowMajor` iterator
+  * Addition of parallelizable iterators: `iterRowMajorST`,
+    `iterTargetRowMajorAccST`, `iterTargetRowMajorAccST_` needed
+    for `RowMajor` iterator
+  * Addition of `iterF` for using with applicative iterators.
+  * Addition of `stepNextMF` for streaming iteration of multi-dimensional
+    arrays.
+* Addition of `repr`.
+* Addition of `quicksortAs`, `quicksortAsBy`, `quicksortAsByM`
+* Fix backwards compatibility with ghc-8.0
+* Get rid of dependency on `genvalidity`: too many compatibility issues for
+  little gain
+* Introduce `PrefIndex` and `unsafePrefIndex`: a preference when indexing into
+  multidimensional `Source` arrays. Adopt it where possible for left and right
+  folds, traversals, but not zipping
+* Improve multi-dimensional indices for streams. Improve `steps` and `isteps`
+* Get rid of build warnings for ghc-9.4
+  * Make `Monoid` and `Monad` instances forward compatible
+* Compatibility with `vector-0.13`:
+  * Fix `Unbox` instance for `IxN`
+  * Fix safety of boxed vector conversion: `toBoxedVector`/`fromBoxedVector`
+* Re-export Manifest realetd functions from `Data.Massiv.Array.Manifest`
+  as a migration strategy for the move in the next minor version bump.
+
 # 1.0.1
 
 * Relax constraint on `computeInto` by removing requirement for `Size`
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright Alexey Kuleshevich (c) 2017-2021
+Copyright Alexey Kuleshevich (c) 2017-2022
 
 All rights reserved.
 
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -7,14 +7,13 @@
 
 ## Status
 
-| Language | Github Actions | Coveralls |Gitter.im |
-|:--------:|:--------------:|:---------:|:--------:|
+| Language | Github Actions | Coveralls | Gitter.im |
+|:--------:|:--------------:|:---------:|:---------:|
 | ![GitHub top language](https://img.shields.io/github/languages/top/lehins/massiv.svg) | [![Build Status](https://github.com/lehins/massiv/workflows/massiv-CI/badge.svg)](https://github.com/lehins/massiv/actions) | [![Coverage Status](https://coveralls.io/repos/github/lehins/massiv/badge.svg?branch=master)](https://coveralls.io/github/lehins/massiv?branch=master) | [![Join the chat at https://gitter.im/haskell-massiv/Lobby](https://badges.gitter.im/haskell-massiv/Lobby.svg)](https://gitter.im/haskell-massiv/Lobby?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
 
 |      Package       | Hackage | Nightly | LTS |
 |:-------------------|:-------:|:-------:|:---:|
-|  [`massiv`](https://github.com/lehins/massiv/tree/master/massiv)|                                       [![Hackage](https://img.shields.io/hackage/v/massiv.svg)](https://hackage.haskell.org/package/massiv)|                                                                                                        [![Nightly](https://www.stackage.org/package/massiv/badge/nightly)](https://www.stackage.org/nightly/package/massiv)|                                                                                         [![Nightly](https://www.stackage.org/package/massiv/badge/lts)](https://www.stackage.org/lts/package/massiv)|
-|  [`massiv-io`](https://github.com/lehins/massiv-io)|                                [![Hackage](https://img.shields.io/hackage/v/massiv-io.svg)](https://hackage.haskell.org/package/massiv-io)|                                                                                                  [![Nightly](https://www.stackage.org/package/massiv-io/badge/nightly)](https://www.stackage.org/nightly/package/massiv-io)|                                                                                   [![Nightly](https://www.stackage.org/package/massiv-io/badge/lts)](https://www.stackage.org/lts/package/massiv-io)|
+|  [`massiv`](https://github.com/lehins/massiv/tree/master/massiv)|                                       [![Hackage](https://img.shields.io/hackage/v/massiv.svg)](https://hackage.haskell.org/package/massiv)|                                                                                                        [![Nightly](https://www.stackage.org/package/massiv/badge/nightly)](https://www.stackage.org/nightly/package/massiv)| ](https://www.stackage.org/package/massiv-io/badge/lts)](https://www.stackage.org/lts/package/massiv-io)|
 |  [`massiv-test`](https://github.com/lehins/massiv/tree/master/massiv-test)|                            [![Hackage](https://img.shields.io/hackage/v/massiv-test.svg)](https://hackage.haskell.org/package/massiv-test)|                                                                                              [![Nightly](https://www.stackage.org/package/massiv-test/badge/nightly)](https://www.stackage.org/nightly/package/massiv-test)|                                                                               [![Nightly](https://www.stackage.org/package/massiv-test/badge/lts)](https://www.stackage.org/lts/package/massiv-test)|
 |  [`haskell-scheduler`](https://github.com/lehins/haskell-scheduler)|                                   [![Hackage](https://img.shields.io/hackage/v/scheduler.svg)](https://hackage.haskell.org/package/scheduler)|                                                                                          [![Nightly](https://www.stackage.org/package/scheduler/badge/nightly)](https://www.stackage.org/nightly/package/scheduler)|                                                                   [![Nightly](https://www.stackage.org/package/scheduler/badge/lts)](https://www.stackage.org/lts/package/scheduler)|
 
@@ -666,5 +665,6 @@
 
 # Further resources on learning `massiv`:
 
+* [2021 - Haskell eXchange - Multi-dimensional Arrays that Do Not Exist](#2021---haskell-exchange---multi-dimensional-arrays-that-do-not-exist)
 * [2019 - Monadic Party - Haskell arrays with Massiv](https://github.com/lehins/talks#2019---monadic-party---haskell-arrays-with-massiv)
 * [2018 - Monadic Warsaw #14 - Haskell arrays that are easy and fast](https://github.com/lehins/talks#2018---monadic-warsaw-14---haskell-arrays-that-are-easy-and-fast)
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             1.0.1.1
+version:             1.0.2.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
@@ -7,7 +7,7 @@
 license-file:        LICENSE
 author:              Alexey Kuleshevich
 maintainer:          alexey@kuleshevi.ch
-copyright:           2018-2021 Alexey Kuleshevich
+copyright:           2018-2022 Alexey Kuleshevich
 category:            Array, Data, Data Structures, Parallelism
 build-type:          Simple
 extra-source-files:  README.md
@@ -19,11 +19,12 @@
                     , GHC == 8.6.5
                     , GHC == 8.8.4
                     , GHC == 8.10.7
-                    , GHC == 9.0.1
+                    , GHC == 9.0.2
+                    , GHC == 9.2.4
 
 flag unsafe-checks
   description: Enable all the bounds checks for unsafe functions at the cost of
-               significant performance penalty
+               performance penalty
   default: False
   manual: True
 
@@ -75,17 +76,18 @@
                      , Data.Massiv.Core.Index.Ix
                      , Data.Massiv.Core.Index.Stride
                      , Data.Massiv.Core.Index.Tuple
-                     , Data.Massiv.Core.Iterator
+                     , Data.Massiv.Core.Index.Iterator
+                     , Data.Massiv.Core.Loop
                      , Data.Massiv.Vector.Unsafe
   build-depends:       base >= 4.9 && < 5
                      , bytestring
                      , deepseq
                      , exceptions
-                     , scheduler >= 2.0.0 && < 3.0.0
+                     , scheduler >= 2.0.0 && < 3
                      , primitive >= 0.7.1.0
                      , random >= 1.2.0
                      , unliftio-core
-                     , vector
+                     , vector >= 0.12
 
   include-dirs: include
   install-includes: massiv.h
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
@@ -1,7 +1,7 @@
 {-# OPTIONS_GHC -fno-warn-duplicate-exports #-}
 -- |
 -- Module      : Data.Massiv.Array
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Delayed.hs b/src/Data/Massiv/Array/Delayed.hs
--- a/src/Data/Massiv/Array/Delayed.hs
+++ b/src/Data/Massiv/Array/Delayed.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      : Data.Massiv.Array.Delayed
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Delayed/Interleaved.hs b/src/Data/Massiv/Array/Delayed/Interleaved.hs
--- a/src/Data/Massiv/Array/Delayed/Interleaved.hs
+++ b/src/Data/Massiv/Array/Delayed/Interleaved.hs
@@ -7,7 +7,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Interleaved
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -27,6 +27,8 @@
 
 -- | Delayed array that will be loaded in an interleaved fashion during parallel
 -- computation.
+--
+-- /Warning/ - Will be deprecated in the next major version update.
 data DI = DI
 
 newtype instance Array DI ix e = DIArray
@@ -42,6 +44,7 @@
   {-# INLINE setComp #-}
   getComp = dComp . diArray
   {-# INLINE getComp #-}
+  repr = DI
 
 instance Index ix => Shape DI ix where
   maxLinearSize = Just . SafeSz . elemsCount
@@ -58,24 +61,24 @@
 instance Index ix => Load DI ix e where
   makeArray c sz = DIArray . makeArray c sz
   {-# INLINE makeArray #-}
-  iterArrayLinearST_ scheduler (DIArray (DArray _ sz f)) uWrite =
-    loopM_ 0 (< numWorkers scheduler) (+ 1) $ \ !start ->
-      scheduleWork scheduler $
-      iterLinearM_ sz start (totalElem sz) (numWorkers scheduler) (<) $ \ !k -> uWrite k . f
+  iterArrayLinearST_ scheduler (DIArray darr@(DArray _ sz _)) uWrite =
+    loopA_ 0 (< numWorkers scheduler) (+ 1) $ \ !start ->
+      scheduleWork_ scheduler $
+        iterLinearM_ sz start (totalElem sz) (numWorkers scheduler) (<) $ \ !k ->
+          uWrite k . unsafeIndex darr
   {-# INLINE iterArrayLinearST_ #-}
 
 instance Index ix => StrideLoad DI ix e where
-  iterArrayLinearWithStrideST_ scheduler stride resultSize arr uWrite =
-    let strideIx = unStride stride
-        DIArray (DArray _ _ f) = arr
-    in loopM_ 0 (< numWorkers scheduler) (+ 1) $ \ !start ->
-          scheduleWork scheduler $
-          iterLinearM_ resultSize start (totalElem resultSize) (numWorkers scheduler) (<) $
-            \ !i ix -> uWrite i (f (liftIndex2 (*) strideIx ix))
+  iterArrayLinearWithStrideST_ scheduler stride resultSize (DIArray arr) uWrite =
+    loopA_ 0 (< numWorkers scheduler) (+ 1) $ \ !start ->
+      scheduleWork_ scheduler $
+        iterLinearM_ resultSize start (totalElem resultSize) (numWorkers scheduler) (<) $
+          \ !i ix -> uWrite i (unsafeIndex arr (liftIndex2 (*) (unStride stride) ix))
   {-# INLINE iterArrayLinearWithStrideST_ #-}
 
 -- | Convert a source array into an array that, when computed, will have its elemets evaluated out
 -- of order (interleaved amongst cores), hence making unbalanced computation better parallelizable.
+--
 toInterleaved :: (Index ix, Source r e) => Array r ix e -> Array DI ix e
 toInterleaved = DIArray . delay
 {-# INLINE toInterleaved #-}
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
@@ -3,17 +3,20 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UndecidableInstances #-}
+
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Pull
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
 -- Portability : non-portable
---
 module Data.Massiv.Array.Delayed.Pull
   ( D(..)
   , Array(..)
@@ -26,28 +29,32 @@
   , unsafeExtract
   , unsafeSlice
   , unsafeInnerSlice
+  , zipWithInternal
   ) where
 
-import           Control.Applicative
+import Control.Applicative
 import qualified Data.Foldable as F
-import           Data.Massiv.Array.Ops.Fold.Internal as A
-import           Data.Massiv.Core.Common as A
-import           Data.Massiv.Core.List (L, showArrayList, showsArrayPrec)
-import           Data.Massiv.Core.Operations
-import           Data.Massiv.Vector.Stream as S (steps)
-import           GHC.Base (build)
-import           Prelude hiding (zipWith)
+import Data.Massiv.Array.Ops.Fold.Internal as A
+import Data.Massiv.Core.Common as A
+import Data.Massiv.Core.List (L, showArrayList, showsArrayPrec)
+import Data.Massiv.Core.Operations
+import qualified Data.Massiv.Vector.Stream as S
+import GHC.Base (build)
+import Prelude hiding (zipWith)
 
 #include "massiv.h"
 
-
 -- | Delayed representation.
-data D = D deriving Show
+data D =
+  D
+  deriving (Show)
 
 
-data instance Array D ix e = DArray { dComp :: !Comp
-                                    , dSize :: !(Sz ix)
-                                    , dIndex :: ix -> e }
+data instance  Array D ix e =
+  DArray { dComp :: !Comp
+         , dSize :: !(Sz ix)
+         , dPrefIndex :: !(PrefIndex ix e)
+         }
 
 instance (Ragged L ix e, Show e) => Show (Array D ix e) where
   showsPrec = showsArrayPrec id
@@ -61,38 +68,51 @@
   size = dSize
   {-# INLINE size #-}
   unsafeResize !sz !arr =
-    DArray (dComp arr) sz $ \ !ix ->
-      unsafeIndex arr (fromLinearIndex (size arr) (toLinearIndex sz ix))
+    makeArrayLinear (dComp arr) sz (unsafeIndex arr . fromLinearIndex (size arr))
   {-# INLINE unsafeResize #-}
 
 instance Strategy D where
-  setComp c arr = arr { dComp = c }
+  setComp c arr = arr {dComp = c}
   {-# INLINE setComp #-}
   getComp = dComp
   {-# INLINE getComp #-}
+  repr = D
 
 instance Source D e where
-  unsafeIndex = INDEX_CHECK("(Source D ix e).unsafeIndex", size, dIndex)
+  unsafeIndex arr =
+    case dPrefIndex arr of
+      PrefIndex f -> f
+      PrefIndexLinear f -> f . toLinearIndex (size arr)
   {-# INLINE unsafeIndex #-}
-
-  unsafeOuterSlice !arr !szL !i = DArray (dComp arr) szL (unsafeIndex arr . consDim i)
+  unsafeLinearIndex arr =
+    case dPrefIndex arr of
+      PrefIndex f -> f . fromLinearIndex (size arr)
+      PrefIndexLinear f -> f
+  {-# INLINE unsafeLinearIndex #-}
+  unsafePrefIndex = dPrefIndex
+  {-# INLINE unsafePrefIndex #-}
+  unsafeOuterSlice !arr !szL !i =
+    makeArray (dComp arr) szL (unsafeIndex arr . consDim i)
   {-# INLINE unsafeOuterSlice #-}
-
   unsafeLinearSlice !o !sz arr =
-    DArray (dComp arr) sz $ \ !i -> unsafeIndex arr (fromLinearIndex (size arr) (i + o))
+    makeArrayLinear (dComp arr) sz $ \ !i -> unsafeLinearIndex arr (i + o)
   {-# INLINE unsafeLinearSlice #-}
 
-
 -- | /O(1)/ - Extract a portion of an array. Staring index and new size are
 -- not validated.
 unsafeExtract :: (Source r e, Index ix) => ix -> Sz ix -> Array r ix e -> Array D ix e
 unsafeExtract !sIx !newSz !arr =
-  DArray (getComp arr) newSz (unsafeIndex arr . liftIndex2 (+) sIx)
+  makeArray (getComp arr) newSz (unsafeIndex arr . liftIndex2 (+) sIx)
 {-# INLINE unsafeExtract #-}
 
 -- | /O(1)/ - Take a slice out of an array from within
-unsafeSlice :: (Source r e, Index ix, Index (Lower ix), MonadThrow m) =>
-  Array r ix e -> ix -> Sz ix -> Dim -> m (Array D (Lower ix) e)
+unsafeSlice ::
+     (Source r e, Index ix, Index (Lower ix), MonadThrow m)
+  => Array r ix e
+  -> ix
+  -> Sz ix
+  -> Dim
+  -> m (Array D (Lower ix) e)
 unsafeSlice arr start cut@(SafeSz cutSz) dim = do
   newSz <- dropDimM cutSz dim
   return $ unsafeResize (SafeSz newSz) (unsafeExtract start cut arr)
@@ -101,11 +121,10 @@
 -- | /O(1)/ - Take a slice out of an array from the inside
 unsafeInnerSlice ::
      (Source r e, Index ix) => Array r ix e -> Sz (Lower ix) -> Int -> Array D (Lower ix) e
-unsafeInnerSlice !arr szL !i = DArray (getComp arr) szL (unsafeIndex arr . (`snocDim` i))
+unsafeInnerSlice !arr szL !i =
+  DArray (getComp arr) szL $ PrefIndex (unsafeIndex arr . (`snocDim` i))
 {-# INLINE unsafeInnerSlice #-}
 
-
-
 instance (Eq e, Index ix) => Eq (Array D ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
@@ -115,12 +134,11 @@
   {-# INLINE compare #-}
 
 instance Functor (Array D ix) where
-  fmap f (DArray c sz g) = DArray c sz (f . g)
+  fmap f (DArray c sz g) = DArray c sz (fmap f g)
   {-# INLINE fmap #-}
-  (<$) e (DArray c sz _) = DArray c sz (const e)
+  (<$) e (DArray c sz g) = DArray c sz (e <$ g)
   {-# INLINE (<$) #-}
 
-
 instance Index ix => Applicative (Array D ix) where
   pure = singleton
   {-# INLINE pure #-}
@@ -151,18 +169,32 @@
   {-# INLINE length #-}
   elem e = A.any (e ==)
   {-# INLINE elem #-}
-  toList arr = build (\ c n -> foldrFB c n arr)
+  toList arr = build (\c n -> foldrFB c n arr)
   {-# INLINE toList #-}
 
-
 instance Index ix => Load D ix e where
-  makeArray = DArray
+  makeArray comp sz = DArray comp sz . PrefIndex
   {-# INLINE makeArray #-}
-  iterArrayLinearST_ !scheduler !arr =
-    splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
+  makeArrayLinear comp sz = DArray comp sz . PrefIndexLinear
+  {-# INLINE makeArrayLinear #-}
+  iterArrayLinearST_ !scheduler DArray {..} uWrite =
+    case dPrefIndex of
+      PrefIndex f ->
+        iterTargetFullST_ defRowMajor scheduler 0 dSize $ \ !i -> uWrite i . f
+      PrefIndexLinear f ->
+        iterTargetFullST_ defRowMajorLinear scheduler 0 dSize $ \ !i _ -> uWrite i (f i)
   {-# INLINE iterArrayLinearST_ #-}
 
-instance Index ix => StrideLoad D ix e
+instance Index ix => StrideLoad D ix e where
+  iterArrayLinearWithStrideST_ !scheduler !stride sz DArray {..} uWrite =
+    case dPrefIndex of
+      PrefIndex f ->
+        iterTargetFullWithStrideST_ defRowMajor scheduler 0 sz stride $ \i ->
+          uWrite i . f
+      PrefIndexLinear f -> do
+        iterTargetFullWithStrideST_ defRowMajor scheduler 0 sz stride $ \i ->
+          uWrite i . f . toLinearIndex dSize
+  {-# INLINE iterArrayLinearWithStrideST_ #-}
 
 instance Index ix => Stream D ix e where
   toStream = S.steps
@@ -173,11 +205,20 @@
 -- | Map an index aware function over an array
 --
 -- @since 0.1.0
-imap :: forall r ix e a. (Index ix, Source r e) => (ix -> e -> a) -> Array r ix e -> Array D ix a
-imap f !arr = DArray (getComp arr) (size arr) (\ !ix -> f ix (unsafeIndex arr ix))
+imap ::
+     forall r ix e a. (Index ix, Source r e)
+  => (ix -> e -> a)
+  -> Array r ix e
+  -> Array D ix a
+imap f !arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix -> DArray (getComp arr) sz $ PrefIndex (\ !ix -> f ix (gix ix))
+    PrefIndexLinear gi ->
+      DArray (getComp arr) sz $ PrefIndex (\ !ix -> f ix (gi (toLinearIndex sz ix)))
+  where
+    !sz = size arr
 {-# INLINE imap #-}
 
-
 instance Num e => FoldNumeric D e where
   unsafeDotProduct = defaultUnsafeDotProduct
   {-# INLINE unsafeDotProduct #-}
@@ -187,23 +228,20 @@
   {-# INLINE foldArray #-}
 
 instance Num e => Numeric D e where
-  unsafeLiftArray f arr = arr {dIndex = f . dIndex arr}
+  unsafeLiftArray f arr = arr {dPrefIndex = f <$> dPrefIndex arr}
   {-# INLINE unsafeLiftArray #-}
-  unsafeLiftArray2 f a1 a2 = -- TODO: possibly use the first size, it is unsafe anyways.
-    DArray (dComp a1 <> dComp a2) (SafeSz (liftIndex2 min (unSz (dSize a1)) (unSz (dSize a2)))) $ \i ->
-      f (dIndex a1 i) (dIndex a2 i)
+  unsafeLiftArray2 f a1 a2 = zipWithInternal (size a1) f a1 a2
   {-# INLINE unsafeLiftArray2 #-}
 
-
 instance Floating e => NumericFloat D e
 
-
-
 -- | /O(1)/ Conversion from a source array to `D` representation.
 delay :: (Index ix, Source r e) => Array r ix e -> Array D ix e
-delay arr = DArray (getComp arr) (size arr) (unsafeIndex arr)
+delay arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix -> makeArray (getComp arr) (size arr) gix
+    PrefIndexLinear gi -> makeArrayLinear (getComp arr) (size arr) gi
 {-# INLINE [1] delay #-}
-
 {-# RULES
 "delay" [~1] forall (arr :: Array D ix e) . delay arr = arr
  #-}
@@ -212,8 +250,12 @@
 -- element. Empty arrays are always equal, regardless of their size.
 --
 -- @since 0.5.7
-eqArrays :: (Index ix, Source r1 e1, Source r2 e2) =>
-            (e1 -> e2 -> Bool) -> Array r1 ix e1 -> Array r2 ix e2 -> Bool
+eqArrays ::
+     (Index ix, Source r1 e1, Source r2 e2)
+  => (e1 -> e2 -> Bool)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Bool
 eqArrays f arr1 arr2 =
   let sz1 = size arr1
       sz2 = size arr2
@@ -221,7 +263,7 @@
        not
          (A.any
             not
-            (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
+            (makeArray @D (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
                f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix)))) ||
       (isZeroSz sz1 && isZeroSz sz2)
 {-# INLINE eqArrays #-}
@@ -231,12 +273,16 @@
 -- you need an ordering but do not care about which one is used.
 --
 -- @since 0.5.7
-compareArrays :: (Index ix, Source r1 e1, Source r2 e2) =>
-       (e1 -> e2 -> Ordering) -> Array r1 ix e1 -> Array r2 ix e2 -> Ordering
+compareArrays ::
+     (Index ix, Source r1 e1, Source r2 e2)
+  => (e1 -> e2 -> Ordering)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Ordering
 compareArrays f arr1 arr2 =
   compare (size arr1) (size arr2) <>
   A.fold
-    (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
+    (makeArray @D (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
        f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
 {-# INLINE compareArrays #-}
 
@@ -244,27 +290,27 @@
 -- sizes.
 --
 -- @since 1.0.0
-liftArray2'
-  :: (HasCallStack, Index ix, Source r1 a, Source r2 b)
-  => (a -> b -> e) -> Array r1 ix a -> Array r2 ix b -> Array D ix e
+liftArray2' ::
+     (HasCallStack, Index ix, Source r1 a, Source r2 b)
+  => (a -> b -> e)
+  -> Array r1 ix a
+  -> Array r2 ix b
+  -> Array D ix e
 liftArray2' f arr1 arr2 = throwEither $ liftArray2M f arr1 arr2
 {-# INLINE liftArray2' #-}
 
-
 -- | Similar to `Data.Massiv.Array.zipWith`, except dimensions of both arrays
 -- have to be the same, otherwise it throws `SizeMismatchException`.
 --
 -- @since 1.0.0
-liftArray2M
-  :: (Index ix, Source r1 a, Source r2 b, MonadThrow m)
-  => (a -> b -> e) -> Array r1 ix a -> Array r2 ix b -> m (Array D ix e)
+liftArray2M ::
+     (Index ix, Source r1 a, Source r2 b, MonadThrow m)
+  => (a -> b -> e)
+  -> Array r1 ix a
+  -> Array r2 ix b
+  -> m (Array D ix e)
 liftArray2M f !arr1 !arr2
-  | sz1 == sz2 =
-    pure $
-    DArray
-      (getComp arr1 <> getComp arr2)
-      sz1
-      (\ !ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
+  | sz1 == sz2 = pure $ zipWithInternal sz1 f arr1 arr2
   | isZeroSz sz1 && isZeroSz sz2 = pure A.empty
   | otherwise = throwM $ SizeMismatchException (size arr1) (size arr2)
   where
@@ -272,3 +318,19 @@
     sz2 = size arr2
 {-# INLINE liftArray2M #-}
 
+zipWithInternal ::
+     (Index ix, Source r1 e1, Source r2 e2)
+  => Sz ix
+  -> (e1 -> e2 -> e3)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Array D ix e3
+zipWithInternal sz f arr1 arr2 =
+  case unsafePrefIndex arr1 of
+    PrefIndexLinear gi1
+      | PrefIndexLinear gi2 <- unsafePrefIndex arr2 ->
+        makeArrayLinear comp sz (\ !i -> f (gi1 i) (gi2 i))
+    _ -> makeArray comp sz (\ !ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
+  where
+    comp = getComp arr1 <> getComp arr2
+{-# INLINE zipWithInternal #-}
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
@@ -11,7 +11,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Push
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -61,6 +61,7 @@
   {-# INLINE getComp #-}
   setComp c arr = arr {dlComp = c}
   {-# INLINE setComp #-}
+  repr = DL
 
 
 instance Index ix => Shape DL ix where
@@ -81,8 +82,10 @@
 instance Monoid (Array DL Ix1 e) where
   mempty = DLArray {dlComp = mempty, dlSize = zeroSz, dlLoad = \_ _ _ _ -> pure ()}
   {-# INLINE mempty #-}
+#if !MIN_VERSION_base(4,11,0)
   mappend = mappendDL
   {-# INLINE mappend #-}
+#endif
   mconcat [] = mempty
   mconcat [x] = x
   mconcat [x, y] = x <> y
@@ -172,18 +175,7 @@
   -> (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 load
-  where
-    load :: forall s.
-      Scheduler s () -> Ix1 -> (Ix1 -> e -> ST s ()) -> (Ix1 -> Sz1 -> e -> ST s ()) -> ST s ()
-    load _scheduler !startAt uWrite uSet = do
-      uSet startAt (toLinearSz sz) defVal
-      let safeWrite !ix !e
-            | isSafeIndex sz ix = uWrite (startAt + toLinearIndex sz ix) e >> pure True
-            | otherwise = pure False
-          {-# INLINE safeWrite #-}
-      writer safeWrite
-    {-# INLINE load #-}
+makeLoadArrayS sz defVal writer = makeLoadArray Seq sz defVal (const writer)
 {-# INLINE makeLoadArrayS #-}
 
 -- | Specify how an array should be loaded into memory. Unlike `makeLoadArrayS`, loading
diff --git a/src/Data/Massiv/Array/Delayed/Stream.hs b/src/Data/Massiv/Array/Delayed/Stream.hs
--- a/src/Data/Massiv/Array/Delayed/Stream.hs
+++ b/src/Data/Massiv/Array/Delayed/Stream.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Stream
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -76,11 +76,11 @@
   isNull = S.unId . S.null . coerce
   {-# INLINE isNull #-}
 
-
---TODO remove
+-- | For now only `Seq` strategy.
 instance Strategy DS where
   getComp _ = Seq
   setComp _ = id
+  repr = DS
 
 
 instance Functor (Array DS Ix1) where
@@ -101,8 +101,6 @@
 #endif
 
 instance Monad (Array DS Ix1) where
-  return = fromSteps . S.singleton
-  {-# INLINE return #-}
   (>>=) arr f = coerce (S.concatMap (coerce . f) (dsArray arr))
   {-# INLINE (>>=) #-}
 
@@ -141,8 +139,10 @@
 instance Monoid (Array DS Ix1 e) where
   mempty = DSArray S.empty
   {-# INLINE mempty #-}
+#if !MIN_VERSION_base(4,11,0)
   mappend = (<>)
   {-# INLINE mappend #-}
+#endif
 
 instance IsList (Array DS Ix1 e) where
   type Item (Array DS Ix1 e) = e
diff --git a/src/Data/Massiv/Array/Delayed/Windowed.hs b/src/Data/Massiv/Array/Delayed/Windowed.hs
--- a/src/Data/Massiv/Array/Delayed/Windowed.hs
+++ b/src/Data/Massiv/Array/Delayed/Windowed.hs
@@ -10,7 +10,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Windowed
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -67,6 +67,7 @@
   {-# INLINE setComp #-}
   getComp = dComp . dwArray
   {-# INLINE getComp #-}
+  repr = DW
 
 
 instance Functor (Array DW ix) where
@@ -197,12 +198,12 @@
   -> Array DW Ix1 e
   -> (Ix1 -> e -> m a)
   -> m (Ix1 -> Ix1 -> m (), Ix1, Ix1)
-loadWithIx1 with (DWArray (DArray _ sz indexB) mWindow) uWrite = do
+loadWithIx1 with (DWArray a@(DArray _ sz _) mWindow) uWrite = do
   let Window it wk indexW _ = fromMaybe zeroWindow mWindow
       wEnd = it + unSz wk
-  with $ iterM_ 0 it 1 (<) $ \ !i -> uWrite i (indexB i)
-  with $ iterM_ wEnd (unSz sz) 1 (<) $ \ !i -> uWrite i (indexB i)
-  return (\from to -> with $ iterM_ from to 1 (<) $ \ !i -> uWrite i (indexW i), it, wEnd)
+  with $ iterA_ 0 it 1 (<) $ \ !i -> uWrite i (unsafeIndex a i)
+  with $ iterA_ wEnd (unSz sz) 1 (<) $ \ !i -> uWrite i (unsafeIndex a i)
+  return (\from to -> with $ iterA_ from to 1 (<) $ \ !i -> uWrite i (indexW i), it, wEnd)
 {-# INLINE loadWithIx1 #-}
 
 instance Index ix => Shape DW ix where
@@ -219,7 +220,7 @@
   iterArrayLinearST_ scheduler arr uWrite = do
     (loadWindow, wStart, wEnd) <- loadWithIx1 (scheduleWork scheduler) arr uWrite
     let (chunkWidth, slackWidth) = (wEnd - wStart) `quotRem` numWorkers scheduler
-    loopM_ 0 (< numWorkers scheduler) (+ 1) $ \ !wid ->
+    loopA_ 0 (< numWorkers scheduler) (+ 1) $ \ !wid ->
       let !it' = wid * chunkWidth + wStart
        in loadWindow it' (it' + chunkWidth)
     when (slackWidth > 0) $
@@ -231,7 +232,7 @@
   iterArrayLinearWithStrideST_ scheduler stride sz arr uWrite = do
       (loadWindow, (wStart, wEnd)) <- loadArrayWithIx1 (scheduleWork scheduler) arr stride sz uWrite
       let (chunkWidth, slackWidth) = (wEnd - wStart) `quotRem` numWorkers scheduler
-      loopM_ 0 (< numWorkers scheduler) (+ 1) $ \ !wid ->
+      loopA_ 0 (< numWorkers scheduler) (+ 1) $ \ !wid ->
         let !it' = wid * chunkWidth + wStart
          in loadWindow (it', it' + chunkWidth)
       when (slackWidth > 0) $
@@ -247,18 +248,18 @@
   -> Sz1
   -> (Ix1 -> e -> m a)
   -> m ((Ix1, Ix1) -> m (), (Ix1, Ix1))
-loadArrayWithIx1 with (DWArray (DArray _ arrSz indexB) mWindow) stride _ uWrite = do
+loadArrayWithIx1 with (DWArray darr@(DArray _ arrSz _) mWindow) stride _ uWrite = do
   let Window it wk indexW _ = fromMaybe zeroWindow mWindow
       wEnd = it + unSz wk
       strideIx = unStride stride
-  with $ iterM_ 0 it strideIx (<) $ \ !i -> uWrite (i `div` strideIx) (indexB i)
+  with $ iterA_ 0 it strideIx (<) $ \ !i -> uWrite (i `div` strideIx) (unsafeIndex darr i)
   with $
-    iterM_ (strideStart stride wEnd) (unSz arrSz) strideIx (<) $ \ !i ->
-      uWrite (i `div` strideIx) (indexB i)
+    iterA_ (strideStart stride wEnd) (unSz arrSz) strideIx (<) $ \ !i ->
+      uWrite (i `div` strideIx) (unsafeIndex darr i)
   return
     ( \(from, to) ->
         with $
-        iterM_ (strideStart stride from) to strideIx (<) $ \ !i ->
+        iterA_ (strideStart stride from) to strideIx (<) $ \ !i ->
           uWrite (i `div` strideIx) (indexW i)
     , (it, wEnd))
 {-# INLINE loadArrayWithIx1 #-}
@@ -272,20 +273,21 @@
   -> (Int -> t1 -> m ())
   -> m (Ix2 -> m (), Ix2)
 loadWithIx2 with arr uWrite = do
-  let DWArray (DArray _ (Sz (m :. n)) indexB) window = arr
-  let Window (it :. jt) (Sz (wm :. wn)) indexW mUnrollHeight = fromMaybe zeroWindow window
-  let ib :. jb = (wm + it) :. (wn + jt)
+  let DWArray darr window = arr
+      Sz (m :. n) = dSize darr
+      Window (it :. jt) (Sz (wm :. wn)) indexW mUnrollHeight = fromMaybe zeroWindow window
+      ib :. jb = (wm + it) :. (wn + jt)
       !blockHeight = maybe 1 (min 7 . max 1) mUnrollHeight
       stride = oneStride
       !sz = strideSize stride $ outerSize arr
-      writeB !ix = uWrite (toLinearIndex sz ix) (indexB ix)
+      writeB !ix = uWrite (toLinearIndex sz ix) (unsafeIndex darr ix)
       {-# INLINE writeB #-}
       writeW !ix = uWrite (toLinearIndex sz ix) (indexW ix)
       {-# INLINE writeW #-}
-  with $ iterM_ (0 :. 0) (it :. n) (1 :. 1) (<) writeB
-  with $ iterM_ (ib :. 0) (m :. n) (1 :. 1) (<) writeB
-  with $ iterM_ (it :. 0) (ib :. jt) (1 :. 1) (<) writeB
-  with $ iterM_ (it :. jb) (ib :. n) (1 :. 1) (<) writeB
+  with $ iterA_ (0 :. 0) (it :. n) (1 :. 1) (<) writeB
+  with $ iterA_ (ib :. 0) (m :. n) (1 :. 1) (<) writeB
+  with $ iterA_ (it :. 0) (ib :. jt) (1 :. 1) (<) writeB
+  with $ iterA_ (it :. jb) (ib :. n) (1 :. 1) (<) writeB
   let f (it' :. ib') = with $ unrollAndJam blockHeight (it' :. jt) (ib' :. jb) 1 writeW
       {-# INLINE f #-}
   return (f, it :. ib)
@@ -300,22 +302,23 @@
   -> (Int -> e -> m ())
   -> m (Ix2 -> m (), Ix2)
 loadArrayWithIx2 with arr stride sz uWrite = do
-  let DWArray (DArray _ (Sz (m :. n)) indexB) window = arr
-  let Window (it :. jt) (Sz (wm :. wn)) indexW mUnrollHeight = fromMaybe zeroWindow window
-  let ib :. jb = (wm + it) :. (wn + jt)
+  let DWArray darr window = arr
+      Sz (m :. n) = dSize darr
+      Window (it :. jt) (Sz (wm :. wn)) indexW mUnrollHeight = fromMaybe zeroWindow window
+      ib :. jb = (wm + it) :. (wn + jt)
       !blockHeight = maybe 1 (min 7 . max 1) mUnrollHeight
       strideIx@(is :. js) = unStride stride
-      writeB !ix = uWrite (toLinearIndexStride stride sz ix) (indexB ix)
+      writeB !ix = uWrite (toLinearIndexStride stride sz ix) (unsafeIndex darr ix)
       {-# INLINE writeB #-}
       writeW !ix = uWrite (toLinearIndexStride stride sz ix) (indexW ix)
       {-# INLINE writeW #-}
-  with $ iterM_ (0 :. 0) (it :. n) strideIx (<) writeB
-  with $ iterM_ (strideStart stride (ib :. 0)) (m :. n) strideIx (<) writeB
-  with $ iterM_ (strideStart stride (it :. 0)) (ib :. jt) strideIx (<) writeB
-  with $ iterM_ (strideStart stride (it :. jb)) (ib :. n) strideIx (<) writeB
+  with $ iterA_ (0 :. 0) (it :. n) strideIx (<) writeB
+  with $ iterA_ (strideStart stride (ib :. 0)) (m :. n) strideIx (<) writeB
+  with $ iterA_ (strideStart stride (it :. 0)) (ib :. jt) strideIx (<) writeB
+  with $ iterA_ (strideStart stride (it :. jb)) (ib :. n) strideIx (<) writeB
   let f (it' :. ib')
         | is > 1 || blockHeight <= 1 -- Turn off unrolling for vertical strides
-         = iterM_ (strideStart stride (it' :. jt)) (ib' :. jb) strideIx (<) writeW
+         = iterA_ (strideStart stride (it' :. jt)) (ib' :. jb) strideIx (<) writeW
         | otherwise =
           unrollAndJam blockHeight (strideStart stride (it' :. jt)) (ib' :. jb) js writeW
       {-# INLINE f #-}
@@ -326,7 +329,7 @@
 loadWindowIx2 :: Monad m => Int -> (Ix2 -> m ()) -> Ix2 -> m ()
 loadWindowIx2 nWorkers loadWindow (it :. ib) = do
   let !(chunkHeight, slackHeight) = (ib - it) `quotRem` nWorkers
-  loopM_ 0 (< nWorkers) (+ 1) $ \ !wid ->
+  loopA_ 0 (< nWorkers) (+ 1) $ \ !wid ->
     let !it' = wid * chunkHeight + it
      in loadWindow (it' :. (it' + chunkHeight))
   when (slackHeight > 0) $
@@ -370,9 +373,8 @@
   -> ST s ()
 loadArrayWithIxN scheduler stride szResult arr uWrite = do
   let DWArray darr window = arr
-      DArray {dSize = szSource, dIndex = indexBorder} = darr
       Window {windowStart, windowSize, windowIndex, windowUnrollIx2} = fromMaybe zeroWindow window
-      !(headSourceSize, lowerSourceSize) = unconsSz szSource
+      !(headSourceSize, lowerSourceSize) = unconsSz (dSize darr)
       !lowerSize = snd $ unconsSz szResult
       !(s, lowerStrideIx) = unconsDim $ unStride stride
       !(curWindowStart, lowerWindowStart) = unconsDim windowStart
@@ -388,7 +390,12 @@
           }
       mkLowerArray mw i =
         DWArray
-          {dwArray = DArray Seq lowerSourceSize (indexBorder . consDim i), dwWindow = ($ i) <$> mw}
+          { dwArray = darr { dComp = Seq
+                           , dSize = lowerSourceSize
+                           , dPrefIndex = PrefIndex (unsafeIndex darr . consDim i)
+                           }
+          , dwWindow = ($ i) <$> mw
+          }
       loadLower mw !i =
         iterArrayLinearWithStrideST_
           scheduler
@@ -397,13 +404,13 @@
           (mkLowerArray mw i)
           (\k -> uWrite (k + pageElements * (i `div` s)))
       {-# NOINLINE loadLower #-}
-  loopM_ 0 (< headDim windowStart) (+ s) (loadLower Nothing)
-  loopM_
+  loopA_ 0 (< headDim windowStart) (+ s) (loadLower Nothing)
+  loopA_
     (strideStart (Stride s) curWindowStart)
     (< curWindowEnd)
     (+ s)
     (loadLower (Just mkLowerWindow))
-  loopM_ (strideStart (Stride s) curWindowEnd) (< unSz headSourceSize) (+ s) (loadLower Nothing)
+  loopA_ (strideStart (Stride s) curWindowEnd) (< unSz headSourceSize) (+ s) (loadLower Nothing)
 {-# INLINE loadArrayWithIxN #-}
 
 
@@ -416,9 +423,8 @@
   -> ST s ()
 loadWithIxN scheduler arr uWrite = do
   let DWArray darr window = arr
-      DArray {dSize = sz, dIndex = indexBorder} = darr
       Window {windowStart, windowSize, windowIndex, windowUnrollIx2} = fromMaybe zeroWindow window
-      !(si, szL) = unconsSz sz
+      !(si, szL) = unconsSz (dSize darr)
       !windowEnd = liftIndex2 (+) windowStart (unSz windowSize)
       !(t, windowStartL) = unconsDim windowStart
       !pageElements = totalElem szL
@@ -430,14 +436,18 @@
           , windowUnrollIx2 = windowUnrollIx2
           }
       mkLowerArray mw i =
-        DWArray {dwArray = DArray Seq szL (indexBorder . consDim i), dwWindow = ($ i) <$> mw}
+        DWArray
+          { dwArray =
+              darr {dComp = Seq, dSize = szL, dPrefIndex = PrefIndex (unsafeIndex darr . consDim i)}
+          , dwWindow = ($ i) <$> mw
+          }
       loadLower mw !i =
         scheduleWork_ scheduler $
         iterArrayLinearST_ scheduler (mkLowerArray mw i) (\k -> uWrite (k + pageElements * i))
       {-# NOINLINE loadLower #-}
-  loopM_ 0 (< headDim windowStart) (+ 1) (loadLower Nothing)
-  loopM_ t (< headDim windowEnd) (+ 1) (loadLower (Just mkLowerWindow))
-  loopM_ (headDim windowEnd) (< unSz si) (+ 1) (loadLower Nothing)
+  loopA_ 0 (< headDim windowStart) (+ 1) (loadLower Nothing)
+  loopA_ t (< headDim windowEnd) (+ 1) (loadLower (Just mkLowerWindow))
+  loopA_ (headDim windowEnd) (< unSz si) (+ 1) (loadLower Nothing)
 {-# INLINE loadWithIxN #-}
 
 
@@ -465,11 +475,11 @@
              6 -> f6
              _ -> f7
   let !ibS = ib - ((ib - it) `mod` bH)
-  loopM_ it (< ibS) (+ bH) $ \ !i ->
-    loopM_ jt (< jb) (+ js) $ \ !j ->
+  loopA_ it (< ibS) (+ bH) $ \ !i ->
+    loopA_ jt (< jb) (+ js) $ \ !j ->
       f' (i :. j)
-  loopM_ ibS (< ib) (+ 1) $ \ !i ->
-    loopM_ jt (< jb) (+ js) $ \ !j ->
+  loopA_ ibS (< ib) (+ 1) $ \ !i ->
+    loopA_ jt (< jb) (+ js) $ \ !j ->
       f (i :. j)
 {-# INLINE unrollAndJam #-}
 
diff --git a/src/Data/Massiv/Array/Manifest.hs b/src/Data/Massiv/Array/Manifest.hs
--- a/src/Data/Massiv/Array/Manifest.hs
+++ b/src/Data/Massiv/Array/Manifest.hs
@@ -8,7 +8,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -17,6 +17,33 @@
 module Data.Massiv.Array.Manifest
   ( -- * Manifest
     Manifest
+  -- ** Generate
+  , generateArray
+  , generateArrayLinear
+  , generateArrayS
+  , generateArrayLinearS
+  , generateSplitSeedArray
+  -- ** Stateful worker threads
+  , generateArrayWS
+  , generateArrayLinearWS
+  -- ** Unfold
+  , unfoldrPrimM_
+  , iunfoldrPrimM_
+  , unfoldrPrimM
+  , iunfoldrPrimM
+  , unfoldlPrimM_
+  , iunfoldlPrimM_
+  , unfoldlPrimM
+  , iunfoldlPrimM
+  -- ** Mapping
+  , forPrimM
+  , forPrimM_
+  , iforPrimM
+  , iforPrimM_
+  , iforLinearPrimM
+  , iforLinearPrimM_
+  , for2PrimM_
+  , ifor2PrimM_
   -- * Boxed
   , B(..)
   , BL(..)
diff --git a/src/Data/Massiv/Array/Manifest/Boxed.hs b/src/Data/Massiv/Array/Manifest/Boxed.hs
--- a/src/Data/Massiv/Array/Manifest/Boxed.hs
+++ b/src/Data/Massiv/Array/Manifest/Boxed.hs
@@ -11,7 +11,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Boxed
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -60,9 +60,10 @@
 import Control.Monad ((>=>))
 import Control.Monad.Primitive
 import qualified Data.Foldable as F (Foldable(..))
+import Data.Massiv.Array.Delayed.Pull (D)
 import Data.Massiv.Array.Delayed.Push (DL)
 import Data.Massiv.Array.Delayed.Stream (DS)
-import Data.Massiv.Array.Manifest.Internal (computeAs)
+import Data.Massiv.Array.Manifest.Internal (compute, computeAs)
 import Data.Massiv.Array.Manifest.List as L
 import Data.Massiv.Array.Mutable
 import Data.Massiv.Array.Ops.Fold
@@ -78,6 +79,9 @@
 import GHC.Exts as GHC
 import Prelude hiding (mapM, replicate)
 import System.IO.Unsafe (unsafePerformIO)
+#if !MIN_VERSION_vector(0,13,0)
+import Unsafe.Coerce (unsafeCoerce)
+#endif
 
 #include "massiv.h"
 
@@ -149,6 +153,7 @@
   {-# INLINE setComp #-}
   getComp = blComp
   {-# INLINE getComp #-}
+  repr = BL
 
 
 instance Source BL e where
@@ -216,6 +221,9 @@
   {-# INLINE maxLinearSize #-}
 
 instance Index ix => Load BL ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
+
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
@@ -334,6 +342,7 @@
   {-# INLINE getComp #-}
   setComp c arr = coerceBoxedArray (coerce arr) { blComp = c }
   {-# INLINE setComp #-}
+  repr = B
 
 
 instance Index ix => Shape B ix where
@@ -383,6 +392,9 @@
   {-# INLINE unsafeLinearWrite #-}
 
 instance Index ix => Load B ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
+
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
@@ -492,11 +504,12 @@
   compare = compareArrays compare
   {-# INLINE compare #-}
 
-instance Strategy N where
+instance Strategy BN where
   setComp c = coerce (setComp c)
   {-# INLINE setComp #-}
   getComp = blComp . coerce
   {-# INLINE getComp #-}
+  repr = BN
 
 instance NFData e => Source BN e where
   unsafeLinearIndex (BNArray arr) = unsafeLinearIndex arr
@@ -553,6 +566,8 @@
   {-# INLINE unsafeLinearWrite #-}
 
 instance (Index ix, NFData e) => Load BN ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
   replicate comp sz e = runST (newMArray sz e >>= unsafeFreeze comp)
@@ -741,7 +756,7 @@
   -> m (MArray (PrimState m) BL Ix1 e)
 fromMutableArraySeq with ma = do
   let !sz = A.sizeofMutableArray ma
-  loopM_ 0 (< sz) (+ 1) (A.readArray ma >=> (`with` return ()))
+  loopA_ 0 (< sz) (+ 1) (A.readArray ma >=> (`with` return ()))
   return $! MBLArray (SafeSz sz) 0 ma
 {-# INLINE fromMutableArraySeq #-}
 
@@ -770,14 +785,27 @@
 evalNormalForm (BArray arr) = arr `deepseqArray` BNArray arr
 {-# INLINE evalNormalForm #-}
 
--- | /O(1)/ - Converts a boxed `Array` into a `VB.Vector`.
+-- | /O(1)/ - Converts a boxed `Array` into a `VB.Vector` without touching any
+-- elements.
 --
 -- @since 0.5.0
-toBoxedVector :: Index ix => Array BL ix a -> VB.Vector a
-toBoxedVector arr = runST $ VB.unsafeFreeze . toBoxedMVector =<< unsafeThaw arr
 {-# INLINE toBoxedVector #-}
+toBoxedVector :: Index ix => Array BL ix a -> VB.Vector a
+toBoxedVector BLArray{blOffset = off, blSize = sz, blData = arr } =
+#if MIN_VERSION_vector(0,13,0)
+  VB.unsafeFromArraySlice arr off (totalElem sz)
+#elif MIN_VERSION_vector(0,12,2)
+  VB.unsafeTake (totalElem sz) (VB.unsafeDrop off (VB.fromArray arr))
+#else
+  fromVectorCast $ VectorCast off (totalElem sz) arr
 
--- | /O(1)/ - Converts a boxed `MArray` into a `VMB.MVector`.
+fromVectorCast :: VectorCast a -> VB.Vector a
+fromVectorCast = unsafeCoerce
+#endif
+
+
+
+-- | /O(1)/ - Converts a boxed `MArray` into a `MVB.MVector`.
 --
 -- @since 0.5.0
 toBoxedMVector :: Index ix => MArray s BL ix a -> MVB.MVector s a
@@ -800,7 +828,7 @@
 evalBoxedMVector :: PrimMonad m => MVB.MVector (PrimState m) a -> m (MArray (PrimState m) B Ix1 a)
 evalBoxedMVector (MVB.MVector o k ma) =
   let marr = MBArray (MBLArray (SafeSz k) o ma)
-   in marr <$ loopM_ o (< k) (+ 1) (A.readArray ma >=> (`seq` pure ()))
+   in marr <$ loopA_ o (< k) (+ 1) (A.readArray ma >=> (`seq` pure ()))
 {-# INLINE evalBoxedMVector #-}
 
 
@@ -808,13 +836,25 @@
 --
 -- @since 0.6.0
 fromBoxedVector :: VB.Vector a -> Vector BL a
-fromBoxedVector v =
-  runST $ do
-    MVB.MVector o k ma <- VB.unsafeThaw v
-    unsafeFreeze Seq $ MBLArray (SafeSz k) o ma
 {-# INLINE fromBoxedVector #-}
+fromBoxedVector v =
+  BLArray {blComp = Seq, blSize = SafeSz n, blOffset = offset, blData = arr}
+  where
+#if MIN_VERSION_vector(0,13,0)
+    (arr, offset, n) = VB.toArraySlice v
+#else
+    VectorCast offset n arr = toVectorCast v
 
+-- This internal type is needed to get into the internals of a boxed vector,
+-- since it is not possible until vector-0.13 version.
+data VectorCast a =
+  VectorCast {-# UNPACK #-}!Int {-# UNPACK #-}!Int {-# UNPACK #-}!(A.Array a)
 
+toVectorCast :: VB.Vector a -> VectorCast a
+toVectorCast = unsafeCoerce
+#endif
+
+
 -- | /O(1)/ - Convert mutable boxed vector to a lazy mutable boxed array. Both keep
 -- pointing to the same memory
 --
@@ -849,7 +889,7 @@
      (NFData a, PrimMonad m) => MVB.MVector (PrimState m) a -> m (MArray (PrimState m) N Ix1 a)
 evalNormalBoxedMVector (MVB.MVector o k ma) =
   let marr = MBNArray (MBLArray (SafeSz k) o ma)
-   in marr <$ loopM_ o (< k) (+ 1) (A.readArray ma >=> (`deepseq` pure ()))
+   in marr <$ loopA_ o (< k) (+ 1) (A.readArray ma >=> pure . rnf)
 {-# INLINE evalNormalBoxedMVector #-}
 
 -- | /O(n)/ - Convert a boxed vector and evaluate all elements to WHNF. Computation
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
@@ -10,7 +10,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Manifest/List.hs b/src/Data/Massiv/Array/Manifest/List.hs
--- a/src/Data/Massiv/Array/Manifest/List.hs
+++ b/src/Data/Massiv/Array/Manifest/List.hs
@@ -6,7 +6,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.List
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Manifest/Primitive.hs b/src/Data/Massiv/Array/Manifest/Primitive.hs
--- a/src/Data/Massiv/Array/Manifest/Primitive.hs
+++ b/src/Data/Massiv/Array/Manifest/Primitive.hs
@@ -11,7 +11,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Primitive
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -112,6 +112,7 @@
   {-# INLINE getComp #-}
   setComp c arr = arr { pComp = c }
   {-# INLINE setComp #-}
+  repr = P
 
 
 instance Index ix => Shape P ix where
@@ -206,6 +207,8 @@
 
 
 instance (Prim e, Index ix) => Load P ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
diff --git a/src/Data/Massiv/Array/Manifest/Storable.hs b/src/Data/Massiv/Array/Manifest/Storable.hs
--- a/src/Data/Massiv/Array/Manifest/Storable.hs
+++ b/src/Data/Massiv/Array/Manifest/Storable.hs
@@ -9,7 +9,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Storable
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -40,7 +40,7 @@
 import Control.Monad
 import Control.Monad.IO.Unlift
 import Control.Monad.Primitive
-import Data.Massiv.Array.Delayed.Pull (compareArrays, eqArrays)
+import Data.Massiv.Array.Delayed.Pull (D, compareArrays, eqArrays)
 import Data.Massiv.Array.Manifest.Internal
 import Data.Massiv.Array.Manifest.List as A
 import Data.Massiv.Array.Mutable
@@ -103,6 +103,7 @@
   {-# INLINE getComp #-}
   setComp c arr = arr { sComp = c }
   {-# INLINE setComp #-}
+  repr = S
 
 plusFp :: ForeignPtr a -> Int -> ForeignPtr b
 plusFp (ForeignPtr addr c) (I# d) = ForeignPtr (plusAddr# addr d) c
@@ -185,8 +186,8 @@
     INDEX_CHECK("(Manifest S ix e).unsafeLinearWrite", const (toLinearSz _sz), (\_ _ -> unsafeWithForeignPtr fp (\p -> pokeElemOff p o e))) fp o
   {-# INLINE unsafeLinearWrite #-}
 
-  unsafeLinearSet (MSArray _ fp) i k =
-    MVG.basicSet (MVS.unsafeFromForeignPtr0 (advanceForeignPtr fp i) (unSz k))
+  unsafeLinearSet (MSArray _ fp) i k e =
+    stToPrim (MVG.basicSet (MVS.unsafeFromForeignPtr0 (advanceForeignPtr fp i) (unSz k)) e)
   {-# INLINE unsafeLinearSet #-}
 
   unsafeLinearCopy (MSArray _ fpFrom) iFrom (MSArray _ fpTo) iTo (Sz k) = do
@@ -218,6 +219,9 @@
   {-# INLINE unsafeLinearShrink #-}
 
 instance (Index ix, Storable e) => Load S ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
+
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
diff --git a/src/Data/Massiv/Array/Manifest/Unboxed.hs b/src/Data/Massiv/Array/Manifest/Unboxed.hs
--- a/src/Data/Massiv/Array/Manifest/Unboxed.hs
+++ b/src/Data/Massiv/Array/Manifest/Unboxed.hs
@@ -8,7 +8,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Unboxed
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -26,15 +26,17 @@
   ) where
 
 import Control.DeepSeq (NFData(..), deepseq)
-import Data.Massiv.Array.Delayed.Pull (eqArrays, compareArrays)
+import Control.Monad.Primitive (stToPrim)
+import Data.Massiv.Array.Delayed.Pull (D, compareArrays, eqArrays)
+import Data.Massiv.Array.Manifest.Internal
 import Data.Massiv.Array.Manifest.List as A
-import Data.Massiv.Vector.Stream as S (steps, isteps)
 import Data.Massiv.Array.Mutable
 import Data.Massiv.Core.Common
 import Data.Massiv.Core.List
 import Data.Massiv.Core.Operations
-import Data.Vector.Unboxed (Unbox)
+import Data.Massiv.Vector.Stream as S (isteps, steps)
 import qualified Data.Vector.Generic.Mutable as VGM
+import Data.Vector.Unboxed (Unbox)
 import qualified Data.Vector.Unboxed as VU
 import qualified Data.Vector.Unboxed.Mutable as MVU
 import GHC.Exts as GHC (IsList(..))
@@ -69,6 +71,7 @@
   {-# INLINE getComp #-}
   setComp c arr = arr { uComp = c }
   {-# INLINE setComp #-}
+  repr = U
 
 
 instance (Unbox e, Eq e, Index ix) => Eq (Array U ix e) where
@@ -104,6 +107,9 @@
   {-# INLINE unsafeResize #-}
 
 instance (Unbox e, Index ix) => Load U ix e where
+  makeArray comp sz f = compute (makeArray comp sz f :: Array D ix e)
+  {-# INLINE makeArray #-}
+
   makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
@@ -144,7 +150,7 @@
   unsafeNew sz = MUArray sz <$> MVU.unsafeNew (totalElem sz)
   {-# INLINE unsafeNew #-}
 
-  initialize (MUArray _ marr) = VGM.basicInitialize marr
+  initialize (MUArray _ marr) = stToPrim (VGM.basicInitialize marr)
   {-# INLINE initialize #-}
 
   unsafeLinearCopy (MUArray _ mvFrom) iFrom (MUArray _ mvTo) iTo (Sz k) =
diff --git a/src/Data/Massiv/Array/Manifest/Vector.hs b/src/Data/Massiv/Array/Manifest/Vector.hs
--- a/src/Data/Massiv/Array/Manifest/Vector.hs
+++ b/src/Data/Massiv/Array/Manifest/Vector.hs
@@ -6,7 +6,7 @@
 {-# LANGUAGE TypeOperators #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Vector
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -7,7 +7,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -61,6 +61,7 @@
   , generateArrayLinear
   , generateArrayS
   , generateArrayLinearS
+  , generateSplitSeedArray
   -- *** Stateful worker threads
   , generateArrayWS
   , generateArrayLinearWS
@@ -108,15 +109,17 @@
 
 -- TODO: add fromListM, et al.
 
-import Data.Maybe (fromMaybe)
-import Control.Monad (void, when, unless, (>=>))
-import Control.Monad.ST
+import Control.Monad (unless, void, when, (>=>))
 import Control.Monad.Primitive
+import Control.Monad.ST
 import Control.Scheduler
-import Data.Massiv.Core.Common
-import Data.Massiv.Array.Mutable.Internal
+import Data.IORef
 import Data.Massiv.Array.Delayed.Pull (D)
+import Data.Massiv.Array.Mutable.Internal
+import Data.Massiv.Core.Common
+import Data.Maybe (fromMaybe)
 import Prelude hiding (mapM, read)
+import System.IO.Unsafe (unsafePerformIO)
 
 -- | /O(1)/ - Change the size of a mutable array. Throws
 -- `SizeElementsMismatchException` if total number of elements does not match
@@ -285,7 +288,7 @@
     marr <- unsafeNew sz
     withMassivScheduler_ (getComp arr) $ \scheduler ->
       splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork_ scheduler $ unsafeArrayLinearCopy arr start marr start (SafeSz chunkLength)
         let slackLength = totalLength - slackStart
         when (slackLength > 0) $
@@ -346,7 +349,7 @@
     tmarr <- unsafeNew sz
     withMassivScheduler_ comp $ \scheduler ->
       splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork_ scheduler $ unsafeLinearCopy smarr start tmarr start (SafeSz chunkLength)
         let slackLength = totalLength - slackStart
         when (slackLength > 0) $
@@ -445,7 +448,7 @@
   -> m (MArray (PrimState m) r ix e)
 makeMArrayLinearS sz f = do
   marr <- unsafeNew sz
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\ !i -> f i >>= unsafeLinearWrite marr i)
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\ !i -> f i >>= unsafeLinearWrite marr i)
   return marr
 {-# INLINE makeMArrayLinearS #-}
 
@@ -615,8 +618,8 @@
 -- @since 0.2.6
 generateArrayS ::
      forall r ix e m. (Manifest r e, Index ix, PrimMonad m)
-  => Sz ix -- ^ Resulting size of the array
-  -> (ix -> m e) -- ^ Element producing generator
+  => Sz ix -- ^ Size of the array
+  -> (ix -> m e) -- ^ Element producing action
   -> m (Array r ix e)
 generateArrayS sz gen = generateArrayLinearS sz (gen . fromLinearIndex sz)
 {-# INLINE generateArrayS #-}
@@ -631,7 +634,7 @@
   -> m (Array r ix e)
 generateArrayLinearS sz gen = do
   marr <- unsafeNew sz
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) $ \i -> gen i >>= unsafeLinearWrite marr i
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) $ \i -> gen i >>= unsafeLinearWrite marr i
   unsafeFreeze Seq marr
 {-# INLINE generateArrayLinearS #-}
 
@@ -657,12 +660,65 @@
      forall r ix e m. (MonadUnliftIO m, Manifest r e, Index ix)
   => Comp
   -> Sz ix
-  -> (Int -> m e)
+  -> (Ix1 -> m e)
   -> m (Array r ix e)
 generateArrayLinear comp sz f = makeMArrayLinear comp sz f >>= liftIO . unsafeFreeze comp
 {-# INLINE generateArrayLinear #-}
 
 
+-- | Similar to `Data.Massiv.Array.makeSplitSeedArray`, except it will produce a
+-- Manifest array and will return back the last unused seed together with all
+-- final seeds produced by each scheduled job. This function can be thought of
+-- as an unfolding done in parallel while iterating in a customizable manner.
+--
+-- @since 1.0.2
+generateSplitSeedArray ::
+     forall r ix e g it. (Iterator it, Manifest r e, Index ix)
+  => it -- ^ Iterator
+  -> g -- ^ Initial seed
+  -> (forall s. g -> ST s (g, g))
+     -- ^ An ST action that can split a seed into two independent seeds. It will
+     -- be called the same number of times as the number of jobs that will get
+     -- scheduled during parallelization. Eg. only once for the sequential case.
+  -> Comp -- ^ Computation strategy.
+  -> Sz ix -- ^ Resulting size of the array.
+  -> (forall s. Ix1 -> ix -> g -> ST s (e, g))
+     -- ^ An ST action that produces a value and the next seed. It takes both
+     -- versions of the index, in linear and in multi-dimensional forms, as well
+     -- as the current seeding value. Returns the element for the array cell
+     -- together with the new seed that will be used for the next element
+     -- generation
+  -> (g, [g], Array r ix e)
+  -- ^ Returned values are:
+  --
+  -- * The final split of the supplied seed.
+  --
+  -- * Results of scheduled jobs in the same order that they where scheduled
+  --
+  -- * Final array that was fully filled using the supplied action and iterator.
+generateSplitSeedArray it seed splitSeed comp sz genFunc =
+  unsafePerformIO $ do
+    marr <- unsafeNew sz
+    ref <- newIORef Nothing
+    res <- withSchedulerR comp $ \ scheduler -> do
+      fin <- stToIO $
+        iterTargetFullAccST it scheduler 0 sz seed splitSeed $ \ !i ix !g ->
+          genFunc i ix g >>= \ (x, g') -> g' <$ unsafeLinearWrite marr i x
+      writeIORef ref $ Just fin
+    mFin <- readIORef ref
+    case res of
+      Finished gs |
+        Just fin <- mFin -> do
+          arr <- unsafeFreeze comp marr
+          pure (fin, gs, arr)
+      -- This case does not make much sence for array filling and can only
+      -- happen with a custom 'Iterator' defined outside massiv, therefore it is
+      -- ok to not support it.
+      _ -> error $ "Parallelized array filling finished prematurely. " ++
+           "This feature is not supported by the 'generateSplitSeedArray' function."
+{-# INLINE generateSplitSeedArray #-}
+
+
 -- | Same as `generateArrayWS`, but use linear indexing instead.
 --
 -- @since 0.3.4
@@ -701,21 +757,20 @@
 --
 -- ====__Examples__
 --
--- Create an array with Fibonacci numbers while performing and `IO` action on the accumulator for
--- each element of the array.
+-- Create an array with Fibonacci numbers while performing an `IO` action at each iteration.
 --
 -- >>> import Data.Massiv.Array
--- >>> unfoldrPrimM_ (Sz1 10) (\a@(f0, f1) -> let fn = f0 + f1 in print a >> return (f0, (f1, fn))) (0, 1) :: IO (Array P Ix1 Int)
--- (0,1)
--- (1,1)
--- (1,2)
--- (2,3)
--- (3,5)
--- (5,8)
--- (8,13)
--- (13,21)
--- (21,34)
--- (34,55)
+-- >>> unfoldrPrimM_ (Sz1 10) (\(f0, f1) -> (f0, (f1, f0 + f1)) <$ print f1) (0, 1) :: IO (Array P Ix1 Int)
+-- 1
+-- 1
+-- 2
+-- 3
+-- 5
+-- 8
+-- 13
+-- 21
+-- 34
+-- 55
 -- Array P Seq (Sz1 10)
 --   [ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34 ]
 --
@@ -865,7 +920,7 @@
 -- @since 0.4.0
 forPrimM_ :: (Manifest r e, Index ix, PrimMonad m) => MArray (PrimState m) r ix e -> (e -> m ()) -> m ()
 forPrimM_ marr f =
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+1) (unsafeLinearRead marr >=> f)
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+1) (unsafeLinearRead marr >=> f)
 {-# INLINE forPrimM_ #-}
 
 -- | Sequentially loop over a mutable array while modifying each element with an action.
@@ -873,7 +928,7 @@
 -- @since 0.4.0
 forPrimM :: (Manifest r e, Index ix, PrimMonad m) => MArray (PrimState m) r ix e -> (e -> m e) -> m ()
 forPrimM marr f =
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+1) (unsafeLinearModify marr f)
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+1) (unsafeLinearModify marr f)
 {-# INLINE forPrimM #-}
 
 
@@ -904,7 +959,7 @@
 iforLinearPrimM_ ::
      (Manifest r e, Index ix, PrimMonad m) => MArray (PrimState m) r ix e -> (Int -> e -> m ()) -> m ()
 iforLinearPrimM_ marr f =
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\i -> unsafeLinearRead marr i >>= f i)
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\i -> unsafeLinearRead marr i >>= f i)
 {-# INLINE iforLinearPrimM_ #-}
 
 -- | Sequentially loop over a mutable array while modifying each element with an index aware action.
@@ -913,7 +968,7 @@
 iforLinearPrimM ::
      (Manifest r e, Index ix, PrimMonad m) => MArray (PrimState m) r ix e -> (Int -> e -> m e) -> m ()
 iforLinearPrimM marr f =
-  loopM_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\i -> unsafeLinearModify marr (f i) i)
+  loopA_ 0 (< totalElem (sizeOfMArray marr)) (+ 1) (\i -> unsafeLinearModify marr (f i) i)
 {-# INLINE iforLinearPrimM #-}
 
 
@@ -943,7 +998,7 @@
   -> m ()
 ifor2PrimM_ m1 m2 f = do
   let sz = liftIndex2 min (unSz (sizeOfMArray m1)) (unSz (sizeOfMArray m2))
-  iterM_ zeroIndex sz oneIndex (<) $ \ix -> do
+  iterA_ zeroIndex sz oneIndex (<) $ \ix -> do
     e1 <- unsafeRead m1 ix
     e2 <- unsafeRead m2 ix
     f ix e1 e2
@@ -1128,7 +1183,7 @@
         MArray (PrimState m) r ix e -> ix -> m e
 readM marr ix =
   read marr ix >>= \case
-    Just e -> pure e
+    Just e  -> pure e
     Nothing -> throwM $ IndexOutOfBoundsException (sizeOfMArray marr) ix
 {-# INLINE readM #-}
 
@@ -1304,7 +1359,7 @@
   let sz1 = sizeOfMArray m1
       sz2 = sizeOfMArray m2
       sz = liftIndex2 min (unSz sz1) (unSz sz2)
-  iterM_ startIx sz oneIndex (<) $ \ix -> do
+  iterA_ startIx sz oneIndex (<) $ \ix -> do
     let i1 = toLinearIndex sz1 ix
         i2 = toLinearIndex sz2 ix
     e1 <- unsafeLinearRead m1 i1
diff --git a/src/Data/Massiv/Array/Mutable/Algorithms.hs b/src/Data/Massiv/Array/Mutable/Algorithms.hs
--- a/src/Data/Massiv/Array/Mutable/Algorithms.hs
+++ b/src/Data/Massiv/Array/Mutable/Algorithms.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Algorithms
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Mutable/Atomic.hs b/src/Data/Massiv/Array/Mutable/Atomic.hs
--- a/src/Data/Massiv/Array/Mutable/Atomic.hs
+++ b/src/Data/Massiv/Array/Mutable/Atomic.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Atomic
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Mutable/Internal.hs b/src/Data/Massiv/Array/Mutable/Internal.hs
--- a/src/Data/Massiv/Array/Mutable/Internal.hs
+++ b/src/Data/Massiv/Array/Mutable/Internal.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE ExplicitForAll #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Numeric.hs b/src/Data/Massiv/Array/Numeric.hs
--- a/src/Data/Massiv/Array/Numeric.hs
+++ b/src/Data/Massiv/Array/Numeric.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Numeric
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -147,9 +147,9 @@
   -> Array r ix e
   -> m (Array r ix e)
 applyExactSize2M f a1 a2
-  | size a1 == size a2 = pure $ f a1 a2
-  | isZeroSz sz1 && isZeroSz sz2 = pure $ unsafeResize zeroSz a1
-  | otherwise = throwM $ SizeMismatchException sz1 sz2
+  | size a1 == size a2 = pure $! f a1 a2
+  | isZeroSz sz1 && isZeroSz sz2 = pure $! unsafeResize zeroSz a1
+  | otherwise = throwM $! SizeMismatchException sz1 sz2
   where
     !sz1 = size a1
     !sz2 = size a2
@@ -180,7 +180,7 @@
 --   [ 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40 ]
 --
 -- @since 0.5.6
-(!+!) :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e -> Array r ix e
+(!+!) :: (HasCallStack, Index ix, Numeric r e) => Array r ix e -> Array r ix e -> Array r ix e
 (!+!) a1 a2 = throwEither (a1 .+. a2)
 {-# INLINE (!+!) #-}
 
@@ -364,7 +364,7 @@
     withScheduler comp $ \scheduler ->
       splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> liftIO $ do
         let n = SafeSz chunkLength
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork scheduler $
           pure $! unsafeDotProduct (unsafeLinearSlice start n v1) (unsafeLinearSlice start n v2)
         when (slackStart < totalLength) $ do
@@ -398,7 +398,7 @@
     withScheduler (getComp v) $ \scheduler ->
       splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> liftIO $ do
         let n = SafeSz chunkLength
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork scheduler $ pure $! powerSumArray (unsafeLinearSlice start n v) p
         when (slackStart < totalLength) $ do
           let k = SafeSz (totalLength - slackStart)
@@ -718,7 +718,7 @@
   | isEmpty arr1 || isEmpty arr2 = pure $ setComp comp empty
   | otherwise =
     pure $
-    DArray comp (SafeSz (m1 :. n2)) $ \(i :. j) ->
+    makeArray comp (SafeSz (m1 :. n2)) $ \(i :. j) ->
       unsafeDotProduct (unsafeLinearSlice (i * n1) n arr1) (unsafeLinearSlice (j * n1) n arr2)
   where
     comp = getComp arr1 <> getComp arr2
@@ -744,7 +744,7 @@
 -- @since 0.3.6
 identityMatrix :: Num e => Sz1 -> Matrix DL e
 identityMatrix (Sz n) =
-  makeLoadArrayS (Sz2 n n) 0 $ \ w -> loopM_ 0 (< n) (+1) $ \ i -> w (i :. i) 1
+  makeLoadArrayS (Sz2 n n) 0 $ \ w -> loopA_ 0 (< n) (+1) $ \ i -> w (i :. i) 1
 {-# INLINE identityMatrix #-}
 
 -- | Create a lower triangular (L in LU decomposition) matrix of size @NxN@
@@ -762,15 +762,16 @@
 --   ]
 --
 -- @since 0.5.2
-lowerTriangular :: Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
-lowerTriangular comp (Sz1 n) f =
-  let sz = Sz2 n n
-   in unsafeMakeLoadArrayAdjusted comp sz (Just 0) $ \scheduler wr ->
-        forM_ (0 ..: n) $ \i ->
-          scheduleWork scheduler $
-          forM_ (0 ... i) $ \j ->
-            let ix = i :. j
-             in wr (toLinearIndex sz ix) (f ix)
+lowerTriangular :: forall e. Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
+lowerTriangular comp (Sz1 n) f = DLArray comp (SafeSz (n :. n)) load
+  where
+    load :: Loader e
+    load scheduler startAt uWrite uSet = do
+      forM_ (0 ..: n) $ \i -> do
+        let !k = startAt + i * n
+        scheduleWork_ scheduler $ do
+          forM_ (0 ... i) $ \j -> uWrite (k + j) (f (i :. j))
+          uSet (k + i + 1) (Sz (n - i - 1)) 0
 {-# INLINE lowerTriangular #-}
 
 -- | Create an upper triangular (U in LU decomposition) matrix of size @NxN@
@@ -788,15 +789,16 @@
 --   ]
 --
 -- @since 0.5.2
-upperTriangular :: Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
-upperTriangular comp (Sz1 n) f =
-  let sz = Sz2 n n
-   in unsafeMakeLoadArrayAdjusted comp sz (Just 0) $ \scheduler wr ->
-        forM_ (0 ..: n) $ \i ->
-          scheduleWork scheduler $
-          forM_ (i ..: n) $ \j ->
-            let ix = i :. j
-             in wr (toLinearIndex sz ix) (f ix)
+upperTriangular :: forall e. Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
+upperTriangular comp (Sz1 n) f = DLArray comp (SafeSz (n :. n)) load
+  where
+    load :: Loader e
+    load scheduler startAt uWrite uSet = do
+      forM_ (0 ..: n) $ \i -> do
+        let !k = startAt + i * n
+        scheduleWork_ scheduler $ do
+          uSet k (SafeSz i) 0
+          forM_ (i ..: n) $ \j -> uWrite (k + j) (f (i :. j))
 {-# INLINE upperTriangular #-}
 
 -- | Negate each element of the array
diff --git a/src/Data/Massiv/Array/Numeric/Integral.hs b/src/Data/Massiv/Array/Numeric/Integral.hs
--- a/src/Data/Massiv/Array/Numeric/Integral.hs
+++ b/src/Data/Massiv/Array/Numeric/Integral.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Numeric.Integral
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -7,7 +7,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Construct
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -29,13 +29,10 @@
   , iiterateN
     -- *** Unfolding
   , unfoldlS_
-  -- , unfoldlS
   , iunfoldlS_
-  --, iunfoldlS
   , unfoldrS_
-  --, unfoldrS
   , iunfoldrS_
-  --, iunfoldrS
+  , makeSplitSeedArray
     -- *** Random
   , uniformArray
   , uniformRangeArray
@@ -134,16 +131,10 @@
   => Sz ix
   -> (ix -> f e)
   -> f (Array r ix e)
-makeArrayA !sz f =
-  let n = totalElem sz
-      go !i
-        | i < n =
-          liftA2
-            (\e (STA st) -> STA (\ma -> unsafeLinearWrite ma i e >> st ma))
-            (f (fromLinearIndex sz i))
-            (go (i + 1))
-        | otherwise = pure (STA (unsafeFreeze Seq))
-   in runSTA sz <$> go 0
+makeArrayA sz@(Sz n) f =
+  fmap (runSTA sz) $
+  iterF zeroIndex n oneIndex (<) (pure (STA (unsafeFreeze Seq))) $ \ix g ->
+    liftA2 (\e (STA st) -> STA (\ma -> unsafeWrite ma ix e >> st ma)) (f ix) g
 {-# INLINE makeArrayA  #-}
 
 -- | Same as `makeArrayA`, but with linear index.
@@ -155,12 +146,9 @@
   -> (Int -> f e)
   -> f (Array r ix e)
 makeArrayLinearA !sz f =
-  let n = totalElem sz
-      go !i
-        | i < n =
-          liftA2 (\e (STA st) -> STA (\ma -> unsafeLinearWrite ma i e >> st ma)) (f i) (go (i + 1))
-        | otherwise = pure (STA (unsafeFreeze Seq))
-   in runSTA sz <$> go 0
+  fmap (runSTA sz) $
+  loopF 0 (< totalElem sz) (+ 1) (pure (STA (unsafeFreeze Seq))) $ \i ->
+    liftA2 (\e (STA st) -> STA (\ma -> unsafeLinearWrite ma i e >> st ma)) (f i)
 {-# INLINE makeArrayLinearA  #-}
 
 
@@ -234,16 +222,16 @@
     load :: Loader e
     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
+      iterTargetM defRowMajor startAt sz zeroIndex oneStride acc0 $ \ !i !ix !acc ->
+        case f acc ix of
+          (e, !acc') -> acc' <$ dlWrite i e
     {-# INLINE load #-}
 {-# INLINE iunfoldrS_ #-}
 
 
 -- | Unfold sequentially from the end. There is no way to save the accumulator after
 -- unfolding is done, since resulting array is delayed, but it's possible to use
--- `Data.Massiv.Array.Mutable.unfoldlPrimM` to achive such effect.
+-- `Data.Massiv.Array.Mutable.unfoldlPrimM` to achieve such effect.
 --
 -- @since 0.3.0
 unfoldlS_ :: Index ix => Sz ix -> (a -> (a, e)) -> a -> Array DL ix e
@@ -280,6 +268,9 @@
 -- to parallelize the random value generation, but also guarantee that it will be
 -- deterministic, granted none of the arguments have changed.
 --
+-- __Note__: Starting with massiv-1.1.0 this function will be deprecated in
+-- favor of a more general `genSplitArray`
+--
 -- ==== __Examples__
 --
 -- >>> import Data.Massiv.Array
@@ -332,12 +323,50 @@
             scheduleWork_ scheduler $
               void $ loopM start (< start + chunkLength) (+ 1) genI0 writeRandom
             pure genI1
-        when (slackStartAt < totalLength + startAt) $
+        when (slackStart < totalLength) $
           scheduleWork_ scheduler $
           void $ loopM slackStartAt (< totalLength + startAt) (+ 1) genForSlack writeRandom
 {-# INLINE randomArray #-}
 
 
+
+
+-- | Create a delayed array with an initial seed and a splitting function. It is
+-- somewhat similar to `iunfoldlS_` function, but it is capable of parallelizing
+-- computation and iterating over the array accoriding to the supplied
+-- `Iterator`. Upon parallelization every job will get the second part of the
+-- result produced by the split function, while the first part will be used for
+-- subsequent splits. This function is similar to
+-- `Data.Massiv.Array.Manifest.generateSplitSeedArray`
+--
+-- @since 1.0.2
+makeSplitSeedArray ::
+     forall ix e g it. (Iterator it, Index ix)
+  => it -- ^ Iterator
+  -> g -- ^ Initial seed
+  -> (g -> (g, g))
+     -- ^ A function that can split a seed into two independent seeds. It will
+     -- be called the same number of times as the number of jobs that will get
+     -- scheduled during parallelization. Eg. only once for the sequential case.
+  -> Comp -- ^ Computation strategy.
+  -> Sz ix -- ^ Resulting size of the array.
+  -> (Ix1 -> ix -> g -> (e, g))
+     -- ^ A function that produces a value and the next seed. It takes both
+     -- versions of the index, in linear and in multi-dimensional forms, as well as
+     -- the current seeding value.
+  -> Array DL ix e
+makeSplitSeedArray it seed splitSeed comp sz genFunc =
+  DLArray {dlComp = comp, dlSize = sz, dlLoad = load}
+  where
+    load :: Loader e
+    load scheduler startAt writeAt _ =
+      iterTargetFullAccST_ it scheduler startAt sz seed (pure . splitSeed) $ \ i ix g ->
+        case genFunc (i - startAt) ix g of
+          (x, g') -> g' <$ writeAt i x
+    {-# INLINE load #-}
+{-# INLINE makeSplitSeedArray #-}
+
+
 -- | Generate a random array where all elements are sampled from a uniform distribution.
 --
 -- @since 1.0.0
@@ -348,6 +377,7 @@
   -> Sz ix -- ^ Resulting size of the array.
   -> Array DL ix e
 uniformArray gen = randomArray gen split uniform
+{-# INLINE uniformArray #-}
 
 -- | Same as `uniformArray`, but will generate values in a supplied range.
 --
@@ -360,6 +390,7 @@
   -> Sz ix -- ^ Resulting size of the array.
   -> Array DL ix e
 uniformRangeArray gen r = randomArray gen split (uniformR r)
+{-# INLINE uniformRangeArray #-}
 
 
 -- | Similar to `randomArray` but performs generation sequentially, which means it doesn't
@@ -516,7 +547,7 @@
      forall ix m. (Index ix, MonadThrow m)
   => Comp -- ^ Computation strategy
   -> ix -- ^ Start
-  -> ix -- ^ Step (Can't have zeros)
+  -> ix -- ^ Step. Negative and positive values are ok, but can't have zeros
   -> ix -- ^ End
   -> m (Array D ix ix)
 rangeStepM comp !from !step !to
diff --git a/src/Data/Massiv/Array/Ops/Fold.hs b/src/Data/Massiv/Array/Ops/Fold.hs
--- a/src/Data/Massiv/Array/Ops/Fold.hs
+++ b/src/Data/Massiv/Array/Ops/Fold.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Fold
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -6,7 +6,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Fold.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -40,6 +40,7 @@
   , foldrP
   , ifoldlP
   , ifoldrP
+  , foldlIO
   , ifoldlIO
   , ifoldrIO
   , splitReduce
@@ -49,11 +50,12 @@
   ) where
 
 import Control.Monad (void, when)
+import Control.Monad.Primitive
 import Control.Scheduler
 import qualified Data.Foldable as F
 import Data.Functor.Identity (runIdentity)
 import Data.Massiv.Core.Common
-import Prelude hiding (foldl, foldr, any)
+import Prelude hiding (any, foldl, foldr)
 import System.IO.Unsafe (unsafePerformIO)
 
 
@@ -86,7 +88,14 @@
 --
 -- @since 0.1.0
 foldlM :: (Index ix, Source r e, Monad m) => (a -> e -> m a) -> a -> Array r ix e -> m a
-foldlM f = ifoldlM (\ a _ b -> f a b)
+foldlM f acc arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterM zeroIndex (unSz sz) (pureIndex 1) (<) acc $ \ !ix !a -> f a (gix ix)
+    PrefIndexLinear gi ->
+      loopM 0 (< totalElem sz) (+ 1) acc $ \ !i !a -> f a (gi i)
+  where
+    sz = size arr
 {-# INLINE foldlM #-}
 
 
@@ -94,7 +103,7 @@
 --
 -- @since 0.1.0
 foldlM_ :: (Index ix, Source r e, Monad m) => (a -> e -> m a) -> a -> Array r ix e -> m ()
-foldlM_ f = ifoldlM_ (\ a _ b -> f a b)
+foldlM_ f acc = void . foldlM f acc
 {-# INLINE foldlM_ #-}
 
 
@@ -103,7 +112,11 @@
 -- @since 0.1.0
 ifoldlM :: (Index ix, Source r e, Monad m) => (a -> ix -> e -> m a) -> a -> Array r ix e -> m a
 ifoldlM f !acc !arr =
-  iterM zeroIndex (unSz (size arr)) (pureIndex 1) (<) acc $ \ !ix !a -> f a ix (unsafeIndex arr ix)
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterM zeroIndex (unSz (size arr)) (pureIndex 1) (<) acc $ \ !ix !a -> f a ix (gix ix)
+    PrefIndexLinear gi ->
+      iterTargetM defRowMajor 0 (size arr) zeroIndex oneStride acc $ \i ix !a -> f a ix (gi i)
 {-# INLINE ifoldlM #-}
 
 
@@ -119,7 +132,14 @@
 --
 -- @since 0.1.0
 foldrM :: (Index ix, Source r e, Monad m) => (e -> a -> m a) -> a -> Array r ix e -> m a
-foldrM f = ifoldrM (\_ e a -> f e a)
+foldrM f acc arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterM (liftIndex (subtract 1) (unSz sz)) zeroIndex (pureIndex (-1)) (>=) acc (f . gix)
+    PrefIndexLinear gi ->
+      loopM (totalElem sz - 1) (>= 0) (subtract 1) acc (f . gi)
+  where
+    !sz = size arr
 {-# INLINE foldrM #-}
 
 
@@ -136,8 +156,8 @@
 -- @since 0.1.0
 ifoldrM :: (Index ix, Source r e, Monad m) => (ix -> e -> a -> m a) -> a -> Array r ix e -> m a
 ifoldrM f !acc !arr =
-  iterM (liftIndex (subtract 1) (unSz (size arr))) zeroIndex (pureIndex (-1)) (>=) acc $ \ !ix !acc0 ->
-    f ix (unsafeIndex arr ix) acc0
+  iterM (liftIndex (subtract 1) (unSz (size arr))) zeroIndex (pureIndex (-1)) (>=) acc $ \ !ix ->
+    f ix (unsafeIndex arr ix)
 {-# INLINE ifoldrM #-}
 
 
@@ -157,7 +177,7 @@
 lazyFoldlS f initAcc arr = go initAcc 0
   where
     len = totalElem (size arr)
-    go acc k
+    go acc !k
       | k < len = go (f acc (unsafeLinearIndex arr k)) (k + 1)
       | otherwise = acc
 {-# INLINE lazyFoldlS #-}
@@ -175,7 +195,7 @@
 --
 -- @since 0.1.0
 foldlS :: (Index ix, Source r e) => (a -> e -> a) -> a -> Array r ix e -> a
-foldlS f = ifoldlS (\ a _ e -> f a e)
+foldlS f acc = runIdentity . foldlM (\ a e -> pure $! f a e) acc
 {-# INLINE foldlS #-}
 
 
@@ -184,7 +204,7 @@
 -- @since 0.1.0
 ifoldlS :: (Index ix, Source r e)
         => (a -> ix -> e -> a) -> a -> Array r ix e -> a
-ifoldlS f acc = runIdentity . ifoldlM (\ a ix e -> return $ f a ix e) acc
+ifoldlS f acc = runIdentity . ifoldlM (\ a ix e -> pure $! f a ix e) acc
 {-# INLINE ifoldlS #-}
 
 
@@ -192,7 +212,7 @@
 --
 -- @since 0.1.0
 foldrS :: (Index ix, Source r e) => (e -> a -> a) -> a -> Array r ix e -> a
-foldrS f = ifoldrS (\_ e a -> f e a)
+foldrS f acc = runIdentity . foldrM (\ e a -> pure $! f e a) acc
 {-# INLINE foldrS #-}
 
 
@@ -200,7 +220,7 @@
 --
 -- @since 0.1.0
 ifoldrS :: (Index ix, Source r e) => (ix -> e -> a -> a) -> a -> Array r ix e -> a
-ifoldrS f acc = runIdentity . ifoldrM (\ ix e a -> return $ f ix e a) acc
+ifoldrS f acc = runIdentity . ifoldrM (\ ix e a -> pure $! f ix e a) acc
 {-# INLINE ifoldrS #-}
 
 
@@ -213,7 +233,7 @@
     !k = totalElem (size arr)
     go !i
       | i == k = n
-      | otherwise = let !v = unsafeLinearIndex arr i in v `c` go (i + 1)
+      | otherwise = let v = unsafeLinearIndex arr i in v `c` go (i + 1)
 {-# INLINE [0] foldrFB #-}
 
 
@@ -242,7 +262,8 @@
        -> (b -> a -> b) -- ^ Chunk results folding function @f@.
        -> b -- ^ Accumulator for results of chunks folding.
        -> Array r ix e -> m b
-foldlP f fAcc g gAcc = liftIO . ifoldlP (\ x _ -> f x) fAcc g gAcc
+foldlP f fAcc g gAcc =
+  liftIO . foldlIO (\acc -> pure . f acc) fAcc (\acc -> pure . g acc) gAcc
 {-# INLINE foldlP #-}
 
 -- | /O(n)/ - Left fold with an index aware function, computed in parallel. Just
@@ -253,7 +274,7 @@
 ifoldlP :: (MonadIO m, Index ix, Source r e) =>
            (a -> ix -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> m b
 ifoldlP f fAcc g gAcc =
-  liftIO . ifoldlIO (\acc ix -> return . f acc ix) fAcc (\acc -> return . g acc) gAcc
+  liftIO . ifoldlIO (\acc ix -> pure . f acc ix) fAcc (\acc -> pure . g acc) gAcc
 {-# INLINE ifoldlP #-}
 
 
@@ -296,16 +317,48 @@
 
 -- | This folding function breaks referential transparency on some functions
 -- @f@, therefore it is kept here for internal use only.
-foldlInternal :: (Index ix, Source r e) => (a -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> b
+foldlInternal ::
+     (Index ix, Source r e) => (a -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> b
 foldlInternal g initAcc f resAcc = unsafePerformIO . foldlP g initAcc f resAcc
 {-# INLINE foldlInternal #-}
 
 
-ifoldlInternal :: (Index ix, Source r e) => (a -> ix -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> b
+ifoldlInternal ::
+     (Index ix, Source r e) => (a -> ix -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> b
 ifoldlInternal g initAcc f resAcc = unsafePerformIO . ifoldlP g initAcc f resAcc
 {-# INLINE ifoldlInternal #-}
 
 
+-- | Similar to `foldlP`, except that folding functions themselves do live in IO
+--
+-- @since 0.1.0
+foldlIO ::
+     (MonadUnliftIO m, Index ix, Source r e)
+  => (a -> e -> m a) -- ^ Index aware folding IO action
+  -> a -- ^ Accumulator
+  -> (b -> a -> m b) -- ^ Folding action that is applied to the results of a parallel fold
+  -> b -- ^ Accumulator for chunks folding
+  -> Array r ix e
+  -> m b
+foldlIO f !initAcc g !tAcc !arr
+  | getComp arr == Seq = foldlM f initAcc arr >>= g tAcc
+  | otherwise = do
+    let splitAcc _ = pure (initAcc, initAcc)
+        !sz = size arr
+    results <-
+      withScheduler (getComp arr) $ \scheduler ->
+        withRunInIO $ \run ->
+          stToPrim $
+          case unsafePrefIndex arr of
+            PrefIndex gix ->
+              iterFullAccST defRowMajor scheduler zeroIndex sz initAcc splitAcc $ \ !ix !acc ->
+                ioToPrim (run (f acc (gix ix)))
+            PrefIndexLinear gi ->
+              iterFullAccST defRowMajor scheduler 0 (toLinearSz sz) initAcc splitAcc $ \ !i !acc ->
+                ioToPrim (run (f acc (gi i)))
+    F.foldlM g tAcc results
+{-# INLINE foldlIO #-}
+
 -- | Similar to `ifoldlP`, except that folding functions themselves do live in IO
 --
 -- @since 0.1.0
@@ -320,21 +373,20 @@
 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 -> do
-          withRunInIO $ \run -> do
-            splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-              loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
-                scheduleWork scheduler $ run $
-                iterLinearM sz start (start + chunkLength) 1 (<) initAcc $ \ !i ix !acc ->
-                  f acc ix (unsafeLinearIndex arr i)
-              when (slackStart < totalLength) $
-                scheduleWork scheduler $ run $
-                iterLinearM sz slackStart totalLength 1 (<) initAcc $ \ !i ix !acc ->
-                  f acc ix (unsafeLinearIndex arr i)
-      F.foldlM g tAcc results
+    let !sz = size arr
+        splitAcc _ = pure (initAcc, initAcc)
+    results <-
+      withScheduler (getComp arr) $ \scheduler ->
+        withRunInIO $ \run ->
+          stToPrim $
+          case unsafePrefIndex arr of
+            PrefIndex gix ->
+              iterFullAccST defRowMajor scheduler zeroIndex sz initAcc splitAcc $ \ !ix !acc ->
+                ioToPrim (run (f acc ix (gix ix)))
+            PrefIndexLinear gi ->
+              iterTargetFullAccST defRowMajor scheduler 0 sz initAcc splitAcc $ \ !i !ix !acc ->
+                ioToPrim (run (f acc ix (gi i)))
+    F.foldlM g tAcc results
 {-# INLINE ifoldlIO #-}
 
 -- | Slice an array into linear row-major vector chunks and apply an action to each of
@@ -356,7 +408,7 @@
     withScheduler (getComp arr) $ \scheduler -> do
       withRunInIO $ \run -> do
         splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-          loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+          loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
             scheduleWork scheduler $ run $ f scheduler $
               unsafeLinearSlice start (SafeSz chunkLength) arr
           when (slackStart < totalLength) $
@@ -383,12 +435,12 @@
           splitLinearly (numWorkers scheduler) totalLength $ \ chunkLength slackStart -> do
             when (slackStart < totalLength) $
               scheduleWork scheduler $ run $
-              iterLinearM sz (totalLength - 1) slackStart (-1) (>=) initAcc $ \ !i ix !acc ->
-                f ix (unsafeLinearIndex arr i) acc
-            loopM_ slackStart (> 0) (subtract chunkLength) $ \ !start ->
+              iterLinearM sz (totalLength - 1) slackStart (-1) (>=) initAcc $ \ !i ix ->
+                f ix (unsafeLinearIndex arr i)
+            loopA_ slackStart (> 0) (subtract chunkLength) $ \ !start ->
               scheduleWork scheduler $ run $
-                iterLinearM sz (start - 1) (start - chunkLength) (-1) (>=) initAcc $ \ !i ix !acc ->
-                  f ix (unsafeLinearIndex arr i) acc
+                iterLinearM sz (start - 1) (start - chunkLength) (-1) (>=) initAcc $ \ !i ix ->
+                  f ix (unsafeLinearIndex arr i)
     F.foldlM (flip g) tAcc results
 {-# INLINE ifoldrIO #-}
 
@@ -448,6 +500,9 @@
 
 -- | Parallelizable implementation of `any` with unrolling
 anyPu :: (Index ix, Source r e) => (e -> Bool) -> Array r ix e -> IO Bool
+-- TODO: switch to splitReduce
+-- anyPu f arr =
+--   splitReduce anySu (\r acc -> pure (r || acc)) False
 anyPu f arr = do
   let !sz = size arr
       !totalLength = totalElem sz
@@ -455,7 +510,7 @@
     withScheduler (getComp arr) $ \scheduler -> do
       batch <- getCurrentBatch scheduler
       splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork scheduler $ anySliceSuM batch start (Sz (start + chunkLength)) f arr
         when (slackStart < totalLength) $
           scheduleWork scheduler $ anySliceSuM batch slackStart (Sz totalLength) f arr
@@ -471,5 +526,5 @@
 any f arr =
   case getComp arr of
     Seq -> anySu f arr
-    _ -> unsafePerformIO $ anyPu f arr
+    _   -> unsafePerformIO $ anyPu f arr
 {-# INLINE any #-}
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE MonoLocalBinds #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Map
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -73,6 +73,8 @@
   , izipWith3A
   ) where
 
+import Data.Traversable (traverse)
+import Data.Massiv.Array.Manifest.List
 import Control.Monad (void)
 import Control.Monad.Primitive
 import Control.Scheduler
@@ -92,7 +94,7 @@
 --
 -- @since 0.1.0
 map :: (Index ix, Source r e') => (e' -> e) -> Array r ix e' -> Array D ix e
-map f = imap (const f)
+map f = fmap f . delay
 {-# INLINE map #-}
 
 
@@ -164,7 +166,17 @@
 -- source arrays in case their dimensions do not match.
 zipWith :: (Index ix, Source r1 e1, Source r2 e2)
         => (e1 -> e2 -> e) -> Array r1 ix e1 -> Array r2 ix e2 -> Array D ix e
-zipWith f = izipWith (\ _ e1 e2 -> f e1 e2)
+zipWith f arr1 arr2 = DArray comp sz prefIndex
+  where
+    sz = SafeSz (liftIndex2 min (coerce (size arr1)) (coerce (size arr2)))
+    comp = getComp arr1 <> getComp arr2
+    prefIndex = PrefIndex (\ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
+      -- Somehow checking for size equality destroys performance
+      --  | PrefIndexLinear gi1 <- unsafePrefIndex arr1,
+      --    PrefIndexLinear gi2 <- unsafePrefIndex arr2,
+      --    size arr1 == size arr2 =
+      --      PrefIndexLinear (\i -> f (gi1 i) (gi2 i))
+      --  | otherwise = PrefIndex (\ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
 {-# INLINE zipWith #-}
 
 
@@ -174,15 +186,30 @@
 izipWith f arr1 arr2 =
   DArray
     (getComp arr1 <> getComp arr2)
-    (SafeSz (liftIndex2 min (coerce (size arr1)) (coerce (size arr2)))) $ \ !ix ->
-    f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix)
+    (SafeSz (liftIndex2 min (coerce (size arr1)) (coerce (size arr2))))
+    (PrefIndex (\ix -> f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix)))
 {-# INLINE izipWith #-}
 
 
 -- | Just like `zipWith`, except zip three arrays with a function.
-zipWith3 :: (Index ix, Source r1 e1, Source r2 e2, Source r3 e3)
-         => (e1 -> e2 -> e3 -> e) -> Array r1 ix e1 -> Array r2 ix e2 -> Array r3 ix e3 -> Array D ix e
-zipWith3 f = izipWith3 (\ _ e1 e2 e3 -> f e1 e2 e3)
+zipWith3 ::
+     (Index ix, Source r1 e1, Source r2 e2, Source r3 e3)
+  => (e1 -> e2 -> e3 -> e)
+  -> Array r1 ix e1
+  -> Array r2 ix e2
+  -> Array r3 ix e3
+  -> Array D ix e
+zipWith3 f arr1 arr2 arr3 = izipWith3 (\_ e1 e2 e3 -> f e1 e2 e3) arr1 arr2 arr3
+  -- See note on zipWith
+  --  | sz1 == size arr2 && sz1 == size arr3
+  --  , PrefIndexLinear gi1 <- unsafePrefIndex arr1
+  --  , PrefIndexLinear gi2 <- unsafePrefIndex arr2
+  --  , PrefIndexLinear gi3 <- unsafePrefIndex arr3 =
+  --    makeArrayLinear comp sz1 (\ !i -> f (gi1 i) (gi2 i) (gi3 i))
+  --  | otherwise = izipWith3 (\_ e1 e2 e3 -> f e1 e2 e3) arr1 arr2 arr3
+  -- where
+  --   comp = getComp arr1 <> getComp arr2 <> getComp arr3
+  --   sz1 = size arr1
 {-# INLINE zipWith3 #-}
 
 
@@ -201,8 +228,8 @@
        (liftIndex2
           min
           (liftIndex2 min (coerce (size arr1)) (coerce (size arr2)))
-          (coerce (size arr3)))) $ \ !ix ->
-    f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix)
+          (coerce (size arr3))))
+    (PrefIndex $ \ !ix -> f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix))
 {-# INLINE izipWith3 #-}
 
 
@@ -218,7 +245,19 @@
   -> Array r3 ix e3
   -> Array r4 ix e4
   -> Array D ix e
-zipWith4 f = izipWith4 (\ _ e1 e2 e3 e4 -> f e1 e2 e3 e4)
+zipWith4 f arr1 arr2 arr3 arr4 =
+  izipWith4 (\ _ e1 e2 e3 e4 -> f e1 e2 e3 e4) arr1 arr2 arr3 arr4
+  -- See note on zipWith
+  --  | sz1 == size arr2 && sz1 == size arr3 && sz1 == size arr4
+  --  , PrefIndexLinear gi1 <- unsafePrefIndex arr1
+  --  , PrefIndexLinear gi2 <- unsafePrefIndex arr2
+  --  , PrefIndexLinear gi3 <- unsafePrefIndex arr3
+  --  , PrefIndexLinear gi4 <- unsafePrefIndex arr4 =
+  --    makeArrayLinear comp sz1 (\ !i -> f (gi1 i) (gi2 i) (gi3 i) (gi4 i))
+  --  | otherwise = izipWith4 (\ _ e1 e2 e3 e4 -> f e1 e2 e3 e4) arr1 arr2 arr3 arr4
+  --  where
+  --    comp = getComp arr1 <> getComp arr2 <> getComp arr3 <> getComp arr4
+  --    sz1 = size arr1
 {-# INLINE zipWith4 #-}
 
 
@@ -234,7 +273,7 @@
   -> Array r4 ix e4
   -> Array D ix e
 izipWith4 f arr1 arr2 arr3 arr4 =
-  DArray
+  makeArray
     (getComp arr1 <> getComp arr2 <> getComp arr3 <> getComp arr4)
     (SafeSz
        (liftIndex2
@@ -243,8 +282,9 @@
              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)
+          (coerce (size arr4))))
+    (\ !ix ->
+       f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix) (unsafeIndex arr4 ix))
 {-# INLINE izipWith4 #-}
 
 
@@ -253,12 +293,19 @@
 --
 -- @since 0.3.0
 zipWithA ::
-     (Source r1 e1, Source r2 e2, Applicative f, Manifest r e, Index ix)
-  => (e1 -> e2 -> f e)
-  -> Array r1 ix e1
-  -> Array r2 ix e2
-  -> f (Array r ix e)
-zipWithA f = izipWithA (const f)
+  (Source r1 e1, Source r2 e2, Applicative f, Manifest r e, Index ix) =>
+  (e1 -> e2 -> f e) ->
+  Array r1 ix e1 ->
+  Array r2 ix e2 ->
+  f (Array r ix e)
+zipWithA f arr1 arr2
+  | sz1 == size arr2
+  , PrefIndexLinear gi1 <- unsafePrefIndex arr1
+  , PrefIndexLinear gi2 <- unsafePrefIndex arr2 =
+    setComp (getComp arr1 <> getComp arr2) <$> makeArrayLinearA sz1 (\ !i -> f (gi1 i) (gi2 i))
+  | otherwise = izipWithA (const f) arr1 arr2
+  where
+    !sz1 = size arr1
 {-# INLINE zipWithA #-}
 
 -- | Similar to `zipWith`, except does it sequentiall and using the `Applicative`. Note that
@@ -326,7 +373,8 @@
   => (a -> f e)
   -> Array r' ix a
   -> f (Array r ix e)
-traverseA f arr = makeArrayLinearA (size arr) (f . unsafeLinearIndex arr)
+traverseA f arr =
+  unsafeResize (size arr) . fromList (getComp arr) <$> traverse f (toList arr)
 {-# INLINE traverseA #-}
 
 -- | Traverse sequentially over a source array, while discarding the result.
@@ -338,7 +386,12 @@
   => (e -> f a)
   -> Array r ix e
   -> f ()
-traverseA_ f arr = loopA_ 0 (< totalElem (size arr)) (+ 1) (f . unsafeLinearIndex arr)
+traverseA_ f arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix -> iterA_ zeroIndex (unSz sz) oneIndex (<) (f . gix)
+    PrefIndexLinear gi -> loopA_ 0 (< totalElem sz) (+ 1) (f . gi)
+  where
+    sz = size arr
 {-# INLINE traverseA_ #-}
 
 -- | Sequence actions in a source array.
@@ -388,7 +441,11 @@
   -> Array r ix a
   -> f ()
 itraverseA_ f arr =
-  loopA_ 0 (< totalElem sz) (+ 1) (\ !i -> f (fromLinearIndex sz i) (unsafeLinearIndex arr i))
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterA_ zeroIndex (unSz sz) oneIndex (<) (\ !ix -> f ix (gix ix))
+    PrefIndexLinear gi ->
+      iterTargetA_ defRowMajor 0 sz zeroIndex oneStride $ \i ix -> f ix (gi i)
   where
     sz = size arr
 {-# INLINE itraverseA_ #-}
@@ -403,7 +460,17 @@
   => (a -> m b)
   -> Array r' ix a
   -> m (Array r ix b)
-traversePrim f = itraversePrim (const f)
+traversePrim f arr = do
+  let sz = size arr
+  marr <- unsafeNew sz
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterTargetA_ defRowMajor 0 sz zeroIndex oneStride $ \i ix ->
+        f (gix ix) >>= unsafeLinearWrite marr i
+    PrefIndexLinear gi ->
+      loopA_ 0 (< totalElem sz) (+ 1) $ \i ->
+        f (gi i) >>= unsafeLinearWrite marr i
+  unsafeFreeze (getComp arr) marr
 {-# INLINE traversePrim #-}
 
 -- | Same as `traversePrim`, but traverse with index aware action.
@@ -415,13 +482,17 @@
   => (ix -> a -> m b)
   -> Array r' ix a
   -> m (Array r ix b)
-itraversePrim f arr =
-  setComp (getComp arr) <$>
-  generateArrayLinearS
-    (size arr)
-    (\ !i ->
-       let ix = fromLinearIndex (size arr) i
-        in f ix (unsafeLinearIndex arr i))
+itraversePrim f arr = do
+  let sz = size arr
+  marr <- unsafeNew sz
+  case unsafePrefIndex arr of
+    PrefIndex gix ->
+      iterTargetA_ defRowMajor 0 sz zeroIndex oneStride $ \i ix ->
+        f ix (gix ix) >>= unsafeLinearWrite marr i
+    PrefIndexLinear gi ->
+      iterTargetA_ defRowMajor 0 sz zeroIndex oneStride $ \i ix ->
+        f ix (gi i) >>= unsafeLinearWrite marr i
+  unsafeFreeze (getComp arr) marr
 {-# INLINE itraversePrim #-}
 
 --------------------------------------------------------------------------------
@@ -490,7 +561,7 @@
 --
 -- @since 0.1.0
 mapM_ :: (Source r a, Index ix, Monad m) => (a -> m b) -> Array r ix a -> m ()
-mapM_ f !arr = iterM_ zeroIndex (unSz (size arr)) (pureIndex 1) (<) (f . unsafeIndex arr)
+mapM_ = traverseA_
 {-# INLINE mapM_ #-}
 
 
@@ -509,13 +580,31 @@
 -- 499500
 --
 forM_ :: (Source r a, Index ix, Monad m) => Array r ix a -> (a -> m b) -> m ()
-forM_ = flip mapM_
+forM_ = flip traverseA_
 {-# INLINE forM_ #-}
 
 
+-- | Map a monadic index aware function over an array sequentially, while discarding the result.
+--
+-- ==== __Examples__
+--
+-- >>> import Data.Massiv.Array
+-- >>> imapM_ (curry print) $ range Seq (Ix1 10) 15
+-- (0,10)
+-- (1,11)
+-- (2,12)
+-- (3,13)
+-- (4,14)
+--
+-- @since 0.1.0
+imapM_ :: (Index ix, Source r a, Monad m) => (ix -> a -> m b) -> Array r ix a -> m ()
+imapM_ = itraverseA_
+{-# INLINE imapM_ #-}
+
+
 -- | Just like `imapM_`, except with flipped arguments.
 iforM_ :: (Source r a, Index ix, Monad m) => Array r ix a -> (ix -> a -> m b) -> m ()
-iforM_ = flip imapM_
+iforM_ = flip itraverseA_
 {-# INLINE iforM_ #-}
 
 
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
@@ -4,7 +4,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Slice
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Ops/Sort.hs b/src/Data/Massiv/Array/Ops/Sort.hs
--- a/src/Data/Massiv/Array/Ops/Sort.hs
+++ b/src/Data/Massiv/Array/Ops/Sort.hs
@@ -4,7 +4,7 @@
 {-# LANGUAGE MonoLocalBinds #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Sort
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -15,6 +15,9 @@
   , quicksort
   , quicksortBy
   , quicksortByM
+  , quicksortAs
+  , quicksortAsBy
+  , quicksortAsByM
   , quicksortM_
   , quicksortByM_
   , unsafeUnstablePartitionRegionM
@@ -98,11 +101,45 @@
 {-# INLINE unsafeUnstablePartitionRegionM #-}
 
 
--- | This is an implementation of [Quicksort](https://en.wikipedia.org/wiki/Quicksort), which is an
--- efficient, but unstable sort that uses Median-of-three for pivot choosing, as such it performs
--- very well not only for random values, but also for common edge cases like already sorted,
--- reversed sorted and arrays with many duplicate elements. It will also respect the computation
--- strategy and will result in a nice speed up for systems with multiple CPUs.
+-- | Same as `quicksort` except it accepts any array that is computable.
+--
+-- @since 1.0.2
+quicksortAs ::
+     (Load r Ix1 e, Manifest r' e, Ord e) => r' -> Vector r e -> Vector r' e
+quicksortAs _ arr = unsafePerformIO $ withLoadMArray_ arr quicksortM_
+{-# INLINE quicksortAs #-}
+
+-- | Same as `quicksortBy` except it accepts any array that is computable.
+--
+-- @since 1.0.2
+quicksortAsBy ::
+     (Load r Ix1 e, Manifest r' e) => r' -> (e -> e -> Ordering) -> Vector r e -> Vector r' e
+quicksortAsBy _ f arr =
+  unsafePerformIO $ withLoadMArray_ arr (quicksortByM_ (\x y -> pure $ f x y))
+{-# INLINE quicksortAsBy #-}
+
+
+-- | Same as `quicksortByM` except it accepts any array that is computable.
+--
+-- @since 1.0.2
+quicksortAsByM ::
+     (Load r Ix1 e, Manifest r' e, MonadUnliftIO m)
+  => r'
+  -> (e -> e -> m Ordering)
+  -> Vector r e
+  -> m (Vector r' e)
+quicksortAsByM _ f arr =
+  withRunInIO $ \run -> withLoadMArray_ arr (quicksortByM_ (\x y -> run (f x y)))
+{-# INLINE quicksortAsByM #-}
+
+
+-- | This is an implementation of
+-- [Quicksort](https://en.wikipedia.org/wiki/Quicksort), which is an efficient,
+-- but unstable sort. This implementation uses Median-of-three for pivot
+-- choosing, as such it performs very well not only for random values, but also
+-- for common edge cases like already sorted, reversed sorted and arrays with
+-- many duplicate elements. It will also respect the computation strategy and
+-- will result in a nice speed up for systems with multiple CPUs.
 --
 -- @since 0.3.2
 quicksort ::
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
@@ -6,7 +6,7 @@
 {-# OPTIONS_GHC -fno-warn-redundant-constraints #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Transform
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -542,10 +542,10 @@
   let load :: Loader e
       load scheduler !startAt dlWrite _dlSet = do
         scheduleWork scheduler $
-          iterM_ zeroIndex (unSz sz1) (pureIndex 1) (<) $ \ix ->
+          iterA_ zeroIndex (unSz sz1) (pureIndex 1) (<) $ \ix ->
             dlWrite (startAt + toLinearIndex newSz ix) (unsafeIndex arr1 ix)
         scheduleWork scheduler $
-          iterM_ zeroIndex (unSz sz2) (pureIndex 1) (<) $ \ix ->
+          iterA_ 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)
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -193,7 +193,7 @@
       DArray
         (getComp arr)
         sz
-        (stencilF (borderIndex border arr) (borderIndex border arr) . liftIndex2 (+) offset)
+        (PrefIndex (stencilF (borderIndex border arr) (borderIndex border arr) . liftIndex2 (+) offset))
     -- Size by which the resulting array will shrink (not accounting for padding)
     !shrinkSz = Sz (liftIndex (subtract 1) (unSz sSz))
     !sz = liftSz2 (-) (SafeSz (liftIndex2 (+) po (liftIndex2 (+) pb (unSz (size arr))))) shrinkSz
@@ -247,8 +247,7 @@
   -> Stencil ix e a
 makeStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
   where
-    stencil _ getVal !ix =
-      inline relStencil $ \ !ixD -> getVal (liftIndex2 (+) ix ixD)
+    stencil _ getVal !ix = inline (relStencil (getVal . liftIndex2 (+) ix))
     {-# INLINE stencil #-}
 {-# INLINE makeStencil #-}
 
diff --git a/src/Data/Massiv/Array/Stencil/Convolution.hs b/src/Data/Massiv/Array/Stencil/Convolution.hs
--- a/src/Data/Massiv/Array/Stencil/Convolution.hs
+++ b/src/Data/Massiv/Array/Stencil/Convolution.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Convolution
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -6,7 +6,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE RecordWildCards #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -39,7 +39,7 @@
 makeUnsafeStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
   where
     stencil unsafeGetVal _getVal !ix =
-      inline $ relStencil ix (unsafeGetVal . liftIndex2 (+) ix)
+      inline (relStencil ix (unsafeGetVal . liftIndex2 (+) ix))
     {-# INLINE stencil #-}
 {-# INLINE makeUnsafeStencil #-}
 
diff --git a/src/Data/Massiv/Array/Unsafe.hs b/src/Data/Massiv/Array/Unsafe.hs
--- a/src/Data/Massiv/Array/Unsafe.hs
+++ b/src/Data/Massiv/Array/Unsafe.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 -- |
 -- Module      : Data.Massiv.Array.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -117,8 +117,7 @@
 
 unsafeBackpermute :: (Index ix', Source r' e, Index ix) =>
                      Sz ix -> (ix -> ix') -> Array r' ix' e -> Array D ix e
-unsafeBackpermute !sz ixF !arr =
-  makeArray (getComp arr) sz $ \ !ix -> unsafeIndex arr (ixF ix)
+unsafeBackpermute !sz ixF !arr = makeArray (getComp arr) sz (unsafeIndex arr . ixF)
 {-# INLINE unsafeBackpermute #-}
 
 -- | Same 'Data.Array.transform'', except no bounds checking is performed, thus making it faster,
diff --git a/src/Data/Massiv/Core.hs b/src/Data/Massiv/Core.hs
--- a/src/Data/Massiv/Core.hs
+++ b/src/Data/Massiv/Core.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      : Data.Massiv.Core
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -16,6 +16,7 @@
   , Load(iterArrayLinearST_, iterArrayLinearWithSetST_)
   , Stream(..)
   , Source
+  , PrefIndex(..)
   , Size
   , Shape(..)
   , LengthHint(..)
@@ -36,7 +37,6 @@
   , initWorkerStates
   , scheduleWork
   , scheduleWork_
-  , withMassivScheduler_
   , module Data.Massiv.Core.Index
   -- * Numeric
   , FoldNumeric
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
@@ -10,7 +10,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Core.Common
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -26,6 +26,7 @@
   , Stream(..)
   , Strategy(..)
   , Source(..)
+  , PrefIndex(..)
   , Load(..)
   , StrideLoad(..)
   , Size(..)
@@ -37,7 +38,6 @@
   , numWorkers
   , scheduleWork
   , scheduleWork_
-  , withMassivScheduler_
   , WorkerStates
   , unsafeRead
   , unsafeWrite
@@ -72,7 +72,6 @@
   , inline2
   , module Data.Massiv.Core.Index
   -- * Common Operations
-  , imapM_
   , Semigroup((<>))
   -- * Exceptions
   , MonadThrow(..)
@@ -99,9 +98,7 @@
 import Control.Monad.Primitive
 import Control.Monad.ST
 import Control.Scheduler (Comp(..), Scheduler, WorkerStates, numWorkers,
-                          scheduleWork, scheduleWork_, trivialScheduler_,
-                          withScheduler_)
-import Control.Scheduler.Global
+                          scheduleWork, scheduleWork_, trivialScheduler_)
 import GHC.Exts (IsList)
 import Data.Massiv.Core.Exception
 import Data.Massiv.Core.Index
@@ -111,8 +108,6 @@
 import qualified Data.Vector.Fusion.Stream.Monadic as S (Stream)
 import Data.Vector.Fusion.Util
 
-#include "massiv.h"
-
 -- | The array family. Representations @r@ describe how data is arranged or computed. All
 -- arrays have a common property that each index @ix@ always maps to the same unique
 -- element @e@, even if that element does not yet exist in memory and the array has to be
@@ -182,7 +177,17 @@
   -- @since 0.1.0
   getComp :: Array r ix e -> Comp
 
+  -- | Array representation. Representation is never evaluated in @massiv@,
+  -- therefore default implementation is bottom. However, it is recommended to
+  -- supply a constructor that doesn't result in an error when evaluated.
+  --
+  -- @since 1.0.2
+  repr :: r
+  repr =
+    error $ "Array representation should never be evaluated: " ++
+            show (typeRep (Proxy :: Proxy r))
 
+
 -- | Size hint
 --
 -- @since 1.0.0
@@ -278,7 +283,18 @@
   -- @since 0.1.0
   unsafeResize :: (Index ix, Index ix') => Sz ix' -> Array r ix e -> Array r ix' e
 
+-- | Prefered indexing function.
+data PrefIndex ix e
+  = PrefIndex (ix -> e)
+  | PrefIndexLinear (Int -> e)
 
+instance Functor (PrefIndex ix) where
+  fmap f = \case
+    PrefIndex ig -> PrefIndex (f . ig)
+    PrefIndexLinear ig -> PrefIndexLinear (f . ig)
+  {-# INLINE fmap #-}
+  (<$) e _ = PrefIndexLinear (const e)
+  {-# INLINE (<$) #-}
 
 -- | Arrays that can be used as source to practically any manipulation function.
 class (Strategy r, Size r) => Source r e where
@@ -289,9 +305,7 @@
   --
   -- @since 0.1.0
   unsafeIndex :: Index ix => Array r ix e -> ix -> e
-  unsafeIndex =
-    INDEX_CHECK("(Source r e).unsafeIndex",
-                size, \ !arr -> unsafeLinearIndex arr . toLinearIndex (size arr))
+  unsafeIndex !arr = unsafeLinearIndex arr . toLinearIndex (size arr)
   {-# INLINE unsafeIndex #-}
 
   -- | Lookup element in the array using flat index in a row-major fashion. No
@@ -302,7 +316,15 @@
   unsafeLinearIndex !arr = unsafeIndex arr . fromLinearIndex (size arr)
   {-# INLINE unsafeLinearIndex #-}
 
+  -- | Alternative indexing function that can choose an index that is most
+  -- efficient for underlying representation
+  --
+  -- @since 1.0.2
+  unsafePrefIndex :: Index ix => Array r ix e -> PrefIndex ix e
+  unsafePrefIndex !arr = PrefIndexLinear (unsafeLinearIndex arr)
+  {-# INLINE unsafePrefIndex #-}
 
+
   -- | /O(1)/ - Take a slice out of an array from the outside
   --
   -- @since 0.1.0
@@ -317,10 +339,9 @@
   -- @since 0.5.0
   unsafeLinearSlice :: Index ix => Ix1 -> Sz1 -> Array r ix e -> Array r Ix1 e
 
-
 -- | Any array that can be computed and loaded into memory
 class (Strategy r, Shape r ix) => Load r ix e where
-  {-# MINIMAL (makeArray | makeArrayLinear), (iterArrayLinearST_ | iterArrayLinearWithSetST_)#-}
+  {-# MINIMAL (makeArray | makeArrayLinear), (iterArrayLinearST_ | iterArrayLinearWithSetST_) #-}
 
   -- | Construct an Array. Resulting type either has to be unambiguously inferred or restricted
   -- manually, like in the example below. Use "Data.Massiv.Array.makeArrayR" if you'd like to
@@ -387,14 +408,12 @@
     -> ST s ()
   iterArrayLinearST_ scheduler arr uWrite =
     iterArrayLinearWithSetST_ scheduler arr uWrite $ \offset sz e ->
-      loopM_ offset (< (offset + unSz sz)) (+1) (`uWrite` e)
+      loopA_ offset (< (offset + unSz sz)) (+1) (`uWrite` e)
   {-# INLINE iterArrayLinearST_ #-}
 
   -- | Similar to `iterArrayLinearST_`. Except it also accepts a function that is
   -- potentially optimized for setting many cells in a region to the same
-  -- value. There is no guarantees, but some array representations, might
-  -- utilize this region setting function, in which case for such regions index
-  -- aware action will not be called.
+  -- value.
   --
   -- @since 1.0.0
   iterArrayLinearWithSetST_
@@ -453,6 +472,7 @@
                else unsafeLinearGrow mvec k
 {-# INLINE resizeMVector #-}
 
+
 class Load r ix e => StrideLoad r ix e where
   -- | Load an array into memory with stride. Default implementation requires an instance of
   -- `Source`.
@@ -472,12 +492,12 @@
     -> (Int -> e -> ST s ())
     -> ST s ()
   iterArrayLinearWithStrideST_ scheduler stride resultSize arr =
-    splitLinearlyWith_ scheduler (totalElem resultSize) unsafeLinearWriteWithStride
+    splitLinearlyWith_ scheduler (totalElem resultSize) unsafeLinearIndexWithStride
     where
       !strideIx = unStride stride
-      unsafeLinearWriteWithStride =
+      unsafeLinearIndexWithStride =
         unsafeIndex arr . liftIndex2 (*) strideIx . fromLinearIndex resultSize
-      {-# INLINE unsafeLinearWriteWithStride #-}
+      {-# INLINE unsafeLinearIndexWithStride #-}
   {-# INLINE iterArrayLinearWithStrideST_ #-}
 
 -- class (Load r ix e) => StrideLoad r ix e where
@@ -490,9 +510,10 @@
   --   -> MArray RealWorld r' ix e
   --   -> m (MArray RealWorld r' ix e)
 
--- | Starting with massiv-1.0 `Mutable` and `Manifest` are synonymous. However,
--- this type class synonym will be deprecated in the next major version.
+-- | Starting with massiv-1.0 `Mutable` and `Manifest` are synonymous. Since massiv-1.1
+-- it is deprecated and will be removed in massiv-1.2
 type Mutable r e = Manifest r e
+{-# DEPRECATED Mutable "In favor of `Manifest`" #-}
 
 -- | Manifest arrays are backed by actual memory and values are looked up versus
 -- computed as it is with delayed arrays. Because manifest arrays are located in
@@ -573,7 +594,7 @@
   unsafeLinearSet :: (Index ix, PrimMonad m) =>
                      MArray (PrimState m) r ix e -> Ix1 -> Sz1 -> e -> m ()
   unsafeLinearSet marr offset len e =
-    loopM_ offset (< (offset + unSz len)) (+1) (\i -> unsafeLinearWrite marr i e)
+    loopA_ offset (< (offset + unSz len)) (+1) (\i -> unsafeLinearWrite marr i e)
   {-# INLINE unsafeLinearSet #-}
 
   -- | Copy part of one mutable array into another
@@ -588,7 +609,7 @@
                    -> m ()
   unsafeLinearCopy marrFrom iFrom marrTo iTo (SafeSz k) = do
     let delta = iTo - iFrom
-    loopM_ iFrom (< k + iFrom) (+1) $ \i ->
+    loopA_ iFrom (< k + iFrom) (+1) $ \i ->
       unsafeLinearRead marrFrom i >>= unsafeLinearWrite marrTo (i + delta)
   {-# INLINE unsafeLinearCopy #-}
 
@@ -604,7 +625,7 @@
                         -> m ()
   unsafeArrayLinearCopy arrFrom iFrom marrTo iTo (SafeSz k) = do
     let delta = iTo - iFrom
-    loopM_ iFrom (< k + iFrom) (+1) $ \i ->
+    loopA_ iFrom (< k + iFrom) (+1) $ \i ->
       unsafeLinearWrite marrTo (i + delta) (unsafeLinearIndex arrFrom i)
   {-# INLINE unsafeArrayLinearCopy #-}
 
@@ -643,19 +664,6 @@
   pure marr'
 {-# INLINE unsafeDefaultLinearShrink #-}
 
-
--- | Selects an optimal scheduler for the supplied strategy, but it works only in `IO`
---
--- @since 1.0.0
-withMassivScheduler_ :: Comp -> (Scheduler RealWorld () -> IO ()) -> IO ()
-withMassivScheduler_ comp f =
-  case comp of
-    Par -> withGlobalScheduler_ globalScheduler f
-    Seq -> f trivialScheduler_
-    _   -> withScheduler_ comp f
-{-# INLINE withMassivScheduler_ #-}
-
-
 -- | Read an array element
 --
 -- @since 0.1.0
@@ -975,25 +983,6 @@
 evaluate' :: (HasCallStack, Index ix, Source r e) => Array r ix e -> ix -> e
 evaluate' arr ix = throwEither (evaluateM arr ix)
 {-# INLINE evaluate' #-}
-
-
--- | Map a monadic index aware function over an array sequentially, while discarding the result.
---
--- ==== __Examples__
---
--- >>> import Data.Massiv.Array
--- >>> imapM_ (curry print) $ range Seq (Ix1 10) 15
--- (0,10)
--- (1,11)
--- (2,12)
--- (3,13)
--- (4,14)
---
--- @since 0.1.0
-imapM_ :: (Index ix, Source r a, Monad m) => (ix -> a -> m b) -> Array r ix a -> m ()
-imapM_ f !arr =
-  iterM_ zeroIndex (unSz (size arr)) (pureIndex 1) (<) $ \ !ix -> f ix (unsafeIndex arr ix)
-{-# INLINE imapM_ #-}
 
 
 
diff --git a/src/Data/Massiv/Core/Exception.hs b/src/Data/Massiv/Core/Exception.hs
--- a/src/Data/Massiv/Core/Exception.hs
+++ b/src/Data/Massiv/Core/Exception.hs
@@ -5,7 +5,7 @@
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Massiv.Core.Exception
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Core/Index.hs b/src/Data/Massiv/Core/Index.hs
--- a/src/Data/Massiv/Core/Index.hs
+++ b/src/Data/Massiv/Core/Index.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE ExplicitNamespaces #-}
 -- |
 -- Module      : Data.Massiv.Core.Index
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
@@ -86,9 +86,12 @@
   , insertDimension
   -- * Iterators
   , iter
+  , iterA_
+  , iterM_
   , iterLinearM
   , iterLinearM_
-  , module Data.Massiv.Core.Iterator
+  , module Data.Massiv.Core.Loop
+  , module Data.Massiv.Core.Index.Iterator
   , module Data.Massiv.Core.Index.Tuple
   -- * Exceptions
   , IndexException(..)
@@ -107,7 +110,8 @@
 import Data.Massiv.Core.Index.Ix
 import Data.Massiv.Core.Index.Stride
 import Data.Massiv.Core.Index.Tuple
-import Data.Massiv.Core.Iterator
+import Data.Massiv.Core.Index.Iterator
+import Data.Massiv.Core.Loop
 import GHC.TypeLits
 
 
@@ -542,7 +546,7 @@
 -- ==== __Examples__
 --
 -- >>> sz = Sz2 3 4
--- >>> iterLinearM sz 0 3 1 (<) 100 $ \ k ix acc -> print (fromLinearIndex sz k == ix) >> pure (acc + k)
+-- >>> iterLinearM sz 0 3 1 (<) 100 $ \ k ix acc -> (acc + k) <$ print (fromLinearIndex sz k == ix)
 -- True
 -- True
 -- True
@@ -582,9 +586,8 @@
              -> (Int -> ix -> m ()) -- ^ Monadic action that takes index in both forms
              -> m ()
 iterLinearM_ sz !k0 !k1 !inc cond f =
-  loopM_ k0 (`cond` k1) (+ inc) $ \ !i -> f i (fromLinearIndex sz i)
+  loopA_ k0 (`cond` k1) (+ inc) $ \ !i -> f i (fromLinearIndex sz i)
 {-# INLINE iterLinearM_ #-}
-
 
 -- | This is used by @INDEX_CHECK@ macro and thus used whenever the @unsafe-checks@ cabal
 -- flag is on.
diff --git a/src/Data/Massiv/Core/Index/Internal.hs b/src/Data/Massiv/Core/Index/Internal.hs
--- a/src/Data/Massiv/Core/Index/Internal.hs
+++ b/src/Data/Massiv/Core/Index/Internal.hs
@@ -10,13 +10,14 @@
 {-# LANGUAGE LambdaCase #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
@@ -51,6 +52,8 @@
   , ReportInvalidDim
   , Lower
   , Index(..)
+  , iterA_
+  , iterM_
   , Ix0(..)
   , type Ix1
   , pattern Ix1
@@ -60,13 +63,15 @@
   , showsPrecWrapped
   ) where
 
+import Control.Monad.ST
+import Control.Scheduler
 import Control.DeepSeq
 import Control.Exception (Exception(..), throw)
-import Control.Monad (when)
+import Control.Monad (when, void)
 import Control.Monad.Catch (MonadThrow(..))
 import Data.Coerce
 import Data.Kind
-import Data.Massiv.Core.Iterator
+import Data.Massiv.Core.Loop
 import Data.Typeable
 import GHC.TypeLits
 import System.Random.Stateful
@@ -123,7 +128,7 @@
 --
 -- @since 0.3.0
 pattern Sz1 :: Ix1 -> Sz Ix1
-pattern Sz1 ix  <- SafeSz ix where
+pattern Sz1 ix <- SafeSz ix where
         Sz1 ix = SafeSz (max 0 ix)
 {-# COMPLETE Sz1 #-}
 
@@ -160,7 +165,9 @@
   negate x
     | x == zeroSz = x
     | otherwise =
-      error $ "Attempted to negate: " ++ show x ++ ", this can lead to unexpected behavior. See https://github.com/lehins/massiv/issues/114"
+      error $
+      "Attempted to negate: " ++ show x ++
+      ", this can lead to unexpected behavior. See https://github.com/lehins/massiv/issues/114"
   {-# INLINE negate #-}
   signum x = SafeSz (signum (coerce x))
   {-# INLINE signum #-}
@@ -179,6 +186,7 @@
         when (acc' /= 0 && acc' < acc) $ throwM $ SizeOverflowException (SafeSz ix)
         pure acc'
   Sz ix <$ foldlIndex (\acc i -> acc >>= guardNegativeOverflow i) (pure 1) ix
+{-# INLINE mkSzM #-}
 
 
 
@@ -565,7 +573,7 @@
   {-# INLINE [1] toLinearIndex #-}
 
   -- | Convert linear index from size and index with an accumulator. Currently is useless and will
-  -- likley be removed in future versions.
+  -- likely be removed in future versions.
   --
   -- @since 0.1.0
   toLinearIndexAcc :: Ix1 -> ix -> ix -> Ix1
@@ -581,9 +589,9 @@
   -- @since 0.1.0
   fromLinearIndex :: Sz ix -> Ix1 -> ix
   default fromLinearIndex :: Index (Lower ix) => Sz ix -> Ix1 -> ix
-  fromLinearIndex (SafeSz sz) k = consDim q ixL
+  fromLinearIndex (SafeSz sz) !k = consDim q ixL
     where
-      !(q, ixL) = fromLinearIndexAcc (snd (unconsDim sz)) k
+      !(!q, !ixL) = fromLinearIndexAcc (snd (unconsDim sz)) k
   {-# INLINE [1] fromLinearIndex #-}
 
   -- | Compute an index from size and linear index using an accumulator, thus trying to optimize for
@@ -592,11 +600,11 @@
   -- @since 0.1.0
   fromLinearIndexAcc :: ix -> Ix1 -> (Int, ix)
   default fromLinearIndexAcc :: Index (Lower ix) => ix -> Ix1 -> (Ix1, ix)
-  fromLinearIndexAcc ix' !k = (q, consDim r ixL)
+  fromLinearIndexAcc !ix' !k = (q, consDim r ixL)
     where
-      !(m, ix) = unconsDim ix'
-      !(kL, ixL) = fromLinearIndexAcc ix k
-      !(q, r) = quotRem kL m
+      !(!m, !ix) = unconsDim ix'
+      !(!kL, !ixL) = fromLinearIndexAcc ix k
+      !(!q, !r) = quotRem kL m
   {-# INLINE [1] fromLinearIndexAcc #-}
 
   -- | A way to make sure index is withing the bounds for the supplied size. Takes two functions
@@ -641,21 +649,229 @@
       !(inc, incIxL) = unconsDim incIx
   {-# INLINE iterM #-}
 
-  -- TODO: Implement in terms of iterM, benchmark it and remove from `Index`
-  -- | Same as `iterM`, but don't bother with accumulator and return value.
+  iterRowMajorST :: Int -- ^ Scheduler multiplying factor. Must be positive
+                 -> Scheduler s a -- ^ Scheduler to use
+                 -> ix -- ^ Start index
+                 -> ix -- ^ Stride
+                 -> Sz ix -- ^ Size
+                 -> a -- ^ Initial accumulator
+                 -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+                 -> (ix -> a -> ST s a) -- ^ Action
+                 -> ST s a
+  default iterRowMajorST :: Index (Lower ix)
+                         => Int
+                         -> Scheduler s a
+                         -> ix
+                         -> ix
+                         -> Sz ix
+                         -> a
+                         -> (a -> ST s (a, a))
+                         -> (ix -> a -> ST s a)
+                         -> ST s a
+  iterRowMajorST !fact scheduler ixStart ixStride sz initAcc splitAcc f = do
+    let !(SafeSz n, szL@(SafeSz nL)) = unconsSz sz
+    if n > 0
+      then do
+        let !(!start, !ixL) = unconsDim ixStart
+            !(!stride, !sL) = unconsDim ixStride
+        if numWorkers scheduler > 1 && fact > 1 && n < numWorkers scheduler * fact
+          then do
+            let !newFact = 1 + (fact `quot` n)
+            loopM start (< start + n * stride) (+ stride) initAcc $ \j acc ->
+              iterRowMajorST newFact scheduler ixL sL szL acc splitAcc (f . consDim j)
+          else
+            splitWorkWithFactorST fact scheduler start stride n initAcc splitAcc $
+              \ _ _ chunkStartAdj chunkStopAdj acc ->
+                loopM chunkStartAdj (< chunkStopAdj) (+ stride) acc $ \j a ->
+                  iterM ixL nL sL (<) a (f . consDim j)
+      else pure initAcc
+  {-# INLINE iterRowMajorST #-}
+
+  -- | Similar to `iterM`, but no restriction on a Monad.
   --
+  -- @since 1.0.2
+  iterF :: ix -> ix -> ix -> (Int -> Int -> Bool) -> f a -> (ix -> f a -> f a) -> f a
+  default iterF :: (Index (Lower ix)) =>
+    ix -> ix -> ix -> (Int -> Int -> Bool) -> f a -> (ix -> f a -> f a) -> f a
+  iterF !sIx !eIx !incIx cond initAct f =
+    loopF s (`cond` e) (+ inc) initAct $ \ !i g ->
+      iterF sIxL eIxL incIxL cond g (\ !ix -> f (consDim i ix))
+    where
+      !(s, sIxL) = unconsDim sIx
+      !(e, eIxL) = unconsDim eIx
+      !(inc, incIxL) = unconsDim incIx
+  {-# INLINE iterF #-}
+
+  -- | A single step in iteration
+  --
   -- @since 0.1.0
-  iterM_ :: Monad m => ix -> ix -> ix -> (Int -> Int -> Bool) -> (ix -> m a) -> m ()
-  default iterM_ :: (Index (Lower ix), Monad m) =>
-    ix -> ix -> ix -> (Int -> Int -> Bool) -> (ix -> m a) -> m ()
-  iterM_ !sIx eIx !incIx cond f =
-    loopM_ s (`cond` e) (+ inc) $ \ !i -> iterM_ sIxL eIxL incIxL cond $ \ !ix -> f (consDim i ix)
+  stepNextMF :: ix -> ix -> ix -> (Int -> Int -> Bool) -> (Maybe ix -> f a) -> f a
+  default stepNextMF :: (Index (Lower ix)) =>
+    ix -> ix -> ix -> (Int -> Int -> Bool) -> (Maybe ix -> f a) -> f a
+  stepNextMF !sIx !eIx !incIx cond f =
+    nextMaybeF s (`cond` e) (+ inc) $ \ !mni ->
+      stepNextMF sIxL eIxL incIxL cond $ \ mIxN ->
+        f $!
+          case mIxN of
+            Just ixN -> Just $! consDim s ixN
+            Nothing ->
+              case mni of
+                Just ni -> Just $! consDim ni (pureIndex 0)
+                Nothing -> Nothing
     where
       !(s, sIxL) = unconsDim sIx
       !(e, eIxL) = unconsDim eIx
       !(inc, incIxL) = unconsDim incIx
-  {-# INLINE iterM_ #-}
+  {-# INLINE stepNextMF #-}
 
+  iterTargetRowMajorA_ :: Applicative f
+                       => Int -- ^ Target linear index accumulator
+                       -> Int -- ^ Target linear index start
+                       -> Sz ix -- ^ Target size
+                       -> ix -- ^ Source start index
+                       -> ix -- ^ Source stride
+                       -> (Ix1 -> ix -> f a)
+                       -- ^ Action that accepts a linear index of the target,
+                       -- multi-dimensional index of the source and accumulator
+                       -> f ()
+  default iterTargetRowMajorA_ :: (Applicative f, Index (Lower ix))
+                        => Int
+                        -> Int
+                        -> Sz ix
+                        -> ix
+                        -> ix
+                        -> (Ix1 -> ix -> f a)
+                        -> f ()
+  iterTargetRowMajorA_ !iAcc !iStart szRes ixStart ixStride f = do
+    let !(SafeSz nRes, !szL) = unconsSz szRes
+        !(!start, !ixL) = unconsDim ixStart
+        !(!stride, !sL) = unconsDim ixStride
+    iloopA_ (iAcc * nRes) start (< start + nRes * stride) (+ stride) $ \k j ->
+      iterTargetRowMajorA_ k iStart szL ixL sL $ \i jl -> f i (consDim j jl)
+  {-# INLINE iterTargetRowMajorA_ #-}
+
+  iterTargetRowMajorAccM :: Monad m =>
+                            Int -- ^ Target linear index accumulator
+                         -> Int -- ^ Target linear index start
+                         -> Sz ix -- ^ Target size
+                         -> ix -- ^ Source start index
+                         -> ix -- ^ Source stride
+                         -> a -- ^ Accumulator
+                         -> (Ix1 -> ix -> a -> m a)
+                         -- ^ Action that accepts a linear index of the target,
+                         -- multi-dimensional index of the source and accumulator
+                         -> m a
+  default iterTargetRowMajorAccM :: (Monad m, Index (Lower ix))
+                        => Int
+                        -> Int
+                        -> Sz ix
+                        -> ix
+                        -> ix
+                        -> a
+                        -> (Ix1 -> ix -> a -> m a)
+                        -> m a
+  iterTargetRowMajorAccM !iAcc !iStart szRes ixStart ixStride initAcc f = do
+    let !(SafeSz nRes, !szL) = unconsSz szRes
+        !(!start, !ixL) = unconsDim ixStart
+        !(!stride, !sL) = unconsDim ixStride
+    iloopM (iAcc * nRes) start (< start + nRes * stride) (+ stride) initAcc $ \k j acc ->
+      iterTargetRowMajorAccM k iStart szL ixL sL acc $ \i jl -> f i (consDim j jl)
+  {-# INLINE iterTargetRowMajorAccM #-}
+
+  iterTargetRowMajorAccST ::
+       Int -- ^ Linear index accumulator
+    -> Int -- ^ Scheduler multiplying factor. Must be positive
+    -> Scheduler s a -- ^ Scheduler to use
+    -> Int -- ^ Target linear index start
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Source start index
+    -> ix -- ^ Source stride
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s a
+  default iterTargetRowMajorAccST :: Index (Lower ix)
+                                  => Int
+                                  -> Int
+                                  -> Scheduler s a
+                                  -> Int
+                                  -> Sz ix
+                                  -> ix
+                                  -> ix
+                                  -> a
+                                  -> (a -> ST s (a, a))
+                                  -> (Ix1 -> ix -> a -> ST s a)
+                                  -> ST s a
+  iterTargetRowMajorAccST !iAcc !fact scheduler iStart sz ixStart ixStride initAcc splitAcc f = do
+    let !(SafeSz n, nL) = unconsSz sz
+    if n > 0
+      then do
+        let !(!start, !ixL) = unconsDim ixStart
+            !(!stride, !sL) = unconsDim ixStride
+            !iAccL = iAcc * n
+        if numWorkers scheduler > 1 && fact > 1 && n < numWorkers scheduler * fact
+          then do
+            let newFact = 1 + (fact `quot` n)
+            iloopM iAccL start (< start + n * stride) (+ stride) initAcc $ \k j acc -> do
+              iterTargetRowMajorAccST k newFact scheduler iStart nL ixL sL acc splitAcc $ \i ->
+                f i . consDim j
+          else
+            splitWorkWithFactorST fact scheduler start stride n initAcc splitAcc $
+              \ chunkStart _ chunkStartAdj chunkStopAdj acc ->
+                iloopM (iAccL + chunkStart) chunkStartAdj (< chunkStopAdj) (+ stride) acc $ \k j a ->
+                  iterTargetRowMajorAccM k iStart nL ixL sL a $ \i -> f i . consDim j
+      else pure initAcc
+  {-# INLINE iterTargetRowMajorAccST #-}
+
+
+  iterTargetRowMajorAccST_
+    :: Int -- ^ Index accumulator
+    -> Int -- ^ Scheduler multiplying factor. Must be positive
+    -> Scheduler s () -- ^ Scheduler to use
+    -> Int -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Source start index
+    -> ix -- ^ Source stride
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s ()
+  default iterTargetRowMajorAccST_
+    :: Index (Lower ix)
+    => Int
+    -> Int
+    -> Scheduler s ()
+    -> Int
+    -> Sz ix
+    -> ix
+    -> ix
+    -> a
+    -> (a -> ST s (a, a))
+    -> (Ix1 -> ix -> a -> ST s a)
+    -> ST s ()
+  iterTargetRowMajorAccST_ !iAcc !fact scheduler iStart sz ixStart ixStride initAcc splitAcc f = do
+    let !(SafeSz n, szL) = unconsSz sz
+    when (n > 0) $ do
+      let !(!start, !ixL) = unconsDim ixStart
+          !(!stride, !sL) = unconsDim ixStride
+          !iAccL = iAcc * n
+      if numWorkers scheduler > 1 && fact > 1 && n < numWorkers scheduler * fact
+        then do
+          let !newFact = 1 + (fact `quot` n)
+          void $ iloopM iAccL start (< n * stride) (+ stride) initAcc $ \k j acc -> do
+            (accCur, accNext) <- splitAcc acc
+            scheduleWork_ scheduler $
+              iterTargetRowMajorAccST_ k newFact scheduler iStart szL ixL sL accCur splitAcc $ \i ->
+                f i . consDim j
+            pure accNext
+         else
+           void $ splitWorkWithFactorST fact scheduler start stride n initAcc splitAcc $
+             \ chunkStart _ chunkStartAdj chunkStopAdj acc ->
+               void $
+               iloopM (iAccL + chunkStart) chunkStartAdj (< chunkStopAdj) (+ stride) acc $ \k j a ->
+                 iterTargetRowMajorAccM k iStart szL ixL sL a $ \i -> f i . consDim j
+  {-# INLINE iterTargetRowMajorAccST_ #-}
+
 -- | Zero-dimension, i.e. a scalar. Can't really be used directly as there is no instance of
 -- `Index` for it, and is included for completeness.
 data Ix0 = Ix0 deriving (Eq, Ord, Show)
@@ -697,7 +913,7 @@
   {-# INLINE [1] isSafeIndex #-}
   toLinearIndex _ = id
   {-# INLINE [1] toLinearIndex #-}
-  toLinearIndexAcc !acc m i  = acc * m + i
+  toLinearIndexAcc !acc m i = acc * m + i
   {-# INLINE [1] toLinearIndexAcc #-}
   fromLinearIndex _ = id
   {-# INLINE [1] fromLinearIndex #-}
@@ -720,8 +936,8 @@
   getDimM ix 1 = pure ix
   getDimM ix d = throwM $ IndexDimensionException ix d
   {-# INLINE [1] getDimM #-}
-  setDimM _  1 ix = pure ix
-  setDimM ix d _  = throwM $ IndexDimensionException ix d
+  setDimM _ 1 ix = pure ix
+  setDimM ix d _ = throwM $ IndexDimensionException ix d
   {-# INLINE [1] setDimM #-}
   modifyDimM ix 1 f = pure (ix, f ix)
   modifyDimM ix d _ = throwM $ IndexDimensionException ix d
@@ -730,7 +946,7 @@
   pullOutDimM ix d = throwM $ IndexDimensionException ix d
   {-# INLINE [1] pullOutDimM #-}
   insertDimM Ix0 1 i = pure i
-  insertDimM ix  d _ = throwM $ IndexDimensionException ix d
+  insertDimM ix d _ = throwM $ IndexDimensionException ix d
   {-# INLINE [1] insertDimM #-}
   pureIndex i = i
   {-# INLINE [1] pureIndex #-}
@@ -740,10 +956,66 @@
   {-# INLINE [1] liftIndex2 #-}
   foldlIndex f = f
   {-# INLINE [1] foldlIndex #-}
-  iterM k0 k1 inc cond = loopM k0 (`cond` k1) (+inc)
+  iterM k0 k1 inc cond = loopM k0 (`cond` k1) (+ inc)
   {-# INLINE iterM #-}
-  iterM_ k0 k1 inc cond = loopM_ k0 (`cond` k1) (+inc)
-  {-# INLINE iterM_ #-}
+  iterF k0 k1 inc cond = loopF k0 (`cond` k1) (+ inc)
+  {-# INLINE iterF #-}
+  stepNextMF k0 k1 inc cond = nextMaybeF k0 (`cond` k1) (+ inc)
+  {-# INLINE stepNextMF #-}
+
+  iterRowMajorST fact scheduler start step n =
+    iterLinearAccST fact scheduler start step (unSz n)
+  {-# INLINE iterRowMajorST #-}
+
+  iterTargetRowMajorA_ iAcc iStart (SafeSz nRes) start stride =
+    iloopA_ (iAcc * nRes + iStart) start (< start + nRes * stride) (+ stride)
+  {-# INLINE iterTargetRowMajorA_ #-}
+
+  iterTargetRowMajorAccM iAcc iStart (SafeSz nRes) start stride =
+    iloopM (iAcc * nRes + iStart) start (< start + nRes * stride) (+ stride)
+  {-# INLINE iterTargetRowMajorAccM #-}
+
+  iterTargetRowMajorAccST iAcc fact scheduler iStart sz start stride initAcc splitAcc action = do
+    let !n = unSz sz
+        !iAccL = iStart + iAcc * n
+    splitWorkWithFactorST fact scheduler start stride n initAcc splitAcc $
+      \ chunkStart _ chunkStartAdj chunkStopAdj acc ->
+        iloopM (iAccL + chunkStart) chunkStartAdj (< chunkStopAdj) (+ stride) acc action
+  {-# INLINE iterTargetRowMajorAccST #-}
+
+  iterTargetRowMajorAccST_ iAcc fact scheduler iStart sz start stride initAcc splitAcc action = do
+    let !n = unSz sz
+        !iAccL = iStart + iAcc * n
+    void $ splitWorkWithFactorST fact scheduler start stride n initAcc splitAcc $
+      \ chunkStart _ chunkStartAdj chunkStopAdj acc ->
+        void $ iloopM (iAccL + chunkStart) chunkStartAdj (< chunkStopAdj) (+ stride) acc action
+  {-# INLINE iterTargetRowMajorAccST_ #-}
+
+
+-- | Same as `iterM`, but don't bother with accumulator and return value.
+--
+-- @since 0.1.0
+iterM_ :: (Index ix, Monad m) => ix -> ix -> ix -> (Int -> Int -> Bool) -> (ix -> m a) -> m ()
+iterM_ sIx eIx incIx cond f = iterM sIx eIx incIx cond () $ \ !ix !a -> f ix >> pure a
+{-# INLINE iterM_ #-}
+{-# DEPRECATED iterM_ "In favor of more lax `iterA_`" #-}
+
+-- | Same as `iterM`, Iterate over a region with specific step, but using
+-- `Applicative` instead of a `Monad` and don't bother with accumulator or return value.
+--
+-- @since 1.0.2
+iterA_ ::
+     forall ix f a. (Index ix, Applicative f)
+  => ix -- ^ Starting index
+  -> ix -- ^ Ending index (not included)
+  -> ix -- ^ Stepping index
+  -> (Int -> Int -> Bool) -- ^ Continuation function. Loop will stop on `False`
+  -> (ix -> f a) -- ^ Action applied to an index. Result is ignored.
+  -> f ()
+iterA_ sIx eIx incIx cond f =
+  iterF sIx eIx incIx cond (pure ()) $ \ix go -> f ix *> go
+{-# INLINE iterA_ #-}
+
 
 
 -- | Exceptions that get thrown when there is a problem with an index, size or dimension.
diff --git a/src/Data/Massiv/Core/Index/Iterator.hs b/src/Data/Massiv/Core/Index/Iterator.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Massiv/Core/Index/Iterator.hs
@@ -0,0 +1,454 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ExplicitForAll #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MonoLocalBinds #-}
+{-# LANGUAGE PatternSynonyms #-}
+-- |
+-- Module      : Data.Massiv.Core.Index.Iterator
+-- Copyright   : (c) Alexey Kuleshevich 2021-2022
+-- License     : BSD3
+-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
+-- Stability   : experimental
+-- Portability : non-portable
+--
+module Data.Massiv.Core.Index.Iterator
+  ( Iterator(..)
+  -- * Extra iterator functions
+  , iterTargetAccST
+  , iterTargetAccST_
+  , iterTargetFullWithStrideAccST
+  , iterTargetFullWithStrideAccST_
+  , iterTargetST_
+  , iterTargetFullWithStrideST_
+  -- * Iterator implementations
+  , RowMajor(RowMajor)
+  , defRowMajor
+  , RowMajorLinear(RowMajorLinear)
+  , defRowMajorLinear
+  , RowMajorUnbalanced(RowMajorUnbalanced)
+  , defRowMajorUnbalanced
+  ) where
+
+import Control.Monad
+import Control.Monad.ST
+import Control.Scheduler
+import Data.Massiv.Core.Index.Internal
+import Data.Massiv.Core.Index.Stride
+import Data.Massiv.Core.Loop
+
+
+class Iterator it where
+  {-# MINIMAL (iterTargetM, iterTargetA_, iterTargetWithStrideAccST, iterTargetWithStrideAccST_) #-}
+
+
+  -- | Iterate over a target region using linear index with access to the source
+  -- index, which adjusted according to the stride. Use `iterTargetM` if you
+  -- need an accumulator.
+  --
+  -- @since 1.0.2
+  iterTargetA_ ::
+       (Index ix, Applicative f)
+    => it
+    -> Int -- ^ Target linear index start
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Source start index
+    -> Stride ix -- ^ Source stride
+    -> (Ix1 -> ix -> f a)
+    -- ^ Action that accepts a linear index of the target and multi-dimensional
+    -- index of the source.
+    -> f ()
+
+  -- | Iterate over a target region using linear index with access to the source
+  -- index, which adjusted according to the stride.
+  --
+  -- @since 1.0.2
+  iterTargetM ::
+       (Index ix, Monad m)
+    => it
+    -> Ix1 -- ^ Target linear index start
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Source start index
+    -> Stride ix -- ^ Source stride
+    -> a -- ^ Accumulator
+    -> (Ix1 -> ix -> a -> m a)
+    -- ^ Action that accepts a linear index of the target,
+    -- multi-dimensional index of the source and accumulator
+    -> m a
+
+  iterTargetWithStrideAccST ::
+       Index ix
+    => it
+    -> Scheduler s a -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Source start index
+    -> Stride ix -- ^ Source stride
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a))
+    -- ^ Splitting action that produces new accumulators for separate worker threads.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s a
+
+  iterTargetWithStrideAccST_ ::
+       Index ix
+    => it
+    -> Scheduler s () -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Start
+    -> Stride ix -- ^ Stride
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a))
+    -- ^ Splitting action that produces new accumulators for separate worker threads.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s ()
+
+  -- | Iterate over a region with a monadic action and accumulator.
+  --
+  -- @since 1.0.2
+  iterFullM ::
+       (Index ix, Monad m)
+    => it
+    -> ix -- ^ Source start index
+    -> Sz ix -- ^ Source size
+    -> a -- ^ Accumulator
+    -> (ix -> a -> m a)
+    -- ^ Action that accepts a linear index of the target,
+    -- multi-dimensional index of the source and accumulator
+    -> m a
+  iterFullM it start sz acc f =
+    iterTargetM it 0 sz start oneStride acc (const f)
+  {-# INLINE iterFullM #-}
+
+  -- | Iterate over a region with an applicative action ignoring the result.
+  --
+  -- @since 1.0.2
+  iterFullA_ ::
+       (Index ix, Applicative f)
+    => it
+    -> ix -- ^ Source start index
+    -> Sz ix -- ^ Source size
+    -> (ix -> f a)
+    -- ^ Action that accepts a linear index of the target,
+    -- multi-dimensional index of the source and accumulator
+    -> f ()
+  iterFullA_ it start sz f =
+    iterTargetA_ it 0 sz start oneStride (const f)
+  {-# INLINE iterFullA_ #-}
+
+  -- | Iterate over a region in a ST monad with access to `Scheduler`.
+  iterFullAccST ::
+       Index ix
+    => it -- ^ Scheduler multiplying factor. Must be positive
+    -> Scheduler s a -- ^ Scheduler to use
+    -> ix -- ^ Start index
+    -> Sz ix -- ^ Size
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+    -> (ix -> a -> ST s a) -- ^ Action
+    -> ST s a
+  iterFullAccST it scheduler start sz acc splitAcc f =
+    iterTargetAccST it scheduler 0 sz start acc splitAcc (const f)
+  {-# INLINE iterFullAccST #-}
+
+  iterTargetFullAccST ::
+       Index ix
+    => it
+    -> Scheduler s a -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s a
+  iterTargetFullAccST it scheduler iStart sz =
+    iterTargetFullWithStrideAccST it scheduler iStart sz oneStride
+  {-# INLINE iterTargetFullAccST #-}
+
+  iterTargetFullAccST_ ::
+       Index ix
+    => it
+    -> Scheduler s () -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> a -- ^ Initial accumulator
+    -> (a -> ST s (a, a)) -- ^ Function that splits accumulator for each scheduled job.
+    -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+    -> ST s ()
+  iterTargetFullAccST_ it scheduler iStart sz =
+    iterTargetFullWithStrideAccST_ it scheduler iStart sz oneStride
+  {-# INLINE iterTargetFullAccST_ #-}
+
+  iterTargetFullST_ ::
+       Index ix
+    => it
+    -> Scheduler s () -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> (Ix1 -> ix -> ST s ()) -- ^ Action
+    -> ST s ()
+  iterTargetFullST_ it scheduler iStart sz =
+    iterTargetST_ it scheduler iStart sz (pureIndex 0)
+  {-# INLINE iterTargetFullST_ #-}
+
+  -- NOTE: this function does not have to be part of the class, but for some
+  -- reason it creates a severe regression when moved outside.
+  -- | Iterate over a target array with a stride without an accumulator
+  iterTargetWithStrideST_ ::
+       Index ix
+    => it
+    -> Scheduler s () -- ^ Scheduler to use
+    -> Ix1 -- ^ Target linear start index
+    -> Sz ix -- ^ Target size
+    -> ix -- ^ Start
+    -> Stride ix -- ^ Stride
+    -> (Ix1 -> ix -> ST s a) -- ^ Action
+    -> ST s ()
+  iterTargetWithStrideST_ it scheduler i sz ix stride action =
+    iterTargetWithStrideAccST_ it scheduler i sz ix stride () noSplit $ \j jx _ ->
+      void $ action j jx
+  {-# INLINE iterTargetWithStrideST_ #-}
+
+
+-- | Default iterator that parallelizes work in linear chunks. Supplied factor
+-- will be used to schedule that many jobs per capability.
+--
+-- @since 1.0.2
+newtype RowMajor = RowMajorInternal Int
+
+-- | Default row major iterator with multiplying factor set to @8@.
+defRowMajor :: RowMajor
+defRowMajor = RowMajorInternal 8
+
+pattern RowMajor :: Int
+                 -- ^ Multiplier that will be used to scale number of jobs.
+                 -> RowMajor
+pattern RowMajor f <- RowMajorInternal f
+  where RowMajor = RowMajorInternal . max 1
+{-# COMPLETE RowMajor #-}
+
+instance Iterator RowMajor where
+  iterFullM _ start (Sz sz) = iterM start sz (pureIndex 1) (<)
+  {-# INLINE iterFullM #-}
+  iterFullA_ _ start (Sz sz) = iterA_ start sz (pureIndex 1) (<)
+  {-# INLINE iterFullA_ #-}
+  iterFullAccST (RowMajorInternal fact) scheduler startIx =
+    iterRowMajorST fact scheduler startIx (pureIndex 1)
+  {-# INLINE iterFullAccST #-}
+  iterTargetA_ _ i sz start (Stride stride) =
+    iterTargetRowMajorA_ 0 i sz start stride
+  {-# INLINE iterTargetA_ #-}
+  iterTargetM _ i sz start (Stride stride) =
+    iterTargetRowMajorAccM 0 i sz start stride
+  {-# INLINE iterTargetM #-}
+  iterTargetWithStrideAccST (RowMajor fact) scheduler i sz ix (Stride stride) =
+    iterTargetRowMajorAccST 0 fact scheduler i sz ix stride
+  {-# INLINE iterTargetWithStrideAccST #-}
+  iterTargetWithStrideAccST_ (RowMajor fact) scheduler i sz ix (Stride stride) =
+    iterTargetRowMajorAccST_ 0 fact scheduler i sz ix stride
+  {-# INLINE iterTargetWithStrideAccST_ #-}
+
+
+newtype RowMajorLinear = RowMajorLinear Int
+
+defRowMajorLinear :: RowMajorLinear
+defRowMajorLinear = RowMajorLinear 8
+
+instance Iterator RowMajorLinear where
+  iterTargetM _ iStart sz start (Stride stride) acc action =
+    loopM 0 (< totalElem sz) (+ 1) acc $ \i ->
+      action (iStart + i) (liftIndex2 (+) start (liftIndex2 (*) stride (fromLinearIndex sz i)))
+  {-# INLINE iterTargetM #-}
+  iterTargetA_ _ iStart sz start (Stride stride) action =
+    loopA_ 0 (< totalElem sz) (+ 1) $ \i ->
+      action (iStart + i) (liftIndex2 (+) start (liftIndex2 (*) stride (fromLinearIndex sz i)))
+  {-# INLINE iterTargetA_ #-}
+  iterTargetFullAccST it scheduler iStart sz acc splitAcc action =
+    let !(RowMajorLinear fact) = it
+    in iterLinearAccST fact scheduler iStart 1 (totalElem sz) acc splitAcc $ \ !i ->
+      action i (fromLinearIndex sz i)
+  {-# INLINE iterTargetFullAccST #-}
+  iterTargetFullAccST_ it scheduler iStart sz acc splitAcc action =
+    let !(RowMajorLinear fact) = it
+    in iterLinearAccST_ fact scheduler iStart 1 (totalElem sz) acc splitAcc $ \ !i ->
+      action i (fromLinearIndex sz i)
+  {-# INLINE iterTargetFullAccST_ #-}
+  iterTargetFullST_ it scheduler iStart sz action =
+    let !(RowMajorLinear fact) = it
+    in iterLinearST_ fact scheduler iStart 1 (totalElem sz) $ \ !i ->
+      action i (fromLinearIndex sz i)
+  {-# INLINE iterTargetFullST_ #-}
+  iterTargetWithStrideAccST it scheduler iStart sz start (Stride stride) acc spliAcc action =
+    let RowMajorLinear fact = it
+     in iterLinearAccST fact scheduler 0 1 (totalElem sz) acc spliAcc $ \i ->
+          action (iStart + i) $
+            liftIndex2 (+) start (liftIndex2 (*) stride (fromLinearIndex sz i))
+  {-# INLINE iterTargetWithStrideAccST #-}
+  iterTargetWithStrideAccST_ it scheduler iStart sz start (Stride stride) acc spliAcc action =
+    let RowMajorLinear fact = it
+     in iterLinearAccST_ fact scheduler 0 1 (totalElem sz) acc spliAcc $ \i ->
+          action (iStart + i) $
+            liftIndex2 (+) start (liftIndex2 (*) stride (fromLinearIndex sz i))
+  {-# INLINE iterTargetWithStrideAccST_ #-}
+
+
+
+-- | Parallelizing unbalanced computation (i.e. computing some elements of the
+-- array is much more expensive then the others) it can be benefitial to
+-- interleave iteration. Perfect example of this would be a ray tracer or the
+-- Mandelbrot set.
+--
+-- iteration without parallelization is equivalent to `RowMajor`
+--
+-- @since 1.0.2
+newtype RowMajorUnbalanced = RowMajorUnbalancedInternal Int
+
+defRowMajorUnbalanced :: RowMajorUnbalanced
+defRowMajorUnbalanced = RowMajorUnbalancedInternal 8
+
+pattern RowMajorUnbalanced :: Int
+                 -- ^ Multiplier that will be used to scale number of jobs.
+                 -> RowMajorUnbalanced
+pattern RowMajorUnbalanced f <- RowMajorUnbalancedInternal f
+  where RowMajorUnbalanced = RowMajorUnbalancedInternal . max 1
+{-# COMPLETE RowMajorUnbalanced #-}
+
+
+instance Iterator RowMajorUnbalanced where
+  iterFullM (RowMajorUnbalanced fact) = iterFullM (RowMajor fact)
+  {-# INLINE iterFullM #-}
+  iterFullA_ (RowMajorUnbalanced fact) = iterFullA_ (RowMajor fact)
+  {-# INLINE iterFullA_ #-}
+  iterTargetM (RowMajorUnbalanced fact) = iterTargetM (RowMajor fact)
+  {-# INLINE iterTargetM #-}
+  iterTargetA_ (RowMajorUnbalanced fact) = iterTargetA_ (RowMajor fact)
+  {-# INLINE iterTargetA_ #-}
+  iterTargetWithStrideAccST = iterUnbalancedTargetWithStride loopM
+  {-# INLINE iterTargetWithStrideAccST #-}
+  iterTargetWithStrideAccST_ it scheduler iStart sz start stride acc splitAcc' action =
+    void $
+    iterUnbalancedTargetWithStride innerLoop it scheduler iStart sz start stride acc splitAcc' action
+    where
+      innerLoop initial condition increment initAcc f =
+        void $ loopM initial condition increment initAcc f
+      {-# INLINE innerLoop #-}
+  {-# INLINE iterTargetWithStrideAccST_ #-}
+
+iterUnbalancedTargetWithStride ::
+     Index ix
+  => (Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> t) -> ST s b)
+  -> RowMajorUnbalanced
+  -> Scheduler s b
+  -> Int
+  -> Sz ix
+  -> ix
+  -> Stride ix
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Int -> ix -> t)
+  -> ST s a
+iterUnbalancedTargetWithStride innerLoop it scheduler iStart sz start stride acc splitAcc action =
+  let RowMajorUnbalanced fact = it
+      !n = totalElem sz
+      !step = min (fact * numWorkers scheduler) n
+  in loopM 0 (< step) (+ 1) acc $ \ !istep !a -> do
+       (curAcc, nextAcc) <- splitAcc a
+       scheduleMassivWork scheduler $
+         innerLoop istep (< n) (+ step) curAcc $ \i ->
+           action (iStart + i) $
+             liftIndex2 (+) start (liftIndex2 (*) (unStride stride) (fromLinearIndex sz i))
+       pure nextAcc
+{-# INLINE iterUnbalancedTargetWithStride #-}
+
+
+noSplit :: Applicative m => () -> m ((), ())
+noSplit _ = pure ((), ())
+
+
+iterTargetAccST ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s a -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> ix -- ^ Source start
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+  -> ST s a
+iterTargetAccST it scheduler iStart sz ix =
+  iterTargetWithStrideAccST it scheduler iStart sz ix oneStride
+{-# INLINE iterTargetAccST #-}
+
+iterTargetAccST_ ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s () -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> ix -- ^ Source start
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+  -> ST s ()
+iterTargetAccST_ it scheduler iStart sz ix =
+  iterTargetWithStrideAccST_ it scheduler iStart sz ix oneStride
+{-# INLINE iterTargetAccST_ #-}
+
+
+iterTargetFullWithStrideST_ ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s () -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> Stride ix -- ^ Stride
+  -> (Ix1 -> ix -> ST s ()) -- ^ Action
+  -> ST s ()
+iterTargetFullWithStrideST_ it scheduler iStart sz =
+  iterTargetWithStrideST_ it scheduler iStart sz (pureIndex 0)
+{-# INLINE iterTargetFullWithStrideST_ #-}
+
+iterTargetST_ ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s () -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> ix -- ^ Start
+  -> (Ix1 -> ix -> ST s ()) -- ^ Action
+  -> ST s ()
+iterTargetST_ it scheduler iStart sz ix =
+  iterTargetWithStrideST_ it scheduler iStart sz ix oneStride
+{-# INLINE iterTargetST_ #-}
+
+
+iterTargetFullWithStrideAccST ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s a -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> Stride ix -- ^ Stride
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+  -> ST s a
+iterTargetFullWithStrideAccST it scheduler iStart sz =
+  iterTargetWithStrideAccST it scheduler iStart sz (pureIndex 0)
+{-# INLINE iterTargetFullWithStrideAccST #-}
+
+iterTargetFullWithStrideAccST_ ::
+     (Iterator it, Index ix)
+  => it
+  -> Scheduler s () -- ^ Scheduler to use
+  -> Ix1 -- ^ Target linear start index
+  -> Sz ix -- ^ Target size
+  -> Stride ix -- ^ Stride
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Ix1 -> ix -> a -> ST s a) -- ^ Action
+  -> ST s ()
+iterTargetFullWithStrideAccST_ it scheduler iStart sz =
+  iterTargetWithStrideAccST_ it scheduler iStart sz (pureIndex 0)
+{-# INLINE iterTargetFullWithStrideAccST_ #-}
diff --git a/src/Data/Massiv/Core/Index/Ix.hs b/src/Data/Massiv/Core/Index/Ix.hs
--- a/src/Data/Massiv/Core/Index/Ix.hs
+++ b/src/Data/Massiv/Core/Index/Ix.hs
@@ -7,12 +7,12 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE PatternSynonyms #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE TypeFamilyDependencies #-}
+{-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Ix
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -40,16 +40,17 @@
   , HighIxN
   ) where
 
-import Control.Monad.Catch (MonadThrow(..))
 import Control.DeepSeq
+import Control.Monad.Catch (MonadThrow(..))
 import Data.Massiv.Core.Index.Internal
 import Data.Proxy
-import qualified GHC.Arr as I
 import qualified Data.Vector.Generic as V
 import qualified Data.Vector.Generic.Mutable as VM
 import qualified Data.Vector.Unboxed as VU
+import qualified GHC.Arr as I
 import GHC.TypeLits
 import System.Random.Stateful
+import Data.Massiv.Core.Loop
 #if !MIN_VERSION_base(4,11,0)
 import Data.Semigroup
 #endif
@@ -310,8 +311,9 @@
   {-# INLINE [1] isSafeIndex #-}
   toLinearIndex (SafeSz (_ :. k1)) (i2 :. i1) = k1 * i2 + i1
   {-# INLINE [1] toLinearIndex #-}
-  fromLinearIndex (SafeSz (_ :. k1)) i = case i `quotRem` k1 of
-                                           (i2, i1) -> i2 :. i1
+  fromLinearIndex (SafeSz (_ :. k1)) i =
+    case i `quotRem` k1 of
+      (i2, i1) -> i2 :. i1
   {-# INLINE [1] fromLinearIndex #-}
   consDim = (:.)
   {-# INLINE [1] consDim #-}
@@ -321,21 +323,21 @@
   {-# INLINE [1] snocDim #-}
   unsnocDim (i2 :. i1) = (i2, i1)
   {-# INLINE [1] unsnocDim #-}
-  getDimM (i2 :.  _) 2 = pure i2
-  getDimM ( _ :. i1) 1 = pure i1
-  getDimM ix         d = throwM $ IndexDimensionException ix d
+  getDimM (i2 :. _) 2 = pure i2
+  getDimM (_ :. i1) 1 = pure i1
+  getDimM ix d        = throwM $ IndexDimensionException ix d
   {-# INLINE [1] getDimM #-}
-  setDimM ( _ :. i1) 2 i2 = pure (i2 :. i1)
-  setDimM (i2 :.  _) 1 i1 = pure (i2 :. i1)
-  setDimM ix         d _  = throwM $ IndexDimensionException ix d
+  setDimM (_ :. i1) 2 i2 = pure (i2 :. i1)
+  setDimM (i2 :. _) 1 i1 = pure (i2 :. i1)
+  setDimM ix d _         = throwM $ IndexDimensionException ix d
   {-# INLINE [1] setDimM #-}
   pullOutDimM (i2 :. i1) 2 = pure (i2, i1)
   pullOutDimM (i2 :. i1) 1 = pure (i1, i2)
-  pullOutDimM ix         d = throwM $ IndexDimensionException ix d
+  pullOutDimM ix d         = throwM $ IndexDimensionException ix d
   {-# INLINE [1] pullOutDimM #-}
   insertDimM i1 2 i2 = pure (i2 :. i1)
   insertDimM i2 1 i1 = pure (i2 :. i1)
-  insertDimM ix d  _ = throwM $ IndexDimensionException ix d
+  insertDimM ix d _  = throwM $ IndexDimensionException ix d
   {-# INLINE [1] insertDimM #-}
   pureIndex i = i :. i
   {-# INLINE [1] pureIndex #-}
@@ -346,6 +348,10 @@
   repairIndex (SafeSz (k :. szL)) (i :. ixL) rBelow rOver =
     repairIndex (SafeSz k) i rBelow rOver :. repairIndex (SafeSz szL) ixL rBelow rOver
   {-# INLINE [1] repairIndex #-}
+  iterF (s :. sIxL) (e :. eIxL) (inc :. incIxL) cond initAct f =
+    loopF s (`cond` e) (+ inc) initAct $ \ !i g ->
+      loopF sIxL (`cond` eIxL) (+ incIxL) g $ \ !j -> f (i :. j)
+  {-# INLINE iterF #-}
 
 
 instance {-# OVERLAPPING #-} Index (IxN 3) where
@@ -357,10 +363,15 @@
   isSafeIndex (SafeSz (k3 :> k2 :. k1)) (i3 :> i2 :. i1) =
     0 <= i3 && 0 <= i2 && 0 <= i1 && i3 < k3 && i2 < k2 && i1 < k1
   {-# INLINE [1] isSafeIndex #-}
-  toLinearIndex (SafeSz (_ :> k2 :. k1)) (i3 :> i2 :. i1) = (k2 * i3 + i2) * k1 + i1
+  toLinearIndex (SafeSz (_ :> k2 :. k1)) (i3 :> i2 :. i1) = (i3 * k2 + i2) * k1 + i1
   {-# INLINE [1] toLinearIndex #-}
   fromLinearIndex (SafeSz (_ :> ix)) i = let !(q, ixL) = fromLinearIndexAcc ix i in q :> ixL
-  {-# INLINE [1] fromLinearIndex #-}
+  {-# INLINE fromLinearIndex #-}
+  fromLinearIndexAcc (m :> ix) !k = (q, r :> ixL)
+    where
+      !(!kL, !ixL) = fromLinearIndexAcc ix k
+      !(!q, !r) = quotRem kL m
+  {-# INLINE fromLinearIndexAcc #-}
   consDim = (:>)
   {-# INLINE [1] consDim #-}
   unconsDim (i3 :> ix) = (i3, ix)
@@ -398,6 +409,27 @@
   repairIndex (SafeSz (n :> szL)) (i :> ixL) rBelow rOver =
     repairIndex (SafeSz n) i rBelow rOver :> repairIndex (SafeSz szL) ixL rBelow rOver
   {-# INLINE [1] repairIndex #-}
+  iterTargetRowMajorAccM iAcc iStart sz (b3 :> b2 :. b1) (s3 :> s2 :. s1) initAcc action =
+    let n = totalElem sz
+        iShift = iStart + iAcc * n
+     in loopM 0 (< n) (+ 1) initAcc $ \ !i !acc ->
+          let (i3 :> i2 :. i1) = fromLinearIndex sz i
+          in action (iShift + i) ((b3 + s3 * i3) :> (b2 + s2 * i2) :. (b1 + s1 * i1)) acc
+  {-# INLINE iterTargetRowMajorAccM #-}
+  iterTargetRowMajorAccST_ iAcc fact scheduler iStart sz (b3 :> b2 :. b1) (s3 :> s2 :. s1) acc splitAcc action =
+    let n = totalElem sz
+        iShift = iStart + iAcc * n
+     in iterLinearAccST_ fact scheduler 0 1 n acc splitAcc $ \ !i ->
+          let (i3 :> i2 :. i1) = fromLinearIndex sz i
+          in action (iShift + i) ((b3 + s3 * i3) :> (b2 + s2 * i2) :. (b1 + s1 * i1))
+  {-# INLINE iterTargetRowMajorAccST_ #-}
+  iterTargetRowMajorAccST iAcc fact scheduler iStart sz (b3 :> b2 :. b1) (s3 :> s2 :. s1) acc splitAcc action =
+    let n = totalElem sz
+        iShift = iStart + iAcc * n
+     in iterLinearAccST fact scheduler 0 1 n acc splitAcc $ \ !i ->
+          let (i3 :> i2 :. i1) = fromLinearIndex sz i
+          in action (iShift + i) ((b3 + s3 * i3) :> (b2 + s2 * i2) :. (b1 + s1 * i1))
+  {-# INLINE iterTargetRowMajorAccST #-}
 
 -- | Constraint synonym that encapsulates all constraints needed for dimension 4 and higher.
 --
diff --git a/src/Data/Massiv/Core/Index/Stride.hs b/src/Data/Massiv/Core/Index/Stride.hs
--- a/src/Data/Massiv/Core/Index/Stride.hs
+++ b/src/Data/Massiv/Core/Index/Stride.hs
@@ -3,7 +3,7 @@
 {-# LANGUAGE PatternSynonyms #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Stride
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -19,9 +19,9 @@
   , strideSize
   ) where
 
-import Control.DeepSeq
+import Control.DeepSeq (NFData)
 import Data.Massiv.Core.Index.Internal
-import System.Random.Stateful
+import System.Random.Stateful (Random, Uniform(..), UniformRange(..))
 
 -- | Stride provides a way to ignore elements of an array if an index is divisible by a
 -- corresponding value in a stride. So, for a @Stride (i :. j)@ only elements with indices will be
diff --git a/src/Data/Massiv/Core/Index/Tuple.hs b/src/Data/Massiv/Core/Index/Tuple.hs
--- a/src/Data/Massiv/Core/Index/Tuple.hs
+++ b/src/Data/Massiv/Core/Index/Tuple.hs
@@ -6,7 +6,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Tuple
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Core/Iterator.hs b/src/Data/Massiv/Core/Iterator.hs
deleted file mode 100644
--- a/src/Data/Massiv/Core/Iterator.hs
+++ /dev/null
@@ -1,204 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE MonoLocalBinds #-}
--- |
--- Module      : Data.Massiv.Core.Iterator
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
--- License     : BSD3
--- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
--- Stability   : experimental
--- Portability : non-portable
---
-module Data.Massiv.Core.Iterator
-  ( loop
-  , loopA_
-  , loopM
-  , loopM_
-  , loopDeepM
-  , splitLinearly
-  , splitLinearlyM_
-  , splitLinearlyWith_
-  , splitLinearlyWithM_
-  , splitLinearlyWithStartAtM_
-  , splitLinearlyWithStatefulM_
-  ) where
-
-import Control.Scheduler
-import Control.Monad
-import Control.Monad.Primitive
-import Control.Monad.IO.Unlift
-
--- | Efficient loop with an accumulator
---
--- @since 0.1.0
-loop :: Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
-loop !init' condition increment !initAcc f = go init' initAcc
-  where
-    go !step !acc
-      | condition step = go (increment step) (f step acc)
-      | otherwise = acc
-{-# INLINE loop #-}
-
-
--- | Efficient monadic loop with an accumulator
---
--- >>> loopM 1 (< 20) (+ 2) [] (\i a -> Just (i:a))
--- Just [19,17,15,13,11,9,7,5,3,1]
---
--- @since 0.1.0
-loopM :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> m a) -> m a
-loopM !init' condition increment !initAcc f = go init' initAcc
-  where
-    go !step !acc
-      | condition step = f step acc >>= go (increment step)
-      | otherwise = return acc
-{-# INLINE loopM #-}
-
-
--- | Efficient monadic loop. Result of each iteration is discarded.
---
--- @since 0.1.0
-loopM_ :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> m a) -> m ()
-loopM_ !init' condition increment f = go init'
-  where
-    go !step
-      | condition step = f step >> go (increment step)
-      | otherwise = pure ()
-{-# INLINE loopM_ #-}
-
--- | Similar to `loopM_` axcept the action accepts not only the value for current step,
--- but also for the next one as well.
---
--- @since 0.5.7
-loopNextM_ :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> Int -> m a) -> m ()
-loopNextM_ !init' condition increment f = go init'
-  where
-    go step =
-      when (condition step) $
-      let !next = increment step
-       in f step next >> go next
-{-# INLINE loopNextM_ #-}
-
-
--- | Efficient Applicative loop. Result of each iteration is discarded.
---
--- @since 0.3.0
-loopA_ :: Applicative f => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> f a) -> f ()
-loopA_ !init' condition increment f = go init'
-  where
-    go !step
-      | condition step = f step *> go (increment step)
-      | otherwise = pure ()
-{-# INLINE loopA_ #-}
-
-
--- | Similar to `loopM`, but slightly less efficient monadic loop with an accumulator that reverses
--- the direction of action application. eg:
---
--- >>> loopDeepM 1 (< 20) (+ 2) [] (\i a -> Just (i:a))
--- Just [1,3,5,7,9,11,13,15,17,19]
---
--- Equivalent to:
---
--- >>> loopM 19 (>= 1) (subtract 2) [] (\i a -> Just (i:a))
--- Just [1,3,5,7,9,11,13,15,17,19]
---
--- @since 0.1.0
-loopDeepM :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> m a) -> m a
-loopDeepM !init' condition increment !initAcc f = go init' initAcc
-  where
-    go !step !acc
-      | condition step = go (increment step) acc >>= f step
-      | otherwise = return acc
-{-# INLINE loopDeepM #-}
-
-
--- | Divide length in chunks and apply a function to the computed results
---
--- @since 0.2.1
-splitLinearly :: Int -- ^ Number of chunks
-              -> Int -- ^ Total length
-              -> (Int -> Int -> a) -- ^ Function that accepts a chunk length and slack start index
-              -> a
-splitLinearly numChunks totalLength action = action chunkLength slackStart
-  where
-    !chunkLength = totalLength `quot` numChunks
-    !slackStart = chunkLength * numChunks
-{-# INLINE splitLinearly #-}
-
--- | Iterator that expects an action that accepts starting linear index as well as the ending
---
--- @since 0.5.7
-splitLinearlyM_ ::
-     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> Int -> m ()) -> m ()
-splitLinearlyM_ scheduler totalLength action =
-  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-    loopNextM_ 0 (< slackStart) (+ chunkLength) $ \ start next ->
-      scheduleWork_ scheduler $ action start next
-    when (slackStart < totalLength) $
-      scheduleWork_ scheduler $ action slackStart totalLength
-{-# INLINE splitLinearlyM_ #-}
-
--- | Interator that can be used to split computation amongst different workers. For monadic
--- generator see `splitLinearlyWithM_`.
---
--- @since 0.2.1
-splitLinearlyWith_ ::
-     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> b) -> (Int -> b -> m ()) -> m ()
-splitLinearlyWith_ scheduler totalLength index =
-  splitLinearlyWithM_ scheduler totalLength (pure . index)
-{-# INLINE splitLinearlyWith_ #-}
-
-
--- | Interator that can be used to split computation jobs
---
--- @since 0.2.6
-splitLinearlyWithM_ ::
-     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m ()
-splitLinearlyWithM_ scheduler totalLength make write =
-  splitLinearlyM_ scheduler totalLength go
-  where
-    go start end = loopM_ start (< end) (+ 1) $ \ k -> make k >>= write k
-    {-# INLINE go #-}
-{-# INLINE splitLinearlyWithM_ #-}
-
-
--- | Interator that can be used to split computation jobs
---
--- @since 0.3.0
-splitLinearlyWithStartAtM_ ::
-     MonadPrimBase s m => Scheduler s () -> Int -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m ()
-splitLinearlyWithStartAtM_ scheduler startAt totalLength make write =
-  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-    loopM_ startAt (< (slackStart + startAt)) (+ chunkLength) $ \ !start ->
-      scheduleWork_ scheduler $
-      loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> make k >>= write k
-    when (slackStart < totalLength) $
-      scheduleWork_ scheduler $
-        loopM_ (slackStart + startAt) (< (totalLength + startAt)) (+ 1) $ \ !k -> make k >>= write k
-{-# INLINE splitLinearlyWithStartAtM_ #-}
-
-
-
--- | Interator that can be used to split computation jobs, while using a stateful scheduler.
---
--- @since 0.3.4
-splitLinearlyWithStatefulM_ ::
-     MonadUnliftIO m
-  => SchedulerWS ws ()
-  -> Int -- ^ Total linear length
-  -> (Int -> ws -> m b) -- ^ Element producing action
-  -> (Int -> b -> m c) -- ^ Element storing action
-  -> m ()
-splitLinearlyWithStatefulM_ schedulerWS totalLength make store =
-  let nWorkers = numWorkers (unwrapSchedulerWS schedulerWS)
-   in withRunInIO $ \run ->
-      splitLinearly nWorkers totalLength $ \chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
-          scheduleWorkState_ schedulerWS $ \s ->
-            loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k ->
-              run (make k s >>= store k)
-        scheduleWorkState_ schedulerWS $ \s ->
-          loopM_ slackStart (< totalLength) (+ 1) $ \ !k ->
-            run (make k s >>= store k)
-{-# INLINE splitLinearlyWithStatefulM_ #-}
diff --git a/src/Data/Massiv/Core/List.hs b/src/Data/Massiv/Core/List.hs
--- a/src/Data/Massiv/Core/List.hs
+++ b/src/Data/Massiv/Core/List.hs
@@ -4,8 +4,6 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE LambdaCase #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE RankNTypes #-}
-{-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
@@ -13,7 +11,7 @@
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Massiv.Core.List
--- Copyright   : (c) Alexey Kuleshevich 2018-2021
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -123,6 +121,7 @@
 
 outerLength :: Array L ix e -> Sz Int
 outerLength = SafeSz . length . unList . lData
+{-# INLINE outerLength #-}
 
 
 instance Ragged L Ix1 e where
@@ -227,7 +226,7 @@
         !k = unSz ksz
     withScheduler comp $ \ scheduler ->
       splitLinearly (numWorkers scheduler) k $ \ chunkLength slackStart -> do
-        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
           scheduleWork scheduler $ do
             res <- loopDeepM start (< (start + chunkLength)) (+ 1) [] $ \i acc ->
               return (fmap lData (generateRaggedM Seq szL (\ !ixL -> f (consDim i ixL))):acc)
@@ -245,6 +244,7 @@
   {-# INLINE setComp #-}
   getComp = lComp
   {-# INLINE getComp #-}
+  repr = L
 
 -- -- TODO: benchmark against using unsafeGenerateM directly
 -- unsafeGenerateN ::
diff --git a/src/Data/Massiv/Core/Loop.hs b/src/Data/Massiv/Core/Loop.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Massiv/Core/Loop.hs
@@ -0,0 +1,455 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MonoLocalBinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-- |
+-- Module      : Data.Massiv.Core.Loop
+-- Copyright   : (c) Alexey Kuleshevich 2018-2022
+-- License     : BSD3
+-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
+-- Stability   : experimental
+-- Portability : non-portable
+--
+module Data.Massiv.Core.Loop
+  ( loop
+  , loopF
+  , nextMaybeF
+  , loopA
+  , loopA_
+  , loopM
+  , loopM_
+  , iloopM
+  , iloopA_
+  , loopNextM
+  , loopNextA_
+  , loopDeepM
+  , splitLinearly
+  , splitLinearlyM
+  , splitLinearlyM_
+  , splitLinearlyWith_
+  , splitLinearlyWithM_
+  , splitLinearlyWithStartAtM_
+  , splitLinearlyWithStatefulM_
+  , iterLinearST_
+  , iterLinearAccST_
+  , iterLinearAccST
+  , splitNumChunks
+  , stepStartAdjust
+  -- * Experimental
+  , splitWorkWithFactorST
+  , scheduleMassivWork
+  , withMassivScheduler_
+  ) where
+
+import Control.Monad (void, when)
+import Control.Monad.IO.Unlift (MonadUnliftIO(..))
+import Control.Monad.Primitive
+import Control.Monad.ST (ST)
+import Control.Scheduler (Comp(..), Scheduler, SchedulerWS,
+                          numWorkers, scheduleWork, scheduleWorkState_,
+                          scheduleWork_, trivialScheduler_, unwrapSchedulerWS,
+                          withScheduler_)
+import Control.Scheduler.Global (globalScheduler, withGlobalScheduler_)
+import Data.Coerce
+import Data.Functor.Identity
+
+-- | Efficient loop with an accumulator
+--
+-- @since 0.1.0
+loop :: Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
+loop initial condition increment initAcc f =
+  runIdentity (loopM initial condition increment initAcc (coerce f))
+{-# INLINE loop #-}
+
+-- | Efficient monadic loop with an accumulator
+--
+-- >>> loopM 1 (< 20) (+ 2) [] (\i a -> Just (i:a))
+-- Just [19,17,15,13,11,9,7,5,3,1]
+--
+-- @since 0.1.0
+loopM :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> m a) -> m a
+loopM !initial condition increment !initAcc f =
+  go initial initAcc
+  where
+    go !step !acc
+      | condition step = f step acc >>= go (increment step)
+      | otherwise = pure acc
+{-# INLINE loopM #-}
+
+
+-- | Efficient monadic loop with an accumulator and extra linear index incremented by 1.
+--
+-- >>> iloopM 100 1 (< 20) (+ 2) [] (\i ix a -> Just ((i, ix) : a))
+-- Just [(109,19),(108,17),(107,15),(106,13),(105,11),(104,9),(103,7),(102,5),(101,3),(100,1)]
+--
+-- @since 1.0.2
+iloopM ::
+     Monad m => Int -> Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> Int -> a -> m a) -> m a
+iloopM !istart !initIx condition increment !initAcc f = go istart initIx initAcc
+  where
+    go !i !step !acc
+      | condition step = f i step acc >>= go (i + 1) (increment step)
+      | otherwise = pure acc
+{-# INLINE iloopM #-}
+
+-- | Efficient monadic loop. Result of each iteration is discarded.
+--
+-- @since 0.1.0
+loopM_ :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> m a) -> m ()
+loopM_ !initial condition increment f = go initial
+  where
+    go !step
+      | condition step = f step >> go (increment step)
+      | otherwise = pure ()
+  --loopF initial condition increment (pure ()) (\i ma -> f i >> ma)
+{-# INLINE loopM_ #-}
+{-# DEPRECATED loopM_ "In favor of `loopA_`" #-}
+
+-- | Efficient monadic loop with extra linear index incremented by 1.
+--
+-- >>> iloopA_ 100 1 (< 10) (+ 2) (\i ix -> print (i, ix))
+-- (100,1)
+-- (101,3)
+-- (102,5)
+-- (103,7)
+-- (104,9)
+--
+-- @since 1.0.2
+iloopA_ ::
+     Applicative f => Int -> Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> Int -> f a) -> f ()
+iloopA_ !istart !initIx condition increment f = go istart initIx
+  where
+    go !i !step
+      | condition step = f i step *> go (i + 1) (increment step)
+      | otherwise = pure ()
+{-# INLINE iloopA_ #-}
+
+-- | Similar to `loopM_` except the action accepts not only the value for current step,
+-- but also for the next one as well.
+--
+-- @since 1.0.2
+loopNextA_ :: Applicative f => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> Int -> f a) -> f ()
+loopNextA_ !initial condition increment f = go initial
+  where
+    go !step
+      | condition step =
+        let !next = increment step
+        in f step next *> go next
+      | otherwise = pure ()
+{-# INLINE loopNextA_ #-}
+
+-- | Similar to `loopM_` except the action accepts not only the value for current step,
+-- but also for the next one as well.
+--
+-- @since 1.0.2
+loopNextM :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> Int -> a -> m a) -> m a
+loopNextM !initial condition increment !initAcc f = go initial initAcc
+  where
+    go !step !acc
+      | condition step =
+        let !next = increment step
+        in f step next acc >>= go next
+      | otherwise = pure acc
+{-# INLINE loopNextM #-}
+
+-- | Efficient Applicative loop. Result of each iteration is discarded.
+--
+-- > loopA_ initial cond incr f === loopA initial cond incr (pure ()) (\i -> id <$ f i)
+--
+-- @since 1.0.2
+loopA_ :: Applicative f => Int -> (Int -> Bool) -> (Int -> Int) -> (Int -> f a) -> f ()
+loopA_ !initial condition increment f =
+  loopF initial condition increment (pure ()) (\i ma -> f i *> ma)
+{-# INLINE loopA_ #-}
+
+-- | Applicative loop. Use monadic `loopM` when possible, since it will be more efficient.
+--
+-- @since 0.3.0
+loopA :: Applicative f => Int -> (Int -> Bool) -> (Int -> Int) -> f b -> (Int -> f (b -> b)) -> f b
+loopA !initial condition increment lastAction f =
+  loopF initial condition increment lastAction (\i ma -> f i <*> ma)
+{-# INLINE loopA #-}
+
+
+loopF :: Int -> (Int -> Bool) -> (Int -> Int) -> f a -> (Int -> f a -> f a) -> f a
+loopF !initial condition increment lastAction f = go initial
+  where
+    go !step
+      | condition step = f step (go (increment step))
+      | otherwise = lastAction
+{-# INLINE loopF #-}
+
+
+nextMaybeF :: Int -> (Int -> Bool) -> (Int -> Int) -> (Maybe Int -> f a) -> f a
+nextMaybeF !cur condition increment f =
+  let !i = increment cur
+  in f $! if condition i then Just i else Nothing
+{-# INLINE nextMaybeF #-}
+
+-- | Similar to `loopM`, but way less efficient monadic loop with an accumulator that reverses
+-- the direction of action application. eg:
+--
+-- >>> loopDeepM 1 (< 20) (+ 2) [] (\i a -> Just (i:a))
+-- Just [1,3,5,7,9,11,13,15,17,19]
+--
+-- Equivalent to:
+--
+-- >>> loopM 19 (>= 1) (subtract 2) [] (\i a -> Just (i:a))
+-- Just [1,3,5,7,9,11,13,15,17,19]
+--
+-- @since 0.1.0
+loopDeepM :: Monad m => Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> m a) -> m a
+loopDeepM !initial condition increment !initAcc f =
+  loopF initial condition increment (pure initAcc) (\i ma -> ma >>= f i)
+{-# INLINE loopDeepM #-}
+
+
+-- | Divide length in chunks and apply a function to the computed results
+--
+-- @since 0.2.1
+splitLinearly :: Int -- ^ Number of chunks
+              -> Int -- ^ Total length
+              -> (Int -> Int -> a) -- ^ Function that accepts a chunk length and slack start index
+              -> a
+splitLinearly numChunks totalLength action = action chunkLength slackStart
+  where
+    !chunkLength = totalLength `quot` numChunks
+    !slackStart = chunkLength * numChunks
+{-# INLINE splitLinearly #-}
+
+-- | Iterator that expects an action that accepts starting linear index as well as the ending
+--
+-- @since 0.5.7
+splitLinearlyM_ ::
+     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> Int -> m ()) -> m ()
+splitLinearlyM_ scheduler totalLength action =
+  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+    loopNextA_ 0 (< slackStart) (+ chunkLength) $ \ start next ->
+      scheduleWork_ scheduler $ action start next
+    when (slackStart < totalLength) $
+      scheduleWork_ scheduler $ action slackStart totalLength
+{-# INLINE splitLinearlyM_ #-}
+
+-- | Iterator that expects an action that accepts starting linear index as well as the ending
+--
+-- @since 1.0.2
+splitLinearlyM ::
+     MonadPrimBase s m => Scheduler s a -> Int -> (Int -> Int -> m a) -> m ()
+splitLinearlyM scheduler totalLength action =
+  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+    loopNextA_ 0 (< slackStart) (+ chunkLength) $ \ start next ->
+      scheduleWork scheduler (action start next)
+    when (slackStart < totalLength) $
+      scheduleWork scheduler (action slackStart totalLength)
+{-# INLINE splitLinearlyM #-}
+
+-- | Iterator that can be used to split computation amongst different workers. For monadic
+-- generator see `splitLinearlyWithM_`.
+--
+-- @since 0.2.1
+splitLinearlyWith_ ::
+     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> b) -> (Int -> b -> m ()) -> m ()
+splitLinearlyWith_ scheduler totalLength index =
+  splitLinearlyWithM_ scheduler totalLength (pure . index)
+{-# INLINE splitLinearlyWith_ #-}
+
+
+-- | Iterator that can be used to split computation jobs
+--
+-- @since 0.2.6
+splitLinearlyWithM_ ::
+     MonadPrimBase s m => Scheduler s () -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m ()
+splitLinearlyWithM_ scheduler totalLength make write =
+  splitLinearlyM_ scheduler totalLength go
+  where
+    go start end = loopM_ start (< end) (+ 1) $ \ k -> make k >>= write k
+    {-# INLINE go #-}
+{-# INLINE splitLinearlyWithM_ #-}
+
+
+-- | Iterator that can be used to split computation jobs
+--
+-- @since 0.3.0
+splitLinearlyWithStartAtM_ ::
+     MonadPrimBase s m => Scheduler s () -> Int -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m ()
+splitLinearlyWithStartAtM_ scheduler startAt totalLength make write =
+  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+    loopM_ startAt (< (slackStart + startAt)) (+ chunkLength) $ \ !start ->
+      scheduleWork_ scheduler $
+      loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> make k >>= write k
+    when (slackStart < totalLength) $
+      scheduleWork_ scheduler $
+        loopM_ (slackStart + startAt) (< (totalLength + startAt)) (+ 1) $ \ !k -> make k >>= write k
+{-# INLINE splitLinearlyWithStartAtM_ #-}
+
+
+
+-- | Iterator that can be used to split computation jobs, while using a stateful scheduler.
+--
+-- @since 0.3.4
+splitLinearlyWithStatefulM_ ::
+     MonadUnliftIO m
+  => SchedulerWS ws ()
+  -> Int -- ^ Total linear length
+  -> (Int -> ws -> m b) -- ^ Element producing action
+  -> (Int -> b -> m c) -- ^ Element storing action
+  -> m ()
+splitLinearlyWithStatefulM_ schedulerWS totalLength make store =
+  let nWorkers = numWorkers (unwrapSchedulerWS schedulerWS)
+   in withRunInIO $ \run ->
+      splitLinearly nWorkers totalLength $ \chunkLength slackStart -> do
+        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+          scheduleWorkState_ schedulerWS $ \s ->
+            loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k ->
+              run (make k s >>= store k)
+        scheduleWorkState_ schedulerWS $ \s ->
+          loopM_ slackStart (< totalLength) (+ 1) $ \ !k ->
+            run (make k s >>= store k)
+{-# INLINE splitLinearlyWithStatefulM_ #-}
+
+
+-- | This is a major helper function for fair splitting and parallelization of
+-- work with ability to use some arbitrary accumulator and splittable seed
+--
+-- @since 1.0.2
+splitWorkWithFactorST ::
+     Int -- ^ Multiplying factor to be applied to number of workers for number
+         -- of jobs to schedule. Higher the factor, more jobs will be
+         -- scheduled. Only positive values are valid.
+  -> Scheduler s a
+  -> Int -- ^ Starting index
+  -> Int -- ^ Stepping value. Can be negative, but must not be zero.
+  -> Int -- ^ Total number of steps to be taken
+  -> b -- ^ Initial value for an accumulator
+  -> (b -> ST s (b, b)) -- ^ An action to split accumulator for multiple threads
+  -> (Int -> Int -> Int -> Int -> b -> ST s a)
+  -- ^ A job to be scheduled. Accepts:
+  --
+  -- * Chunk index start
+  -- * Chunk length
+  -- * Chunk start index adjusted for supplied start and stepping value
+  -- * Chunk stop index adjusted for supplied start and stepping value
+  -> ST s b
+splitWorkWithFactorST fact scheduler start step totalLength initAcc splitAcc action = do
+  let !(chunkLength, slackStart) = splitNumChunks fact (numWorkers scheduler) totalLength
+  slackAcc <-
+    loopM 0 (< slackStart) (+ chunkLength) initAcc $ \ !chunkStart !acc -> do
+      (accCur, accNext) <- splitAcc acc
+      scheduleMassivWork scheduler $ do
+        let !chunkStartAdj = start + chunkStart * step
+            !chunkStopAdj = chunkStartAdj + chunkLength * step
+        action chunkStart chunkLength chunkStartAdj chunkStopAdj accCur
+      pure accNext
+  let !slackLength = totalLength - slackStart
+  if slackLength > 0
+    then do
+      (curAcc, nextAcc) <- splitAcc slackAcc
+      scheduleMassivWork scheduler $ do
+        let !slackStartAdj = start + slackStart * step
+            !slackStopAdj = slackStartAdj + slackLength * step
+        action slackStart slackLength slackStartAdj slackStopAdj curAcc
+      pure nextAcc
+    else pure slackAcc
+{-# INLINE splitWorkWithFactorST #-}
+
+-- | Linear iterator that supports multiplying factor
+--
+-- @since 1.0.2
+iterLinearST_ ::
+     Int
+  -> Scheduler s ()
+  -> Int
+  -> Int
+  -> Int
+  -> (Int -> ST s a)
+  -> ST s ()
+iterLinearST_ fact scheduler start step n action = do
+  let totalLength = (n - start) `quot` step
+  splitWorkWithFactorST fact scheduler start step totalLength () (\_ -> pure ((), ()))
+    $ \ _ _ chunkStartAdj chunkStopAdj _ ->
+    loopA_ chunkStartAdj (< chunkStopAdj) (+ step) action
+{-# INLINE iterLinearST_ #-}
+
+-- | Linear iterator that supports multiplying factor and accumulator, but the results are discarded.
+--
+-- @since 1.0.2
+iterLinearAccST_ ::
+     Int
+  -> Scheduler s ()
+  -> Int
+  -> Int
+  -> Int
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Int -> a -> ST s a)
+  -> ST s ()
+iterLinearAccST_ fact scheduler start step n initAcc splitAcc action = do
+  let totalLength = (n - start) `quot` step
+  void $ splitWorkWithFactorST fact scheduler start step totalLength initAcc splitAcc
+    $ \ _ _ chunkStartAdj chunkStopAdj accCur ->
+    void $ loopM chunkStartAdj (< chunkStopAdj) (+ step) accCur action
+{-# INLINE iterLinearAccST_ #-}
+
+-- | Linear iterator that supports multiplying factor and accumulator. Results
+-- of actions are stored in the scheduler.
+--
+-- @since 1.0.2
+iterLinearAccST ::
+     Int
+  -> Scheduler s a
+  -> Int
+  -> Int -- ^ Step. Must be non-zero
+  -> Int
+  -> a
+  -> (a -> ST s (a, a))
+  -> (Int -> a -> ST s a)
+  -> ST s a
+iterLinearAccST fact scheduler start step n initAcc splitAcc action = do
+  let totalLength = (n - start) `quot` step
+  splitWorkWithFactorST fact scheduler start step totalLength initAcc splitAcc
+    $ \ _ _ chunkStartAdj chunkStopAdj accCur ->
+    loopM chunkStartAdj (< chunkStopAdj) (+ step) accCur action
+{-# INLINE iterLinearAccST #-}
+
+
+-- | Helper for figuring out the chunk length and slack start
+splitNumChunks :: Int -> Int -> Int -> (Int, Int)
+splitNumChunks fact nw totalLength =
+  let maxNumChunks = nw * max 1 fact
+      !numChunks
+        | nw == 1 || totalLength <= 0 = 1 -- Optimize for Seq and avoid `quot` by 0.
+        | totalLength <= nw = totalLength
+        | totalLength >= maxNumChunks = maxNumChunks
+        | otherwise = nw
+      !chunkLength = totalLength `quot` numChunks
+      !slackStart = chunkLength * numChunks
+  in (chunkLength, slackStart)
+
+
+-- | Helper for adjusting stride of a chunk
+stepStartAdjust :: Int -> Int -> Int
+stepStartAdjust step ix = ix + ((step - (ix `mod` step)) `mod` step)
+{-# INLINE stepStartAdjust #-}
+
+
+-- | Internal version of a `scheduleWork` that will be replaced by
+-- `scheduleWork_` by the compiler whenever action produces `()`
+scheduleMassivWork :: PrimBase m => Scheduler (PrimState m) a -> m a -> m ()
+scheduleMassivWork = scheduleWork
+{-# INLINE[0] scheduleMassivWork #-}
+
+{-# RULES
+"scheduleWork/scheduleWork_/ST" forall (scheduler :: Scheduler s ()) (action :: ST s ()) . scheduleMassivWork scheduler action = scheduleWork_ scheduler action
+"scheduleWork/scheduleWork_/IO" forall (scheduler :: Scheduler RealWorld ()) (action :: IO ()) . scheduleMassivWork scheduler action = scheduleWork_ scheduler action
+ #-}
+
+-- | Selects an optimal scheduler for the supplied strategy, but it works only in `IO`
+--
+-- @since 1.0.0
+withMassivScheduler_ :: Comp -> (Scheduler RealWorld () -> IO ()) -> IO ()
+withMassivScheduler_ comp f =
+  case comp of
+    Par -> withGlobalScheduler_ globalScheduler f
+    Seq -> f trivialScheduler_
+    _   -> withScheduler_ comp f
+{-# INLINE withMassivScheduler_ #-}
diff --git a/src/Data/Massiv/Core/Operations.hs b/src/Data/Massiv/Core/Operations.hs
--- a/src/Data/Massiv/Core/Operations.hs
+++ b/src/Data/Massiv/Core/Operations.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Core.Operations
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -149,9 +149,7 @@
   -> Array r ix e
   -> Array r ix e
 defaultUnsafeLiftArray2 f a1 a2 =
-  makeArrayLinear
-    (getComp a1 <> getComp a2)
-    (SafeSz (liftIndex2 min (unSz (size a1)) (unSz (size a2)))) $ \ !i ->
+  makeArrayLinear (getComp a1 <> getComp a2) (size a1) $ \ !i ->
     f (unsafeLinearIndex a1 i) (unsafeLinearIndex a2 i)
 {-# INLINE defaultUnsafeLiftArray2 #-}
 
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
@@ -4,7 +4,7 @@
 {-# OPTIONS_GHC -fno-warn-duplicate-exports #-}
 -- |
 -- Module      : Data.Massiv.Vector
--- Copyright   : (c) Alexey Kuleshevich 2020-2021
+-- Copyright   : (c) Alexey Kuleshevich 2020-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Vector/Stream.hs b/src/Data/Massiv/Vector/Stream.hs
--- a/src/Data/Massiv/Vector/Stream.hs
+++ b/src/Data/Massiv/Vector/Stream.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE CPP #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE LambdaCase #-}
@@ -11,7 +12,7 @@
 {-# OPTIONS_HADDOCK hide, not-home #-}
 -- |
 -- Module      : Data.Massiv.Vector.Stream
--- Copyright   : (c) Alexey Kuleshevich 2019-2021
+-- Copyright   : (c) Alexey Kuleshevich 2019-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -162,9 +163,10 @@
 instance Monad m => Monoid (Steps m e) where
   mempty = empty
   {-# INLINE mempty #-}
+#if !MIN_VERSION_base(4,11,0)
   mappend = append
   {-# INLINE mappend #-}
-
+#endif
 
 instance GHC.Exts.IsList (Steps Id e) where
   type Item (Steps Id e) = e
@@ -202,33 +204,50 @@
   {-# INLINE minimum #-}
 
 
--- TODO: benchmark: `fmap snd . isteps`
 steps :: forall r ix e m . (Monad m, Index ix, Source r e) => Array r ix e -> Steps m e
-steps arr = k `seq` arr `seq` Steps (S.Stream step 0) (LengthExact (coerce k))
-  where
-    k = totalElem $ size arr
-    step i
-      | i < k =
-        let e = unsafeLinearIndex arr i
-         in e `seq` return $ S.Yield e (i + 1)
-      | otherwise = return S.Done
-    {-# INLINE step #-}
-{-# INLINE steps #-}
+steps !arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix -> gix <$> ixRangeSteps (size arr)
+    PrefIndexLinear gi ->
+      Steps (S.Stream step 0) (LengthExact (coerce k))
+      where
+        k = totalElem $ size arr
+        step i
+          | i < k =
+            let e = gi i
+             in e `seq` pure $ S.Yield e (i + 1)
+          | otherwise = pure S.Done
+        {-# INLINE [0] step #-}
+{-# INLINE [1] steps #-}
 
+ixRangeSteps :: forall ix m . (Monad m, Index ix) => Sz ix -> Steps m ix
+ixRangeSteps sz = Steps (S.Stream step initStep) (LengthExact k)
+  where
+    !k = toLinearSz sz
+    !initStep = if k == zeroSz then Nothing else Just zeroIndex
+    step (Just ix) = stepNextMF ix (unSz sz) oneIndex (<) $ \ mIx -> pure $ S.Yield ix mIx
+    step Nothing = pure S.Done
+    {-# INLINE [0] step #-}
+{-# INLINE [1] ixRangeSteps #-}
 
 isteps :: forall r ix e m . (Monad m, Index ix, Source r e) => Array r ix e -> Steps m (ix, e)
-isteps arr = k `seq` arr `seq` Steps (S.Stream step 0) (LengthExact (coerce k))
+isteps !arr =
+  case unsafePrefIndex arr of
+    PrefIndex gix -> (\ !ix -> let e = gix ix in e `seq` (ix, e)) <$> ixRangeSteps sz
+    PrefIndexLinear gi ->
+      let k = totalElem sz
+          step i
+            | i < k =
+              let e = gi i
+               in e `seq` pure $ S.Yield (fromLinearIndex sz i, e) (i + 1)
+            | otherwise = pure S.Done
+          {-# INLINE [0] step #-}
+      in Steps (S.Stream step 0) (LengthExact (coerce k))
   where
-    sz = size arr
-    k = totalElem sz
-    step i
-      | i < k =
-        let e = unsafeLinearIndex arr i
-         in e `seq` return $ S.Yield (fromLinearIndex sz i, e) (i + 1)
-      | otherwise = return S.Done
-    {-# INLINE step #-}
+    !sz = size arr
 {-# INLINE isteps #-}
 
+
 toBundle :: (Monad m, Index ix, Source r e) => Array r ix e -> B.Bundle m v e
 toBundle arr =
   let Steps str k = steps arr
@@ -309,7 +328,7 @@
           unsafeLinearWrite marr i e'
           stepLoad t' (i + 1)
         S.Skip t' -> stepLoad t' i
-        S.Done -> return i
+        S.Done -> pure i
     {-# INLINE stepLoad #-}
 {-# INLINE unstreamMaxM #-}
 
@@ -695,12 +714,12 @@
       case r of
         S.Yield x s' -> do
           b <- f x
-          return $
+          pure $
             case b of
               Nothing -> S.Skip s'
               Just b' -> S.Yield b' s'
-        S.Skip s' -> return $ S.Skip s'
-        S.Done -> return S.Done
+        S.Skip s' -> pure $ S.Skip s'
+        S.Done -> pure S.Done
     {-# INLINE [0] step' #-}
 {-# INLINE mapMaybeStreamM #-}
 
diff --git a/src/Data/Massiv/Vector/Unsafe.hs b/src/Data/Massiv/Vector/Unsafe.hs
--- a/src/Data/Massiv/Vector/Unsafe.hs
+++ b/src/Data/Massiv/Vector/Unsafe.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Vector.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2020-2021
+-- Copyright   : (c) Alexey Kuleshevich 2020-2022
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -64,7 +64,7 @@
 --
 -- @since 0.5.0
 unsafeLast :: Source r e => Vector r e -> e
-unsafeLast v = unsafeLinearIndex v (max 0 (unSz (size v) - 1))
+unsafeLast v = unsafeLinearIndex v (unSz (size v) - 1)
 {-# INLINE unsafeLast #-}
 
 ----------------------
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -12,6 +12,7 @@
 
 -- TODO: fix doctest support
 main :: IO ()
-main = putStrLn "\nDoctests are not supported for ghc version 8.2 and prior as well as 8.10\n"
+main =
+  putStrLn "\nDoctests are not supported for ghc version 8.2 and prior as well as 8.10 and newer\n"
 
 #endif