diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,11 @@
+# 0.6.1
+
+* Addition of `withLoadMArray_`, `withLoadMArrayS`, `withLoadMArrayS_`,
+  `withLoadMArrayST`, `withLoadMArrayST_`
+* Addition of `replaceSlice` and `replaceOuterSlice`
+* Addition of `quicksortBy`, `quicksortByM` and `quicksortByM_`
+* Fix performance regression for `quicksort` and `quicksortM_` introduced in previous release.
+
 # 0.6.0
 
 * Fix semantics of `Applicative`, `Num` and `Fractional` instance for `D` arrays:
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.6.0.0
+version:             0.6.1.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
diff --git a/src/Data/Massiv/Array.hs b/src/Data/Massiv/Array.hs
--- a/src/Data/Massiv/Array.hs
+++ b/src/Data/Massiv/Array.hs
@@ -133,6 +133,8 @@
   -- * Algorithms
   -- ** Sorting
   , quicksort
+  , quicksortBy
+  , quicksortByM
   , tally
   -- ** Iterations
   , iterateUntil
@@ -162,7 +164,7 @@
 import Data.Massiv.Array.Ops.Fold
 import Data.Massiv.Array.Ops.Map
 import Data.Massiv.Array.Ops.Slice
-import Data.Massiv.Array.Ops.Sort (quicksort, tally)
+import Data.Massiv.Array.Ops.Sort
 import Data.Massiv.Array.Ops.Transform
 import Data.Massiv.Array.Stencil
 import Data.Massiv.Core
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
@@ -91,6 +91,10 @@
   rnf (PArray c sz o a) = c `deepseq` sz `deepseq` o `seq` a `seq` ()
   {-# INLINE rnf #-}
 
+instance NFData ix => NFData (MArray s P ix e) where
+  rnf (MPArray sz _o _mb) = sz `deepseq` ()
+  {-# INLINE rnf #-}
+
 instance (Prim e, Eq e, Index ix) => Eq (Array P ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
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
@@ -79,6 +79,10 @@
   rnf (SArray c sz v) = c `deepseq` sz `deepseq` v `deepseq` ()
   {-# INLINE rnf #-}
 
+instance NFData ix => NFData (MArray s S ix e) where
+  rnf (MSArray sz mv) = sz `deepseq` mv `deepseq` ()
+  {-# INLINE rnf #-}
+
 instance (Storable e, Eq e, Index ix) => Eq (Array S ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
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
@@ -58,6 +58,9 @@
   rnf (UArray c sz v) = c `deepseq` sz `deepseq` v `deepseq` ()
   {-# INLINE rnf #-}
 
+instance NFData ix => NFData (MArray s U ix e) where
+  rnf (MUArray sz mv) = sz `deepseq` mv `deepseq` ()
+  {-# INLINE rnf #-}
 
 instance (VU.Unbox e, Index ix) => Construct U ix e where
   setComp c arr = arr { uComp = c }
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
@@ -81,10 +81,15 @@
   -- *** Modify
   , withMArray
   , withMArray_
+  , withLoadMArray_
   , withMArrayS
+  , withLoadMArrayS
   , withMArrayS_
+  , withLoadMArrayS_
   , withMArrayST
+  , withLoadMArrayST
   , withMArrayST_
+  , withLoadMArrayST_
   -- *** Initialize
   , initialize
   , initializeNew
@@ -851,6 +856,28 @@
 {-# INLINE withMArray_ #-}
 
 
+-- | Same as `withMArray_`, but the array supplied to this function can be any loadable
+-- array. For that reason it will be faster if supplied array is delayed.
+--
+-- @since 0.6.1
+withLoadMArray_ ::
+     forall r ix e r' m b. (Load r' ix e, Mutable r ix e, MonadUnliftIO m)
+  => Array r' ix e
+  -> (Scheduler m () -> MArray RealWorld r ix e -> m b)
+  -> m (Array r ix e)
+withLoadMArray_ arr action = do
+  marr <- liftIO $ unsafeNew (size arr)
+  withScheduler_ (getComp arr) $ \scheduler -> do
+    runBatch_ scheduler $ \_ -> loadArrayM scheduler arr (\i -> liftIO . unsafeLinearWrite marr i)
+    action scheduler marr
+  liftIO $ unsafeFreeze (getComp arr) marr
+{-# INLINE[2] withLoadMArray_ #-}
+{-# RULES
+"withLoadMArray_/withMArray_" [~2] withLoadMArray_ = withMArray_
+"withLoadMArrayS/withMArrayS" [~2] withLoadMArrayS = withMArrayS
+"withLoadMArrayS_/withMArrayS_" [~2] withLoadMArrayS_ = withMArrayS_
+#-}
+
 -- | Create a copy of a pure array, mutate it in place and return its frozen version. The important
 -- benefit over doing a manual `thawS` followed by a `freezeS` is that an array will only be copied
 -- once.
@@ -868,7 +895,7 @@
 {-# INLINE withMArrayS #-}
 
 
--- | Same as `withMArrayS`, but discards rhe element produced by the supplied action
+-- | Same as `withMArrayS`, except it discards the value produced by the supplied action
 --
 -- @since 0.5.0
 withMArrayS_ ::
@@ -880,6 +907,33 @@
 {-# INLINE withMArrayS_ #-}
 
 
+-- | Same as `withMArrayS`, but will work with any loadable array.
+--
+-- @since 0.6.1
+withLoadMArrayS ::
+     forall r ix e r' m a. (Load r' ix e, Mutable r ix e, PrimMonad m)
+  => Array r' ix e
+  -> (MArray (PrimState m) r ix e -> m a)
+  -> m (a, Array r ix e)
+withLoadMArrayS arr action = do
+  marr <- unsafeNew (size arr)
+  loadArrayM trivialScheduler_ arr (unsafeLinearWrite marr)
+  a <- action marr
+  (,) a <$> unsafeFreeze (getComp arr) marr
+{-# INLINE[2] withLoadMArrayS #-}
+
+-- | Same as `withMArrayS_`, but will work with any loadable array.
+--
+-- @since 0.6.1
+withLoadMArrayS_ ::
+     forall r ix e r' m a. (Load r' ix e, Mutable r ix e, PrimMonad m)
+  => Array r' ix e
+  -> (MArray (PrimState m) r ix e -> m a)
+  -> m (Array r ix e)
+withLoadMArrayS_ arr action = snd <$> withLoadMArrayS arr action
+{-# INLINE[2] withLoadMArrayS_ #-}
+
+
 -- | Same as `withMArrayS` but in `ST`. This is not only pure, but also the safest way to do
 -- mutation to the array.
 --
@@ -900,6 +954,29 @@
      Mutable r ix e => Array r ix e -> (forall s. MArray s r ix e -> ST s a) -> Array r ix e
 withMArrayST_ arr f = runST $ withMArrayS_ arr f
 {-# INLINE withMArrayST_ #-}
+
+
+-- | Same as `withMArrayST`, but works with any loadable array.
+--
+-- @since 0.6.1
+withLoadMArrayST ::
+     forall r ix e r' a. (Load r' ix e, Mutable r ix e)
+  => Array r' ix e
+  -> (forall s. MArray s r ix e -> ST s a)
+  -> (a, Array r ix e)
+withLoadMArrayST arr f = runST $ withLoadMArrayS arr f
+{-# INLINE[2] withLoadMArrayST #-}
+
+-- | Same as `withMArrayST_`, but works with any loadable array.
+--
+-- @since 0.6.1
+withLoadMArrayST_ ::
+     forall r ix e r' a. (Load r' ix e, Mutable r ix e)
+  => Array r' ix e
+  -> (forall s. MArray s r ix e -> ST s a)
+  -> Array r ix e
+withLoadMArrayST_ arr f = runST $ withLoadMArrayS_ arr f
+{-# INLINE[2] withLoadMArrayST_ #-}
 
 
 -- | /O(1)/ - Lookup an element in the mutable array. Returns `Nothing` when index is out of bounds.
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
@@ -10,6 +10,7 @@
 --
 module Data.Massiv.Array.Mutable.Algorithms
   ( quicksortM_
+  , quicksortByM_
   , unstablePartitionM
   , iterateUntilM
   ) where
@@ -36,7 +37,7 @@
 -- @since 0.3.2
 unstablePartitionM ::
      forall r e m. (Mutable r Ix1 e, PrimMonad m)
-  => MArray (PrimState m) r Ix1 e
+  => MVector (PrimState m) r e
   -> (e -> Bool) -- ^ Predicate
   -> m Ix1
 unstablePartitionM marr f = unsafeUnstablePartitionRegionM marr f 0 (unSz (msize marr) - 1)
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
@@ -12,10 +12,14 @@
 module Data.Massiv.Array.Ops.Sort
   ( tally
   , quicksort
+  , quicksortBy
+  , quicksortByM
   , quicksortM_
+  , quicksortByM_
   , unsafeUnstablePartitionRegionM
   ) where
 
+import Control.Monad.IO.Unlift
 import Control.Monad (when)
 import Control.Scheduler
 import Data.Massiv.Array.Delayed.Stream
@@ -25,7 +29,7 @@
 import Data.Massiv.Vector (scatMaybes, sunfoldrN)
 import System.IO.Unsafe
 
--- | Count how many occurance of each element there is in the array. Results will be
+-- | Count number of occurrences of each element in the array. Results will be
 -- sorted in ascending order of the element.
 --
 -- ==== __Example__
@@ -40,7 +44,7 @@
 --   [ (1,1), (2,3), (3,1), (4,2), (5,1) ]
 --
 -- @since 0.4.4
-tally :: (Mutable r Ix1 e, Resize r ix, Load r ix e, Ord e) => Array r ix e -> Array DS Ix1 (e, Int)
+tally :: (Mutable r Ix1 e, Resize r ix, Load r ix e, Ord e) => Array r ix e -> Vector DS (e, Int)
 tally arr
   | isEmpty arr = setComp (getComp arr) empty
   | otherwise = scatMaybes $ sunfoldrN (sz + 1) count (0, 0, sorted ! 0)
@@ -58,36 +62,48 @@
 {-# INLINE tally #-}
 
 
-
--- | Partition a segment of a vector. Starting and ending indices are unchecked.
---
--- @since 0.3.2
-unsafeUnstablePartitionRegionM ::
+unsafeUnstablePartitionRegionM' ::
      forall r e m. (Mutable r Ix1 e, PrimMonad m)
   => MArray (PrimState m) r Ix1 e
-  -> (e -> Bool)
+  -> (e -> m Bool)
   -> Ix1 -- ^ Start index of the region
   -> Ix1 -- ^ End index of the region
   -> m Ix1
-unsafeUnstablePartitionRegionM marr f start end = fromLeft start (end + 1)
+unsafeUnstablePartitionRegionM' marr f start end = fromLeft start (end + 1)
   where
     fromLeft i j
       | i == j = pure i
       | otherwise = do
-        x <- unsafeRead marr i
-        if f x
+        e <- f =<< unsafeLinearRead marr i
+        if e
           then fromLeft (i + 1) j
           else fromRight i (j - 1)
     fromRight i j
       | i == j = pure i
       | otherwise = do
-        x <- unsafeRead marr j
-        if f x
+        x <- unsafeLinearRead marr j
+        e <- f x
+        if e
           then do
-            unsafeWrite marr j =<< unsafeRead marr i
-            unsafeWrite marr i x
+            unsafeLinearWrite marr j =<< unsafeLinearRead marr i
+            unsafeLinearWrite marr i x
             fromLeft (i + 1) j
           else fromRight i (j - 1)
+{-# INLINE unsafeUnstablePartitionRegionM' #-}
+
+
+-- TODO: Replace `unsafeUnstablePartitionRegionM` with `unsafeUnstablePartitionRegionM'`
+-- | Partition a segment of a vector. Starting and ending indices are unchecked.
+--
+-- @since 0.3.2
+unsafeUnstablePartitionRegionM ::
+     forall r e m. (Mutable r Ix1 e, PrimMonad m)
+  => MVector (PrimState m) r e
+  -> (e -> Bool)
+  -> Ix1 -- ^ Start index of the region
+  -> Ix1 -- ^ End index of the region
+  -> m Ix1
+unsafeUnstablePartitionRegionM marr f = unsafeUnstablePartitionRegionM' marr (pure . f)
 {-# INLINE unsafeUnstablePartitionRegionM #-}
 
 
@@ -104,39 +120,81 @@
 {-# INLINE quicksort #-}
 
 
+-- | Same as `quicksortBy`, but instead of `Ord` constraint expects a custom `Ordering`.
+--
+-- @since 0.6.1
+quicksortByM ::
+     (Mutable r Ix1 e, MonadUnliftIO m) => (e -> e -> m Ordering) -> Vector r e -> m (Vector r e)
+quicksortByM f arr = withRunInIO $ \run -> withMArray_ arr (quicksortByM_ (\x y -> run (f x y)))
+{-# INLINE quicksortByM #-}
 
+-- | Same as `quicksortBy`, but instead of `Ord` constraint expects a custom `Ordering`.
+--
+-- @since 0.6.1
+quicksortBy ::
+     (Mutable r Ix1 e) => (e -> e -> Ordering) -> Vector r e -> Vector r e
+quicksortBy f arr =
+  unsafePerformIO $ withMArray_ arr (quicksortByM_ (\x y -> pure $ f x y))
+{-# INLINE quicksortBy #-}
+
 -- | Mutable version of `quicksort`
 --
 -- @since 0.3.2
 quicksortM_ ::
      (Ord e, Mutable r Ix1 e, PrimMonad m)
   => Scheduler m ()
-  -> MArray (PrimState m) r Ix1 e
+  -> MVector (PrimState m) r e
   -> m ()
-quicksortM_ scheduler marr =
+quicksortM_ = quicksortInternalM_ (\e1 e2 -> pure $ e1 < e2) (\e1 e2 -> pure $ e1 == e2)
+{-# INLINE quicksortM_ #-}
+
+
+-- | Same as `quicksortM_`, but instead of `Ord` constraint expects a custom `Ordering`.
+--
+-- @since 0.6.1
+quicksortByM_ ::
+     (Mutable r Ix1 e, PrimMonad m)
+  => (e -> e -> m Ordering)
+  -> Scheduler m ()
+  -> MVector (PrimState m) r e
+  -> m ()
+quicksortByM_ compareM =
+  quicksortInternalM_ (\x y -> (LT ==) <$> compareM x y) (\x y -> (EQ ==) <$> compareM x y)
+{-# INLINE quicksortByM_ #-}
+
+
+quicksortInternalM_ ::
+     (Mutable r Ix1 e, PrimMonad m)
+  => (e -> e -> m Bool)
+  -> (e -> e -> m Bool)
+  -> Scheduler m ()
+  -> MVector (PrimState m) r e
+  -> m ()
+quicksortInternalM_ fLT fEQ scheduler marr =
   scheduleWork scheduler $ qsort (numWorkers scheduler) 0 (unSz (msize marr) - 1)
   where
-    leSwap i j = do
-      ei <- unsafeRead marr i
-      ej <- unsafeRead marr j
-      if ei < ej
+    ltSwap i j = do
+      ei <- unsafeLinearRead marr i
+      ej <- unsafeLinearRead marr j
+      lt <- fLT ei ej
+      if lt
         then do
-          unsafeWrite marr i ej
-          unsafeWrite marr j ei
+          unsafeLinearWrite marr i ej
+          unsafeLinearWrite marr j ei
           pure ei
         else pure ej
-    {-# INLINE leSwap #-}
+    {-# INLINE ltSwap #-}
     getPivot lo hi = do
       let !mid = (hi + lo) `div` 2
-      _ <- leSwap mid lo
-      _ <- leSwap hi lo
-      leSwap mid hi
+      _ <- ltSwap mid lo
+      _ <- ltSwap hi lo
+      ltSwap mid hi
     {-# INLINE getPivot #-}
     qsort !n !lo !hi =
       when (lo < hi) $ do
         p <- getPivot lo hi
-        l <- unsafeUnstablePartitionRegionM marr (< p) lo (hi - 1)
-        h <- unsafeUnstablePartitionRegionM marr (== p) l hi
+        l <- unsafeUnstablePartitionRegionM' marr (`fLT` p) lo (hi - 1)
+        h <- unsafeUnstablePartitionRegionM' marr (`fEQ` p) l hi
         if n > 0
           then do
             let !n' = n - 1
@@ -145,4 +203,4 @@
           else do
             qsort n lo (l - 1)
             qsort n h hi
-{-# INLINE quicksortM_ #-}
+{-# INLINE quicksortInternalM_ #-}
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
@@ -49,6 +49,8 @@
   , splitAtM
   , splitAt'
   , splitExtractM
+  , replaceSlice
+  , replaceOuterSlice
   -- ** Upsample/Downsample
   , upsample
   , downsample
@@ -748,6 +750,96 @@
   rightArr <- extractFromToM rightArrStartIx szIx arr
   pure (leftArr, midArr, rightArr)
 {-# INLINE splitExtractM #-}
+
+
+
+-- | Replace a slice of an array with another one
+--
+-- ==== __Example__
+--
+-- >>> import Data.Massiv.Array
+-- >>> arr = makeArrayR U Seq (Sz3 3 4 5) fromIx3
+-- >>> arr' = makeArrayR U Seq (Sz3 3 4 5) (fromIx3 . liftIndex (* 100))
+-- >>> replaceSlice 2 1 (arr' <!> (2, 3)) arr
+-- Array DL Seq (Sz (3 :> 4 :. 5))
+--   [ [ [ (0,0,0), (0,0,1), (0,0,2), (0,0,3), (0,0,4) ]
+--     , [ (0,300,0), (0,300,100), (0,300,200), (0,300,300), (0,300,400) ]
+--     , [ (0,2,0), (0,2,1), (0,2,2), (0,2,3), (0,2,4) ]
+--     , [ (0,3,0), (0,3,1), (0,3,2), (0,3,3), (0,3,4) ]
+--     ]
+--   , [ [ (1,0,0), (1,0,1), (1,0,2), (1,0,3), (1,0,4) ]
+--     , [ (100,300,0), (100,300,100), (100,300,200), (100,300,300), (100,300,400) ]
+--     , [ (1,2,0), (1,2,1), (1,2,2), (1,2,3), (1,2,4) ]
+--     , [ (1,3,0), (1,3,1), (1,3,2), (1,3,3), (1,3,4) ]
+--     ]
+--   , [ [ (2,0,0), (2,0,1), (2,0,2), (2,0,3), (2,0,4) ]
+--     , [ (200,300,0), (200,300,100), (200,300,200), (200,300,300), (200,300,400) ]
+--     , [ (2,2,0), (2,2,1), (2,2,2), (2,2,3), (2,2,4) ]
+--     , [ (2,3,0), (2,3,1), (2,3,2), (2,3,3), (2,3,4) ]
+--     ]
+--   ]
+--
+-- @since 0.6.1
+replaceSlice ::
+     ( MonadThrow m
+     , Extract r ix e
+     , Source (R r) ix e
+     , Load (R r) (Lower ix) e
+     , Resize (R r) (Lower ix)
+     )
+  => Dim
+  -> Ix1
+  -> Array (R r) (Lower ix) e
+  -> Array r ix e
+  -> m (Array DL ix e)
+replaceSlice dim i sl arr = do
+  (l, m, r) <- splitExtractM dim i (SafeSz 1) arr
+  m' <- resizeM (size m) sl
+  concatM dim [l, m', r]
+{-# INLINE replaceSlice #-}
+
+
+-- | Replace an outer slice of an array with another one
+--
+-- ==== __Example__
+--
+-- >>> import Data.Massiv.Array
+-- >>> arr = makeArrayR U Seq (Sz3 3 4 5) fromIx3
+-- >>> arr' = makeArrayR U Seq (Sz3 3 4 5) (fromIx3 . liftIndex (* 100))
+-- >>> replaceOuterSlice 1 (arr' !> 2) arr
+-- Array DL Seq (Sz (3 :> 4 :. 5))
+--   [ [ [ (0,0,0), (0,0,1), (0,0,2), (0,0,3), (0,0,4) ]
+--     , [ (0,1,0), (0,1,1), (0,1,2), (0,1,3), (0,1,4) ]
+--     , [ (0,2,0), (0,2,1), (0,2,2), (0,2,3), (0,2,4) ]
+--     , [ (0,3,0), (0,3,1), (0,3,2), (0,3,3), (0,3,4) ]
+--     ]
+--   , [ [ (200,0,0), (200,0,100), (200,0,200), (200,0,300), (200,0,400) ]
+--     , [ (200,100,0), (200,100,100), (200,100,200), (200,100,300), (200,100,400) ]
+--     , [ (200,200,0), (200,200,100), (200,200,200), (200,200,300), (200,200,400) ]
+--     , [ (200,300,0), (200,300,100), (200,300,200), (200,300,300), (200,300,400) ]
+--     ]
+--   , [ [ (2,0,0), (2,0,1), (2,0,2), (2,0,3), (2,0,4) ]
+--     , [ (2,1,0), (2,1,1), (2,1,2), (2,1,3), (2,1,4) ]
+--     , [ (2,2,0), (2,2,1), (2,2,2), (2,2,3), (2,2,4) ]
+--     , [ (2,3,0), (2,3,1), (2,3,2), (2,3,3), (2,3,4) ]
+--     ]
+--   ]
+--
+-- @since 0.6.1
+replaceOuterSlice ::
+     ( MonadThrow m
+     , Extract r ix e
+     , Source (R r) ix e
+     , Load (R r) (Lower ix) e
+     , Resize (R r) (Lower ix)
+     )
+  => Ix1
+  -> Array (R r) (Lower ix) e
+  -> Array r ix e
+  -> m (Array DL ix e)
+replaceOuterSlice i sl arr = replaceSlice (dimensions (size arr)) i sl arr
+{-# INLINE replaceOuterSlice #-}
+
 
 -- | Delete a region from an array along the specified dimension.
 --
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
@@ -38,6 +38,10 @@
   , WorkerStates
   , initWorkerStates
   , module Data.Massiv.Core.Index
+  -- * Numeric
+  , FoldNumeric
+  , Numeric
+  , NumericFloat
   -- * Exceptions
   , MonadThrow(..)
   , throw
@@ -56,6 +60,7 @@
 import Data.Massiv.Core.Index
 import Data.Massiv.Core.List
 import Data.Massiv.Core.Exception
+import Data.Massiv.Core.Operations (FoldNumeric, Numeric, NumericFloat)
 
 
 -- | Append computation strategy using `Comp`'s `Monoid` instance.
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
@@ -565,7 +565,7 @@
 -- @since 0.1.0
 unsafeRead :: (Mutable r ix e, PrimMonad m) =>
                MArray (PrimState m) r ix e -> ix -> m e
-unsafeRead !marr !ix = unsafeLinearRead marr (toLinearIndex (msize marr) ix)
+unsafeRead marr = unsafeLinearRead marr . toLinearIndex (msize marr)
 {-# INLINE unsafeRead #-}
 
 -- | Write an element into array
@@ -573,7 +573,7 @@
 -- @since 0.1.0
 unsafeWrite :: (Mutable r ix e, PrimMonad m) =>
                MArray (PrimState m) r ix e -> ix -> e -> m ()
-unsafeWrite !marr !ix = unsafeLinearWrite marr (toLinearIndex (msize marr) ix)
+unsafeWrite marr = unsafeLinearWrite marr . toLinearIndex (msize marr)
 {-# INLINE unsafeWrite #-}
 
 
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
@@ -429,7 +429,7 @@
 --
 -- /Related/: `head'`, `shead'`, `sheadM`, `unconsM`.
 --
--- /__Throws Exceptions__/: `SizeEmptyException`
+-- /__Throws Exceptions__/: `SizeEmptyException` when array is empty
 --
 -- ==== __Examples__
 --
@@ -1410,7 +1410,7 @@
 {-# INLINE sfilter #-}
 
 
--- | Similar to `sfilter`, but map with an index aware function.
+-- | Similar to `sfilter`, but filter with an index aware function.
 --
 -- ==== __Examples__
 --
@@ -1462,12 +1462,10 @@
 {-# INLINE sfilterM #-}
 
 
--- | Similar to `filterM`, but map with an index aware function.
+-- | Similar to `filterM`, but filter with an index aware function.
 --
 -- Corresponds to: @`filterM` (uncurry f) . `simap` (,)@
 --
--- ==== __Examples__
---
 -- @since 0.5.0
 sifilterM ::
      (Stream r ix a, Applicative f) => (ix -> a -> f Bool) -> Array r ix a -> f (Vector DS a)
@@ -2458,63 +2456,109 @@
 
 
 
--- |
+-- | Add all elements of the array together
 --
+-- /Related/: `sum`.
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.ssum $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 88
+--
 -- @since 0.5.0
 ssum :: (Num e, Stream r ix e) => Array r ix e -> e
 ssum = sfoldl (+) 0
 {-# INLINE ssum #-}
 
--- |
+-- | Multiply all elements of the array together
 --
+-- /Related/: `product`.
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.sproduct $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 10500
+--
 -- @since 0.5.0
 sproduct :: (Num e, Stream r ix e) => Array r ix e -> e
 sproduct = sfoldl (*) 1
 {-# INLINE sproduct #-}
 
 
--- |
+-- | /O(n)/ - Find the largest value in the array. Throws an error on empty.
 --
+-- /Related/: `smaximumM`, `maximum`, `maximumM`.
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.smaximum' $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 70
+-- >>> V.smaximum' (V.empty :: Vector D Int)
+-- *** Exception: SizeEmptyException: (Sz1 0) corresponds to an empty array
+--
 -- @since 0.5.0
 smaximum' :: (Ord e, Stream r ix e) => Array r ix e -> e
 smaximum' = sfoldl1' max
 {-# INLINE smaximum' #-}
 
--- |
+-- | /O(n)/ - Find the largest value in the array.
 --
+-- /Related/: `smaximum`, `maximum`, `maximumM`.
+--
+-- /__Throws Exceptions__/: `SizeEmptyException` when array is empty
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.smaximumM $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 70
+-- >>> V.smaximumM (V.empty :: Vector D Int) :: Maybe Int
+-- Nothing
+--
 -- @since 0.5.0
 smaximumM :: (Ord e, Stream r ix e, MonadThrow m) => Array r ix e -> m e
 smaximumM = sfoldl1M (\e acc -> pure (max e acc))
 {-# INLINE smaximumM #-}
 
 
-
--- |
+-- | /O(n)/ - Find the smallest value in the array. Throws an error on empty.
 --
+-- /Related/: `sminimumM`, `minimum`, `minimumM`.
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.sminimum' $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 3
+-- >>> V.sminimum' (V.empty :: Array D Ix2 Int)
+-- *** Exception: SizeEmptyException: (Sz (0 :. 0)) corresponds to an empty array
+--
 -- @since 0.5.0
 sminimum' :: (Ord e, Stream r ix e) => Array r ix e -> e
 sminimum' = sfoldl1' min
 {-# INLINE sminimum' #-}
 
--- |
+-- | /O(n)/ - Find the smallest value in the array.
 --
+-- /Related/: `sminimum'`, `minimum`, `minimumM`.
+--
+-- /__Throws Exceptions__/: `SizeEmptyException` when array is empty
+--
 -- ==== __Examples__
 --
+-- >>> import Data.Massiv.Vector as V
+-- >>> V.sminimumM $ V.sfromList [10, 3, 70, 5 :: Int]
+-- 3
+-- >>> V.sminimumM (V.empty :: Array D Ix2 Int) :: Maybe Int
+-- Nothing
+--
 -- @since 0.5.0
 sminimumM :: (Ord e, Stream r ix e, MonadThrow m) => Array r ix e -> m e
 sminimumM = sfoldl1M (\e acc -> pure (min e acc))
 {-# INLINE sminimumM #-}
-
 
 
 -- | See `stake`.
diff --git a/tests/doctests.hs b/tests/doctests.hs
--- a/tests/doctests.hs
+++ b/tests/doctests.hs
@@ -1,6 +1,17 @@
+{-# LANGUAGE CPP #-}
 module Main where
 
+#if __GLASGOW_HASKELL__ >= 802 && __GLASGOW_HASKELL__ < 810
+
 import Test.DocTest (doctest)
 
 main :: IO ()
 main = doctest ["-Iinclude","src"]
+
+#else
+
+-- TODO: fix doctest support
+main :: IO ()
+main = putStrLn "\nDoctests are not supported for older ghc version\n"
+
+#endif