diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.2.5.0
+version:             0.2.6.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
@@ -89,7 +89,7 @@
                     , Data.Massiv.Array.Delayed.WindowedSpec
                     , Data.Massiv.Array.ManifestSpec
                     , Data.Massiv.Array.Manifest.VectorSpec
-                    --, Data.Massiv.Array.MutableSpec
+                    , Data.Massiv.Array.MutableSpec
                     , Data.Massiv.Array.Ops.ConstructSpec
                     , Data.Massiv.Array.Ops.FoldSpec
                     , Data.Massiv.Array.Ops.SliceSpec
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
@@ -152,6 +152,7 @@
                                      , sum
                                      , zip
                                      )
+
 {- $folding
 
 All folding is done in a row-major order.
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
@@ -46,7 +46,7 @@
 
 import           Control.Exception                   (try)
 import           Control.Monad                       (unless)
-import           Control.Monad.ST                    (RealWorld, runST)
+import           Control.Monad.ST                    (runST)
 import           Data.Foldable                       (Foldable (..))
 import           Data.Massiv.Array.Delayed.Internal
 import           Data.Massiv.Array.Ops.Fold.Internal as M
@@ -59,7 +59,6 @@
 import qualified Data.Vector                         as V
 import           GHC.Base                            hiding (ord)
 import           System.IO.Unsafe                    (unsafePerformIO)
-
 
 #if MIN_VERSION_primitive(0,6,2)
 import           Data.Primitive.Array                (sizeofArray,
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
@@ -82,9 +82,10 @@
 
 -- | Same as `fromLists`, but will throw an error on irregular shaped lists.
 --
--- __Note__: This function is almost the same (modulo customizable computation strategy) if you
--- would turn on @{-# LANGUAGE OverloadedLists #-}@. For that reason you can also use
--- `GHC.Exts.fromList`.
+-- __Note__: This function is the same as if you would turn on @{-\# LANGUAGE OverloadedLists #-}@
+-- extension. For that reason you can also use `GHC.Exts.fromList`.
+--
+-- prop> fromLists' Seq xs == fromList xs
 --
 -- ===__Examples__
 --
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
@@ -13,15 +13,8 @@
 -- Portability : non-portable
 --
 module Data.Massiv.Array.Mutable
-  ( Mutable
-  , MArray
-  , msize
-  , new
-  , thaw
-  , freeze
-  , withMArray
-  , withMArrayST
-  , read
+  ( -- * Element-wise mutation
+    read
   , read'
   , write
   , write'
@@ -29,19 +22,31 @@
   , modify'
   , swap
   , swap'
-  -- * Computation
+  -- ** Operate over `MArray`
+  , Mutable
+  , MArray
+  , msize
+  -- *** Convert
+  , new
+  , thaw
+  , freeze
+  -- *** Create
+  , createArray_
+  , createArray
+  , createArrayST_
+  , createArrayST
+  -- *** Generate
+  , generateArray
+  , generateArrayIO
+  -- *** Unfold
+  , unfoldlPrim_
+  , unfoldlPrim
+  -- *** Modify
+  , withMArray
+  , withMArrayST
+  -- ** Computation
   , RealWorld
   , computeInto
-  -- * Generate (experimental)
-
-  -- $generate
-  {-, generateM
-  , generateLinearM
-  , mapM
-  , imapM
-  , forM
-  , iforM
-  , sequenceM -}
   ) where
 
 import           Prelude                             hiding (mapM, read)
@@ -52,6 +57,7 @@
 import           Data.Massiv.Array.Manifest.Internal
 import           Data.Massiv.Array.Unsafe
 import           Data.Massiv.Core.Common
+import           Data.Massiv.Core.Scheduler
 
 -- | Initialize a new mutable array. Negative size will result in an empty array.
 new :: (Mutable r ix e, PrimMonad m) => ix -> m (MArray (PrimState m) r ix e)
@@ -69,6 +75,184 @@
 {-# INLINE freeze #-}
 
 
+-- | Create a new array by supplying an action that will fill the new blank mutable array. Use
+-- `createArray` if you'd like to keep the result of the filling function.
+--
+-- ====__Examples__
+--
+-- >>> createArray_ Seq (Ix1 2) (\ marr -> write marr 0 10 >> write marr 1 11) :: IO (Array P Ix1 Int)
+-- (Array P Seq (2)
+--   [ 10,11 ])
+--
+-- @since 0.2.6
+--
+createArray_ ::
+     (Mutable r ix e, PrimMonad m)
+  => Comp -- ^ Computation strategy to use after `MArray` gets frozen and onward.
+  -> ix -- ^ Size of the newly created array
+  -> (MArray (PrimState m) r ix e -> m a)
+  -- ^ An action that should fill all elements of the brand new mutable array
+  -> m (Array r ix e)
+createArray_ comp sz action = fmap snd $ createArray comp sz action
+{-# INLINE createArray_ #-}
+
+-- | Just like `createArray_`, but together with `Array` it returns the result of the filling action.
+--
+-- @since 0.2.6
+--
+createArray ::
+     (Mutable r ix e, PrimMonad m)
+  => Comp -- ^ Computation strategy to use after `MArray` gets frozen and onward.
+  -> ix -- ^ Size of the newly created array
+  -> (MArray (PrimState m) r ix e -> m a)
+  -- ^ An action that should fill all elements of the brand new mutable array
+  -> m (a, Array r ix e)
+createArray comp sz action = do
+  marr <- new sz
+  a <- action marr
+  arr <- unsafeFreeze comp marr
+  return (a, arr)
+{-# INLINE createArray #-}
+
+-- | Just like `createArray_`, but restricted to `ST`.
+--
+-- @since 0.2.6
+--
+createArrayST_ ::
+     Mutable r ix e => Comp -> ix -> (forall s. MArray s r ix e -> ST s a) -> Array r ix e
+createArrayST_ comp sz action = runST $ createArray_ comp sz action
+{-# INLINE createArrayST_ #-}
+
+
+-- | Just like `createArray`, but restricted to `ST`.
+--
+-- @since 0.2.6
+--
+createArrayST ::
+     Mutable r ix e => Comp -> ix -> (forall s. MArray s r ix e -> ST s a) -> (a, Array r ix e)
+createArrayST comp sz action = runST $ createArray comp sz action
+{-# INLINE createArrayST #-}
+
+
+-- | Sequentially generate a pure array. Much like `makeArray` creates a pure array this function
+-- will use `Mutable` interface to generate a pure `Array` in the end, except that computation
+-- strategy is ignored. Element producing function no longer has to be pure but is a stateful
+-- action, since it is restricted to `PrimMonad` and allows for sharing the state between
+-- computation of each element, which could be arbitrary effects if that monad is `IO`.
+--
+-- @since 0.2.6
+--
+-- ====__Examples__
+--
+-- >>> import Data.IORef
+-- >>> ref <- newIORef (0 :: Int)
+-- >>> generateArray Seq (Ix1 6) (\ i -> modifyIORef' ref (+i) >> print i >> pure i) :: IO (Array U Ix1 Int)
+-- 0
+-- 1
+-- 2
+-- 3
+-- 4
+-- 5
+-- (Array U Seq (6)
+--   [ 0,1,2,3,4,5 ])
+-- >>> readIORef ref
+-- 15
+--
+generateArray ::
+     (Mutable r ix e, PrimMonad m)
+  => Comp -- ^ Computation strategy (ingored during generation)
+  -> ix -- ^ Resulting size of the array
+  -> (ix -> m e) -- ^ Element producing generator
+  -> m (Array r ix e)
+generateArray comp sz' gen = do
+  let sz = liftIndex (max 0) sz'
+  marr <- unsafeNew sz
+  iterM_ zeroIndex (msize marr) (pureIndex 1) (<) $ \ix -> gen ix >>= write marr ix
+  unsafeFreeze comp marr
+{-# INLINE generateArray #-}
+
+
+-- | Just like `generateArray`, except this generator __will__ respect the supplied computation
+-- strategy, and for that reason it is restricted to `IO`.
+--
+-- @since 0.2.6
+generateArrayIO ::
+     (Mutable r ix e)
+  => Comp
+  -> ix
+  -> (ix -> IO e)
+  -> IO (Array r ix e)
+generateArrayIO comp sz' gen = do
+  case comp of
+    Seq -> generateArray comp sz' gen
+    ParOn wids -> do
+      let sz = liftIndex (max 0) sz'
+      marr <- unsafeNew sz
+      withScheduler_ wids $ \scheduler ->
+        splitLinearlyWithM_
+          (numWorkers scheduler)
+          (scheduleWork scheduler)
+          (totalElem sz)
+          (gen . fromLinearIndex sz)
+          (unsafeLinearWrite marr)
+      unsafeFreeze comp marr
+{-# INLINE generateArrayIO #-}
+
+-- | Sequentially unfold an array from the left.
+--
+-- @since 0.2.6
+--
+-- ====__Examples__
+--
+-- Create an array with Fibonacci numbers while performing and `IO` action on the accumulator for
+-- each element of the array.
+--
+-- >>> unfoldlPrim_ Seq  (Ix1 10) (\a@(f0, f1) _ -> let fn = f0 + f1 in print a >> return ((f1, fn), f0)) (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)
+-- (Array P Seq (10)
+--   [ 0,1,1,2,3,5,8,13,21,34 ])
+--
+unfoldlPrim_ ::
+     (Mutable r ix e, PrimMonad m)
+  => Comp -- ^ Computation strategy (ignored during initial creation)
+  -> ix -- ^ Size of the desired array
+  -> (a -> ix -> m (a, e)) -- ^ Unfolding action
+  -> a -- ^ Initial accumulator
+  -> m (Array r ix e)
+unfoldlPrim_ comp sz gen acc0 = fmap snd $ unfoldlPrim comp sz gen acc0
+{-# INLINE unfoldlPrim_ #-}
+
+
+-- | Just like `unfoldlPrim`, but also returns the final value of the accumulator.
+--
+-- @since 0.2.6
+--
+unfoldlPrim ::
+     (Mutable r ix e, PrimMonad m)
+  => Comp -- ^ Computation strategy (ignored during initial creation)
+  -> ix -- ^ Size of the desired array
+  -> (a -> ix -> m (a, e)) -- ^ Unfolding action
+  -> a -- ^ Initial accumulator
+  -> m (a, Array r ix e)
+unfoldlPrim comp sz gen acc0 =
+  createArray comp sz $ \marr ->
+    let sz' = msize marr
+     in iterLinearM sz' 0 (totalElem sz') 1 (<) acc0 $ \i ix acc -> do
+          (acc', e) <- gen acc ix
+          unsafeLinearWrite marr i e
+          return acc'
+{-# INLINE unfoldlPrim #-}
+
+
 -- | Create a copy of a pure array, mutate it in place and return its frozen version.
 --
 -- @since 0.2.2
@@ -186,140 +370,3 @@
       else ix1
 {-# INLINE swap' #-}
 
-{- Disabled until better times. See https://github.com/lehins/massiv/issues/24
-
-unsafeLinearFillM :: (Mutable r ix e, Monad m) =>
-                     MArray RealWorld r ix e -> (Int -> m e) -> WorldState -> m WorldState
-unsafeLinearFillM ma f (State s_#) = go 0# s_#
-  where
-    !(I# k#) = totalElem (msize ma)
-    go i# s# =
-      case i# <# k# of
-        0# -> return (State s#)
-        _ -> do
-          let i = I# i#
-          res <- f i
-          State s'# <- unsafeLinearWriteA ma i res (State s#)
-          go (i# +# 1#) s'#
-{-# INLINE unsafeLinearFillM #-}
-
-
--- | /O(n)/ - Same as `generateM` but using a flat index.
---
--- @since 0.1.1
-generateLinearM :: (Monad m, Mutable r ix e) => Comp -> ix -> (Int -> m e) -> m (Array r ix e)
-generateLinearM comp sz f = do
-  (s, mba) <- unsafeNewA (liftIndex (max 0) sz) (State (noDuplicate# realWorld#))
-  s' <- unsafeLinearFillM mba f s
-  (_, ba) <- unsafeFreezeA comp mba s'
-  return ba
-{-# INLINE generateLinearM #-}
-
--- | /O(n)/ - Generate an array monadically using it's mutable interface. Computation will be done
-  -- sequentially, regardless of `Comp` argument.
---
--- @since 0.1.1
-generateM :: (Monad m, Mutable r ix e) => Comp -> ix -> (ix -> m e) -> m (Array r ix e)
-generateM comp sz f = generateLinearM comp sz (f . fromLinearIndex sz)
-{-# INLINE generateM #-}
-
-
--- | /O(n)/ - Map an index aware monadic action over an Array. This operation will force computation
--- sequentially and will result in a manifest Array.
---
--- @since 0.1.1
-imapM
-  :: (Monad m, Source r ix e, Mutable r' ix e') =>
-     r' -> (ix -> e -> m e') -> Array r ix e -> m (Array r' ix e')
-imapM _ f arr =
-  generateLinearM (getComp arr) sz (\ !i -> f (fromLinearIndex sz i) (unsafeLinearIndex arr i))
-  where
-    !sz = size arr
-{-# INLINE imapM #-}
-
--- | /O(n)/ - Map a monadic action over an Array. This operation will force computation sequentially
--- and will result in a manifest Array.
---
--- @since 0.1.1
---
--- ====__Examples__
---
--- >>> mapM P (\i -> Just (i*i)) $ range Seq 0 5
--- Just (Array P Seq (5)
---   [ 0,1,4,9,16 ])
---
-mapM
-  :: (Monad m, Source r ix e, Mutable r' ix e') =>
-     r' -> (e -> m e') -> Array r ix e -> m (Array r' ix e')
-mapM r f = imapM r (const f)
-{-# INLINE mapM #-}
-
-
--- | /O(n)/ - Same as `mapM`, but with its arguments flipped.
---
--- @since 0.1.1
-forM ::
-     (Monad m, Source r ix e, Mutable r' ix e')
-  => r'
-  -> Array r ix e
-  -> (e -> m e')
-  -> m (Array r' ix e')
-forM r = flip (mapM r)
-{-# INLINE forM #-}
-
-
--- | /O(n)/ - Same as `imapM`, but with its arguments flipped.
---
--- @since 0.1.1
-iforM :: (Monad m, Source r ix e, Mutable r' ix e') =>
-         r' -> Array r ix e -> (ix -> e -> m e') -> m (Array r' ix e')
-iforM r = flip (imapM r)
-{-# INLINE iforM #-}
-
-
--- | /O(n)/ - Sequence monadic actions in a source Array. This operation will force the computation
--- sequentially and will result in a manifest Array.
---
--- @since 0.1.1
-sequenceM
-  :: (Monad m, Source r ix (m e), Mutable r' ix e) =>
-     r' -> Array r ix (m e) -> m (Array r' ix e)
-sequenceM r = mapM r id
-{-# INLINE sequenceM #-}
--}
-
-{- $generate
-
-Functions in this section has been removed until better times due to a known bug https://github.com/lehins/massiv/issues/24
-
--}
-
-{- Disabled until better times
-
-Functions in this section can monadically generate manifest arrays using their associated mutable
-interface. Due to the sequential nature of monads generation is done also sequentially regardless of
-supplied computation strategy. All of functions here are very much experimental, so please
-<https://github.com/lehins/massiv/issues/new report an issue> if you see something not working
-properly.
-
-Here is a very imperative like for loop that creates an array while performing a side effect for
-each newly created element:
-
-@
-printSquare :: Int -> IO (Array P Ix1 Int)
-printSquare n = forM P (range Seq 0 n) $ \i -> do
-  let e = i*i
-  putStrLn $ "Element at index: " ++ show i ++ " = " ++ show e ++ ";"
-  return e
-@
-
->>> printSquare 5
-Element at index: 0 = 0;
-Element at index: 1 = 1;
-Element at index: 2 = 4;
-Element at index: 3 = 9;
-Element at index: 4 = 16;
-(Array P Seq (5)
-  [ 0,1,4,9,16 ])
-
--}
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
@@ -1,7 +1,11 @@
 {-# LANGUAGE BangPatterns          #-}
+{-# LANGUAGE DataKinds             #-}
+{-# LANGUAGE ExplicitForAll        #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
 {-# LANGUAGE TypeFamilies          #-}
+{-# LANGUAGE TypeOperators         #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Construct
 -- Copyright   : (c) Alexey Kuleshevich 2018
@@ -11,18 +15,29 @@
 -- Portability : non-portable
 --
 module Data.Massiv.Array.Ops.Construct
-  ( makeArray
+  ( -- ** From a function
+    makeArray
   , makeArrayR
   , makeVectorR
   , singleton
+    -- *** Applicative
+  , makeArrayA
+  , makeArrayAR
+    -- ** Enumeration
   , range
   , rangeStep
   , enumFromN
   , enumFromStepN
+    -- ** Expansion
+  , expandWithin
+  , expandWithin'
+  , expandOuter
+  , expandInner
   ) where
 
 import           Data.Massiv.Array.Delayed.Internal
 import           Data.Massiv.Core.Common
+import           Data.Massiv.Array.Ops.Map          as A
 import           Prelude                            as P
 
 
@@ -51,7 +66,23 @@
 makeVectorR _ = makeArray
 {-# INLINE makeVectorR #-}
 
+-- | Similar to `makeArray`, but construct the array sequentially using an `Applicative` interface
+-- disregarding the supplied `Comp`.
+--
+-- @since 0.2.6
+--
+makeArrayA :: (Mutable r a b, Applicative f) => Comp -> a -> (a -> f b) -> f (Array r a b)
+makeArrayA comp sz f = traverseA f $ makeArrayR D comp sz id
+{-# INLINE makeArrayA #-}
 
+-- | Same as `makeArrayA`, but with ability to supply result array representation.
+--
+-- @since 0.2.6
+--
+makeArrayAR :: (Mutable r a b, Applicative f) => r -> Comp -> a -> (a -> f b) -> f (Array r a b)
+makeArrayAR _ = makeArrayA
+{-# INLINE makeArrayAR #-}
+
 -- | Create a vector with a range of @Int@s incremented by 1.
 -- @range k0 k1 == rangeStep k0 k1 1@
 --
@@ -120,3 +151,115 @@
 {-# INLINE enumFromStepN #-}
 
 
+-- | Function that expands an array to one with a higher dimension.
+--
+-- This is useful for constructing arrays where there is shared computation
+-- between multiple cells.  The makeArray method of constructing arrays:
+--
+-- > makeArray :: Construct r ix e => Comp -> ix -> (ix -> e) -> Array r ix e
+--
+-- ...runs a function @ix -> e@ at every array index. This is inefficient if
+-- there is a substantial amount of repeated computation that could be shared
+-- while constructing elements on the same dimension. The expand functions make
+-- this possible. First you construct an @Array r (Lower ix) a@ of one fewer
+-- dimensions where @a@ is something like @`Array` r `Ix1` a@ or @`Array` r `Ix2` a@. Then
+-- you use 'expandWithin' and a creation function @a -> Int -> b@ to create an
+-- @`Array` `D` `Ix2` b@ or @`Array` `D` `Ix3` b@ respectfully.
+--
+-- @since 0.2.6
+--
+-- ====__Examples__
+--
+-- >>> a = makeArrayR U Seq (Ix1 6) (+10) -- Imagine (+10) is some expensive function
+-- >>> a
+-- (Array U Seq (6)
+--   [ 10,11,12,13,14,15 ])
+-- >>> expandWithin Dim1 5 (\ e j -> (j + 1) * 100 + e) a :: Array D Ix2 Int
+-- (Array D Seq (6 :. 5)
+--   [ [ 110,210,310,410,510 ]
+--   , [ 111,211,311,411,511 ]
+--   , [ 112,212,312,412,512 ]
+--   , [ 113,213,313,413,513 ]
+--   , [ 114,214,314,414,514 ]
+--   , [ 115,215,315,415,515 ]
+--   ])
+-- >>> expandWithin Dim2 5 (\ e j -> (j + 1) * 100 + e) a :: Array D Ix2 Int
+-- (Array D Seq (5 :. 6)
+--   [ [ 110,111,112,113,114,115 ]
+--   , [ 210,211,212,213,214,215 ]
+--   , [ 310,311,312,313,314,315 ]
+--   , [ 410,411,412,413,414,415 ]
+--   , [ 510,511,512,513,514,515 ]
+--   ])
+--
+expandWithin
+  :: (IsIndexDimension ix n, Manifest r (Lower ix) a)
+  => Dimension n
+  -> Int
+  -> (a -> Int -> b)
+  -> Array r (Lower ix) a
+  -> Array D ix b
+expandWithin dim k f arr = do
+  makeArray (getComp arr) sz $ \ix ->
+    let (i, ixl) = pullOutDimension ix dim
+     in f (unsafeIndex arr ixl) i
+  where
+    szl = size arr
+    sz = insertDimension szl dim k
+{-# INLINE expandWithin #-}
+
+-- | Similar to `expandWithin`, except that dimension is specified at a value level, which means it
+-- will throw an exception on an invalid dimension.
+--
+-- @since 0.2.6
+expandWithin'
+  :: (Index ix, Manifest r (Lower ix) a)
+  => Dim
+  -> Int
+  -> (a -> Int -> b)
+  -> Array r (Lower ix) a
+  -> Array D ix b
+expandWithin' dim k f arr =
+  makeArray (getComp arr) sz $ \ix ->
+    let (i, ixl) = pullOutDim' ix dim
+     in f (unsafeIndex arr ixl) i
+  where
+    szl = size arr
+    sz = insertDim' szl dim k
+{-# INLINE expandWithin' #-}
+
+-- | Similar to `expandWithin`, except it uses the outermost dimension.
+--
+-- @since 0.2.6
+expandOuter
+  :: (Index ix, Manifest r (Lower ix) a)
+  => Int
+  -> (a -> Int -> b)
+  -> Array r (Lower ix) a
+  -> Array D ix b
+expandOuter k f arr =
+  makeArray (getComp arr) sz $ \ix ->
+    let (i, ixl) = unconsDim ix
+     in f (unsafeIndex arr ixl) i
+  where
+    szl = size arr
+    sz = consDim k szl
+{-# INLINE expandOuter #-}
+
+-- | Similar to `expandWithin`, except it uses the innermost dimension.
+--
+-- @since 0.2.6
+expandInner
+  :: (Index ix, Manifest r (Lower ix) a)
+  => Int
+  -> (a -> Int -> b)
+  -> Array r (Lower ix) a
+  -> Array D ix b
+expandInner k f arr =
+  makeArray (getComp arr) sz $ \ix ->
+    let (ixl, i) = unsnocDim ix
+     in f (unsafeIndex arr ixl) i
+  where
+    szl = size arr
+    sz = snocDim szl k
+{-# INLINE expandInner #-}
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
@@ -347,13 +347,13 @@
 {-# INLINE or #-}
 
 
--- | Determines whether all element of the array satisfy the predicate.
+-- | /O(n)/ - Determines whether all element of the array satisfy the predicate.
 all :: Source r ix e =>
        (e -> Bool) -> Array r ix e -> Bool
 all f = foldlInternal (\acc el -> acc && f el) True (&&) True
 {-# INLINE all #-}
 
--- | Determines whether any element of the array satisfies the predicate.
+-- | /O(n)/ - Determines whether any element of the array satisfies the predicate.
 any :: Source r ix e =>
        (e -> Bool) -> Array r ix e -> Bool
 any f = foldlInternal (\acc el -> acc || f el) False (||) False
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
@@ -341,7 +341,7 @@
 {-# INLINE ifoldrP #-}
 
 
--- | This folding function breaks referencial transparency on some functions
+-- | This folding function breaks referential transparency on some functions
 -- @f@, therefore it is kept here for internal use only.
 foldlInternal :: Source r ix e =>
          (a -> e -> a) -> a -> (b -> a -> b) -> b -> Array r ix e -> b
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
@@ -1,6 +1,7 @@
 {-# LANGUAGE BangPatterns          #-}
 {-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables   #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Map
 -- Copyright   : (c) Alexey Kuleshevich 2018
@@ -12,11 +13,34 @@
 module Data.Massiv.Array.Ops.Map
   ( map
   , imap
+  -- ** Traversing
+  , traverseA
+  , itraverseA
+  , traverseAR
+  , itraverseAR
   -- ** Monadic
+  -- *** Sequential
+  , mapM
+  , mapMR
+  , forM
+  , forMR
+  , imapM
+  , imapMR
+  , iforM
+  , iforMR
   , mapM_
   , forM_
   , imapM_
   , iforM_
+  -- *** Parallelizable
+  , mapIO
+  , mapIO_
+  , imapIO
+  , imapIO_
+  , forIO
+  , forIO_
+  , iforIO
+  , iforIO_
   , mapP_
   , imapP_
   -- ** Zipping
@@ -32,15 +56,26 @@
   ) where
 
 
-import           Control.Monad                      (void, when)
+import           Control.Monad                       (void, when)
+import           Control.Monad.ST                    (runST)
+import           Data.Foldable                       (foldlM)
 import           Data.Massiv.Array.Delayed.Internal
+import           Data.Massiv.Array.Mutable
+import           Data.Massiv.Array.Ops.Fold.Internal (foldrFB)
 import           Data.Massiv.Core.Common
 import           Data.Massiv.Core.Scheduler
-import           Data.Monoid                        ((<>))
-import           Prelude                            hiding (map, mapM, mapM_,
-                                                     unzip, unzip3, zip, zip3,
-                                                     zipWith, zipWith3)
+import           Data.Monoid                         ((<>))
+import           GHC.Base                            (build)
+import           Prelude                             hiding (map, mapM, mapM_,
+                                                      traverse, unzip, unzip3,
+                                                      zip, zip3, zipWith,
+                                                      zipWith3)
+import qualified Prelude                             as Prelude (traverse)
 
+--------------------------------------------------------------------------------
+-- map -------------------------------------------------------------------------
+--------------------------------------------------------------------------------
+
 -- | Map a function over an array
 map :: Source r ix e' => (e' -> e) -> Array r ix e' -> Array D ix e
 map f = imap (const f)
@@ -51,6 +86,10 @@
 imap f !arr = DArray (getComp arr) (size arr) (\ !ix -> f ix (unsafeIndex arr ix))
 {-# INLINE imap #-}
 
+--------------------------------------------------------------------------------
+-- zip -------------------------------------------------------------------------
+--------------------------------------------------------------------------------
+
 -- | Zip two arrays
 zip :: (Source r1 ix e1, Source r2 ix e2)
     => Array r1 ix e1 -> Array r2 ix e2 -> Array D ix (e1, e2)
@@ -74,7 +113,9 @@
 unzip3 arr = (map (\ (e, _, _) -> e) arr, map (\ (_, e, _) -> e) arr, map (\ (_, _, e) -> e) arr)
 {-# INLINE unzip3 #-}
 
-
+--------------------------------------------------------------------------------
+-- zipWith ---------------------------------------------------------------------
+--------------------------------------------------------------------------------
 
 -- | Zip two arrays with a function. Resulting array will be an intersection of
 -- source arrays in case their dimensions do not match.
@@ -115,8 +156,195 @@
     f ix (unsafeIndex arr1 ix) (unsafeIndex arr2 ix) (unsafeIndex arr3 ix)
 {-# INLINE izipWith3 #-}
 
+--------------------------------------------------------------------------------
+-- traverse --------------------------------------------------------------------
+--------------------------------------------------------------------------------
 
+-- | Traverse with an `Applicative` action over an array sequentially.
+--
+-- @since 0.2.6
+--
+traverseA ::
+     (Source r' ix a, Mutable r ix b, Applicative f)
+  => (a -> f b)
+  -> Array r' ix a
+  -> f (Array r ix b)
+traverseA f arr = loadList <$> Prelude.traverse f (build (\c n -> foldrFB c n arr))
+  where
+    loadList xs =
+      runST $ do
+        marr <- unsafeNew (size arr)
+        _ <- foldlM (\i e -> unsafeLinearWrite marr i e >> return (i + 1)) 0 xs
+        unsafeFreeze (getComp arr) marr
+    {-# INLINE loadList #-}
+{-# INLINE traverseA #-}
 
+-- | Traverse with an `Applicative` index aware action over an array sequentially.
+--
+-- @since 0.2.6
+--
+itraverseA ::
+     (Source r' ix a, Mutable r ix b, Applicative f)
+  => (ix -> a -> f b)
+  -> Array r' ix a
+  -> f (Array r ix b)
+itraverseA f arr =
+  fmap loadList $ Prelude.traverse (uncurry f) $ build (\c n -> foldrFB c n (zipWithIndex arr))
+  where
+    loadList xs =
+      runST $ do
+        marr <- unsafeNew (size arr)
+        _ <- foldlM (\i e -> unsafeLinearWrite marr i e >> return (i + 1)) 0 xs
+        unsafeFreeze (getComp arr) marr
+    {-# INLINE loadList #-}
+{-# INLINE itraverseA #-}
+
+
+
+-- | Same as `traverseA`, except with ability to specify representation.
+--
+-- @since 0.2.6
+--
+traverseAR ::
+     (Source r' ix a, Mutable r ix b, Applicative f)
+  => r
+  -> (a -> f b)
+  -> Array r' ix a
+  -> f (Array r ix b)
+traverseAR _ = traverseA
+{-# INLINE traverseAR #-}
+
+-- | Same as `itraverseA`, except with ability to specify representation.
+--
+-- @since 0.2.6
+--
+itraverseAR ::
+     (Source r' ix a, Mutable r ix b, Applicative f)
+  => r
+  -> (ix -> a -> f b)
+  -> Array r' ix a
+  -> f (Array r ix b)
+itraverseAR _ = itraverseA
+{-# INLINE itraverseAR #-}
+
+zipWithIndex :: forall r ix e . Source r ix e => Array r ix e -> Array D ix (ix, e)
+zipWithIndex arr = zip (makeArray mempty (size arr) id :: Array D ix ix) arr
+{-# INLINE zipWithIndex #-}
+
+--------------------------------------------------------------------------------
+-- mapM ------------------------------------------------------------------------
+--------------------------------------------------------------------------------
+
+-- | Map a monadic action over an array sequentially.
+--
+-- @since 0.2.6
+--
+mapM ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => (a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+mapM = traverseA
+{-# INLINE mapM #-}
+
+
+-- | Same as `mapM`, except with ability to specify result representation.
+--
+-- @since 0.2.6
+--
+mapMR ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => r
+  -> (a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+mapMR _ = traverseA
+{-# INLINE mapMR #-}
+
+
+-- | Same as `mapM` except with arguments flipped.
+--
+-- @since 0.2.6
+--
+forM ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => Array r' ix a
+  -> (a -> m b)
+  -> m (Array r ix b)
+forM = flip traverseA
+{-# INLINE forM #-}
+
+
+-- | Same as `forM`, except with ability to specify result representation.
+--
+-- @since 0.2.6
+--
+forMR ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => r
+  -> Array r' ix a
+  -> (a -> m b)
+  -> m (Array r ix b)
+forMR _ = flip traverseA
+{-# INLINE forMR #-}
+
+
+
+-- | Map a monadic action over an array sequentially.
+--
+-- @since 0.2.6
+--
+imapM ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => (ix -> a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+imapM = itraverseA
+{-# INLINE imapM #-}
+
+
+-- | Same as `imapM`, except with ability to specify result representation.
+--
+-- @since 0.2.6
+--
+imapMR ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => r
+  -> (ix -> a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+imapMR _ = itraverseA
+{-# INLINE imapMR #-}
+
+
+
+-- | Same as `forM`, except map an index aware action.
+--
+-- @since 0.2.6
+--
+iforM ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => (ix -> a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+iforM = itraverseA
+{-# INLINE iforM #-}
+
+
+-- | Same as `iforM`, except with ability to specify result representation.
+--
+-- @since 0.2.6
+--
+iforMR ::
+     (Source r' ix a, Mutable r ix b, Monad m)
+  => r
+  -> (ix -> a -> m b)
+  -> Array r' ix a
+  -> m (Array r ix b)
+iforMR _ = itraverseA
+{-# INLINE iforMR #-}
+
+
 -- | Map a monadic function over an array sequentially, while discarding the result.
 --
 -- ==== __Examples__
@@ -157,11 +385,87 @@
 {-# INLINE iforM_ #-}
 
 
+-- | Map an `IO` action over an `Array`. Underlying computation strategy is respected and will be
+-- parallelized when requested. Unfortunately no fusion is possible and new array will be create
+-- upon each call.
+--
+-- @since 0.2.6
+mapIO ::
+     (Source r' ix a, Mutable r ix b) => (a -> IO b) -> Array r' ix a -> IO (Array r ix b)
+mapIO action = imapIO (const action)
+{-# INLINE mapIO #-}
 
+-- | Similar to `mapIO`, but ignores the result of mapping action and does not create a resulting
+-- array, therefore it is faster. Use this instead of `mapIO` when result is irrelevant.
+--
+-- @since 0.2.6
+mapIO_ :: Source r b e => (e -> IO a) -> Array r b e -> IO ()
+mapIO_ action = imapIO_ (const action)
+{-# INLINE mapIO_ #-}
+
+-- | Same as `mapIO_`, but map an index aware action instead.
+--
+-- @since 0.2.6
+imapIO_ :: Source r ix e => (ix -> e -> IO a) -> Array r ix e -> IO ()
+imapIO_ action arr =
+  case getComp arr of
+    Seq -> imapM_ action arr
+    ParOn wids -> do
+      let sz = size arr
+      withScheduler_ wids $ \scheduler ->
+        splitLinearlyWith_
+          (numWorkers scheduler)
+          (scheduleWork scheduler)
+          (totalElem sz)
+          (unsafeLinearIndex arr)
+          (\i -> void . action (fromLinearIndex sz i))
+{-# INLINE imapIO_ #-}
+
+
+-- | Same as `mapIO` but map an index aware action instead.
+--
+-- @since 0.2.6
+imapIO ::
+     (Source r' ix a, Mutable r ix b) => (ix -> a -> IO b) -> Array r' ix a -> IO (Array r ix b)
+imapIO action arr = generateArrayIO (getComp arr) (size arr) $ \ix -> action ix (unsafeIndex arr ix)
+{-# INLINE imapIO #-}
+
+-- | Same as `mapIO` but with arguments flipped.
+--
+-- @since 0.2.6
+forIO ::
+     (Source r' ix a, Mutable r ix b) => Array r' ix a -> (a -> IO b) -> IO (Array r ix b)
+forIO = flip mapIO
+{-# INLINE forIO #-}
+
+-- | Same as `mapIO_` but with arguments flipped.
+--
+-- @since 0.2.6
+forIO_ :: Source r ix e => Array r ix e -> (e -> IO a) -> IO ()
+forIO_ = flip mapIO_
+{-# INLINE forIO_ #-}
+
+-- | Same as `imapIO` but with arguments flipped.
+--
+-- @since 0.2.6
+iforIO ::
+     (Source r' ix a, Mutable r ix b) => Array r' ix a -> (ix -> a -> IO b) -> IO (Array r ix b)
+iforIO = flip imapIO
+{-# INLINE iforIO #-}
+
+-- | Same as `imapIO_` but with arguments flipped.
+--
+-- @since 0.2.6
+iforIO_ :: Source r ix a => Array r ix a -> (ix -> a -> IO b) -> IO ()
+iforIO_ = flip imapIO_
+{-# INLINE iforIO_ #-}
+
+
 -- | Map an IO action, over an array in parallel, while discarding the result.
 mapP_ :: Source r ix a => (a -> IO b) -> Array r ix a -> IO ()
 mapP_ f = imapP_ (const f)
 {-# INLINE mapP_ #-}
+{-# DEPRECATED mapP_ "In favor of 'mapIO_'" #-}
 
 
 -- | Map an index aware IO action, over an array in parallel, while
@@ -183,3 +487,4 @@
       iterLinearM_ sz slackStart totalLength 1 (<) $ \ !i ix -> do
         void $ f ix (unsafeLinearIndex arr i)
 {-# INLINE imapP_ #-}
+{-# DEPRECATED imapP_ "In favor of 'imapIO_'" #-}
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
@@ -178,18 +178,20 @@
 -- | To be deprecated in favor of `setDim'`.
 setIndex' :: Index ix => ix -> Dim -> Int -> ix
 setIndex' ix dim i =
-  case setIndex ix dim i of
+  case setDim ix dim i of
     Just ix' -> ix'
     Nothing  -> errorDim "setIndex'" dim
 {-# INLINE [1] setIndex' #-}
+{-# DEPRECATED setIndex' "In favor of `setDim'`" #-}
 
 -- | To be deprecated in favor of `getDim'`.
 getIndex' :: Index ix => ix -> Dim -> Int
 getIndex' ix dim =
-  case getIndex ix dim of
+  case getDim ix dim of
     Just ix' -> ix'
     Nothing  -> errorDim "getIndex'" dim
 {-# INLINE [1] getIndex' #-}
+{-# DEPRECATED getIndex' "In favor of `getDim'`" #-}
 
 dropDim' :: Index ix => ix -> Dim -> Lower ix
 dropDim' ix dim =
diff --git a/src/Data/Massiv/Core/Index/Class.hs b/src/Data/Massiv/Core/Index/Class.hs
--- a/src/Data/Massiv/Core/Index/Class.hs
+++ b/src/Data/Massiv/Core/Index/Class.hs
@@ -115,13 +115,11 @@
   setDim = setIndex
   {-# INLINE [1] setDim #-}
 
-  -- TODO: depricate
   -- | Extract the value index has at specified dimension. To be deprecated.
   getIndex :: ix -> Dim -> Maybe Int
   getIndex = getDim
   {-# INLINE [1] getIndex #-}
 
-  -- TODO: depricate
   -- | Set the value for an index at specified dimension. To be deprecated.
   setIndex :: ix -> Dim -> Int -> Maybe ix
   setIndex = setDim
@@ -249,6 +247,8 @@
       !(inc, incIxL) = unconsDim incIx
   {-# INLINE iterM_ #-}
 
+{-# DEPRECATED getIndex "In favor of 'getDim'" #-}
+{-# DEPRECATED setIndex "In favor of 'setDim'" #-}
 
 instance Index Ix1T where
   type Dimensions Ix1T = 1
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
@@ -2,6 +2,7 @@
 {-# LANGUAGE PatternSynonyms            #-}
 
 #if __GLASGOW_HASKELL__ >= 800
+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 #else
 {-# LANGUAGE GADTs                      #-}
@@ -56,6 +57,7 @@
 
 #if __GLASGOW_HASKELL__ >= 800
 newtype Stride ix = SafeStride ix deriving (Eq, Ord, NFData)
+{-# COMPLETE Stride #-}
 #else
 -- There is an issue in GHC 7.10 which prevents from placing `Index` constraint on a pattern.
 data Stride ix where
@@ -68,15 +70,16 @@
 #endif
 
 
-instance Index ix => Show (Stride ix) where
-  show (SafeStride ix) = "Stride (" ++ show ix ++ ")"
-
-
 -- | A safe bidirectional pattern synonym for `Stride` construction that will make sure stride
 -- elements are always positive.
 pattern Stride :: Index ix => ix -> Stride ix
 pattern Stride ix <- SafeStride ix where
         Stride ix = SafeStride (liftIndex (max 1) ix)
+
+
+instance Index ix => Show (Stride ix) where
+  show (SafeStride ix) = "Stride (" ++ show ix ++ ")"
+
 
 -- | Just a helper function for unwrapping `Stride`.
 unStride :: Stride ix -> ix
diff --git a/src/Data/Massiv/Core/Iterator.hs b/src/Data/Massiv/Core/Iterator.hs
--- a/src/Data/Massiv/Core/Iterator.hs
+++ b/src/Data/Massiv/Core/Iterator.hs
@@ -14,6 +14,7 @@
   , loopDeepM
   , splitLinearly
   , splitLinearlyWith_
+  , splitLinearlyWithM_
   ) where
 
 
@@ -72,9 +73,19 @@
 
 
 splitLinearlyWith_ :: Monad m => Int -> (m () -> m a) -> Int -> (Int -> b) -> (Int -> b -> m ()) -> m a
-splitLinearlyWith_ numChunks with totalLength index write =
+splitLinearlyWith_ numChunks with totalLength index =
+  splitLinearlyWithM_ numChunks with totalLength (pure . index)
+{-# INLINE splitLinearlyWith_ #-}
+
+
+-- | Interator tha can be used to split computation jobs
+--
+-- @since 0.2.6.0
+splitLinearlyWithM_ ::
+     Monad m => Int -> (m () -> m a) -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m a
+splitLinearlyWithM_ numChunks with totalLength make write =
   splitLinearly numChunks totalLength  $ \chunkLength slackStart -> do
     loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
-      with $ loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> write k (index k)
-    with $ loopM_ slackStart (< totalLength) (+ 1) $ \ !k -> write k (index k)
-{-# INLINE splitLinearlyWith_ #-}
+      with $ loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> make k >>= write k
+    with $ loopM_ slackStart (< totalLength) (+ 1) $ \ !k -> make k >>= write k
+{-# INLINE splitLinearlyWithM_ #-}
diff --git a/tests/Data/Massiv/Array/MutableSpec.hs b/tests/Data/Massiv/Array/MutableSpec.hs
new file mode 100644
--- /dev/null
+++ b/tests/Data/Massiv/Array/MutableSpec.hs
@@ -0,0 +1,83 @@
+{-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE MonoLocalBinds        #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module Data.Massiv.Array.MutableSpec (spec) where
+
+import           Control.Monad.ST
+import           Data.Functor.Identity
+import           Data.Massiv.CoreArbitrary as A
+import           Data.Proxy
+import           Test.Hspec
+import           Test.QuickCheck
+import           Test.QuickCheck.Monadic
+import           Test.QuickCheck.Function
+
+prop_MapMapM :: (Show (Array r ix Int), Eq (Array r ix Int), Mutable r ix Int) =>
+                r -> Proxy ix -> Fun Int Int -> ArrTiny D ix Int -> Property
+prop_MapMapM r _ f (ArrTiny arr) =
+  computeAs r (A.map (apply f) arr) === runIdentity (A.mapMR r (return . apply f) arr)
+
+prop_iMapiMapM :: (Show (Array r ix Int), Eq (Array r ix Int), Mutable r ix Int) =>
+                r -> Proxy ix -> Fun (ix, Int) Int -> ArrTiny D ix Int -> Property
+prop_iMapiMapM r _ f (ArrTiny arr) =
+  computeAs r (A.imap (curry (apply f)) arr) ===
+  runIdentity (A.imapMR r (\ix e -> return $ apply f (ix, e)) arr)
+
+
+prop_generateMakeST :: (Show (Array r ix Int), Eq (Array r ix Int), Mutable r ix Int) =>
+                             r -> Proxy ix -> Arr r ix Int -> Property
+prop_generateMakeST _ _ (Arr arr) =
+  arr === runST (generateArray (getComp arr) (size arr) (return . evaluateAt arr))
+
+prop_generateMakeIO :: (Show (Array r ix Int), Eq (Array r ix Int), Mutable r ix Int) =>
+                             r -> Proxy ix -> Arr r ix Int -> Property
+prop_generateMakeIO _ _ (Arr arr) = monadicIO $ do
+  arr' <- run $ generateArray (getComp arr) (size arr) (return . evaluateAt arr)
+  return (arr === arr')
+
+mutableSpec ::
+     ( Show r
+     , Show (Array r Ix3 Int)
+     , Show (Array r Ix1 Int)
+     , Show (Array r Ix2 Int)
+     , Eq (Array r Ix3 Int)
+     , Eq (Array r Ix1 Int)
+     , Eq (Array r Ix2 Int)
+     , Mutable r Ix3 Int
+     , Mutable r Ix1 Int
+     , Mutable r Ix2 Int
+     )
+  => r
+  -> SpecWith ()
+mutableSpec r = do
+  describe (show r) $ do
+    describe "map == mapM" $ do
+      it "Ix1" $ property $ prop_MapMapM r (Proxy :: Proxy Ix1)
+      it "Ix2" $ property $ prop_MapMapM r (Proxy :: Proxy Ix2)
+      it "Ix3" $ property $ prop_MapMapM r (Proxy :: Proxy Ix3)
+    describe "imap == imapM" $ do
+      it "Ix1" $ property $ prop_iMapiMapM r (Proxy :: Proxy Ix1)
+      it "Ix2T" $ property $ prop_iMapiMapM r (Proxy :: Proxy Ix2)
+      it "Ix3T" $ property $ prop_iMapiMapM r (Proxy :: Proxy Ix3)
+    describe "makeArray == generateArrayST" $ do
+      it "Ix1" $ property $ prop_generateMakeST r (Proxy :: Proxy Ix1)
+      it "Ix2" $ property $ prop_generateMakeST r (Proxy :: Proxy Ix2)
+      it "Ix3" $ property $ prop_generateMakeST r (Proxy :: Proxy Ix3)
+    describe "makeArray == generateArrayIO" $ do
+      it "Ix1" $ property $ prop_generateMakeIO r (Proxy :: Proxy Ix1)
+      it "Ix2" $ property $ prop_generateMakeIO r (Proxy :: Proxy Ix2)
+      it "Ix3" $ property $ prop_generateMakeIO r (Proxy :: Proxy Ix3)
+
+
+generateSpec :: Spec
+generateSpec = do
+  mutableSpec P
+  mutableSpec S
+  mutableSpec U
+  mutableSpec B
+  mutableSpec N
+
+
+spec :: Spec
+spec = describe "GenerateM" generateSpec
diff --git a/tests/Data/Massiv/Array/Ops/ConstructSpec.hs b/tests/Data/Massiv/Array/Ops/ConstructSpec.hs
--- a/tests/Data/Massiv/Array/Ops/ConstructSpec.hs
+++ b/tests/Data/Massiv/Array/Ops/ConstructSpec.hs
@@ -1,9 +1,10 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE GADTs            #-}
-module Data.Massiv.Array.Ops.ConstructSpec (spec) where
+module Data.Massiv.Array.Ops.ConstructSpec where
 
-import           Data.Massiv.CoreArbitrary as A
+import           Data.Massiv.Array                 as A
+import           Data.Massiv.CoreArbitrary         as A
 import           Data.Proxy
 import qualified GHC.Exts                   as GHC (IsList (..))
 import           Prelude                    as P
@@ -91,9 +92,29 @@
   it "toFromListIsList" $ property (prop_toFromListIsList (Proxy :: Proxy Ix3))
   it "excFromToListIx3" $ property prop_excFromToListIx3
 
+mkIntermediate :: Int -> Array U Ix1 Int
+mkIntermediate t = A.fromList Seq [t + 50, t + 75]
 
+initArr :: Array N Ix1 (Array U Ix1 Int)
+initArr = makeArray Seq 3 (\ x -> mkIntermediate x)
+
+initArr2 :: Array N Ix2 (Array U Ix1 Int)
+initArr2 = makeArray Seq (2 :. 2) (\ (x :. y) -> mkIntermediate (x+y))
+
+specExpand :: Spec
+specExpand = do
+  it "expandOuter" $ compute (expandOuter 2 A.index' initArr :: Array D Ix2 Int) `shouldBe`
+    resize' (2 :. 3) (fromList Seq [50, 51, 52, 75, 76, 77] :: Array U Ix1 Int)
+  it "expandInner" $ compute (expandInner 2 A.index' initArr :: Array D Ix2 Int) `shouldBe`
+    resize' (3 :. 2) (fromList Seq [50, 75, 51, 76, 52, 77] :: Array U Ix1 Int)
+  it "expandwithin" $ compute (expandWithin Dim1 2 A.index' initArr2 :: Array D Ix3 Int) `shouldBe`
+    resize' (2 :> 2 :. 2) (fromList Seq [50, 75, 51, 76, 51, 76, 52, 77] :: Array U Ix1 Int)
+  it "expandwithin'" $ compute (expandWithin' 1 2 A.index' initArr2 :: Array D Ix3 Int) `shouldBe`
+    resize' (2 :> 2 :. 2) (fromList Seq [50, 75, 51, 76, 51, 76, 52, 77] :: Array U Ix1 Int)
+
 spec :: Spec
 spec = do
   describe "Ix1" specIx1
   describe "Ix2" specIx2
   describe "Ix3" specIx3
+  describe "Expand" specExpand
diff --git a/tests/Data/Massiv/Array/Ops/TransformSpec.hs b/tests/Data/Massiv/Array/Ops/TransformSpec.hs
--- a/tests/Data/Massiv/Array/Ops/TransformSpec.hs
+++ b/tests/Data/Massiv/Array/Ops/TransformSpec.hs
@@ -4,7 +4,6 @@
 module Data.Massiv.Array.Ops.TransformSpec (spec) where
 
 import           Data.Massiv.CoreArbitrary as A
-import           Data.Maybe                (fromJust)
 import           Data.Typeable             (Typeable)
 import           Test.Hspec
 import           Test.QuickCheck
@@ -15,7 +14,7 @@
   => proxy (r, ix, e) -> DimIx ix -> ArrIx r ix e -> Bool
 prop_ExtractAppend _ (DimIx dim) (ArrIx arr ix) =
   maybe False ((delay arr ==) . uncurry (append' dim)) $
-  A.splitAt dim (fromJust (getIndex ix dim)) arr
+  A.splitAt dim (getDim' ix dim) arr
 
 
 prop_transposeOuterInner :: Arr D Ix2 Int -> Property
diff --git a/tests/Data/Massiv/Array/StencilSpec.hs b/tests/Data/Massiv/Array/StencilSpec.hs
--- a/tests/Data/Massiv/Array/StencilSpec.hs
+++ b/tests/Data/Massiv/Array/StencilSpec.hs
@@ -47,8 +47,8 @@
     ix' =
       liftIndex (* signum s) $
       fromJust $ do
-        i <- getIndex sz r
-        setIndex zeroIndex r i
+        i <- getDim sz r
+        setDim zeroIndex r i
 
 
 stencilSpec :: Spec
diff --git a/tests/Data/Massiv/Core/IndexSpec.hs b/tests/Data/Massiv/Core/IndexSpec.hs
--- a/tests/Data/Massiv/Core/IndexSpec.hs
+++ b/tests/Data/Massiv/Core/IndexSpec.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE GADTs               #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeOperators       #-}
+{-# LANGUAGE CPP                 #-}
 module Data.Massiv.Core.IndexSpec (Sz(..), SzZ(..), SzIx(..), DimIx(..), spec) where
 
 import           Control.Monad
@@ -10,6 +11,9 @@
 import           Data.Functor.Identity
 import           Test.Hspec
 import           Test.QuickCheck
+#if !MIN_VERSION_QuickCheck(2,10,0)
+import           Test.QuickCheck.Function
+#endif
 
 -- | Size that will result in a non-empty array
 newtype Sz ix = Sz ix deriving Show
@@ -98,6 +102,19 @@
 
 instance CoArbitrary Ix5 where
   coarbitrary (i :> ix) = coarbitrary i . coarbitrary ix
+
+instance Function Ix2 where
+  function = functionMap fromIx2 toIx2
+
+instance Function Ix3 where
+  function = functionMap fromIx3 toIx3
+
+instance Function Ix4 where
+  function = functionMap fromIx4 toIx4
+
+instance Function Ix5 where
+  function = functionMap fromIx5 toIx5
+
 
 
 -- instance Arbitrary Ix2 where
diff --git a/tests/Spec.hs b/tests/Spec.hs
--- a/tests/Spec.hs
+++ b/tests/Spec.hs
@@ -4,7 +4,7 @@
 import           Data.Massiv.Array.DelayedSpec          as Delayed
 import           Data.Massiv.Array.Manifest.VectorSpec  as Vector
 import           Data.Massiv.Array.ManifestSpec         as Manifest
---import           Data.Massiv.Array.MutableSpec         as Mutable
+import           Data.Massiv.Array.MutableSpec         as Mutable
 import           Data.Massiv.Array.Ops.ConstructSpec    as Construct
 import           Data.Massiv.Array.Ops.FoldSpec         as Fold
 import           Data.Massiv.Array.Ops.MapSpec          as Map
@@ -35,6 +35,6 @@
     describe "Delayed" $ Delayed.spec
     describe "Windowed" $ Windowed.spec
     describe "Manifest" $ Manifest.spec
-    --describe "Mutable" $ Mutable.spec
+    describe "Mutable" $ Mutable.spec
     describe "Stencil" $ Stencil.spec
     describe "Vector" $ Vector.spec