diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,10 @@
+# 0.5.7
+
+* Improve performance of `><.` and `><!` while making their constraints a bit more relaxed.
+* Expose `unsafeLoadIntoM` and `unsafeLoadIntoS`
+* Expose `eqArrays` and `compareArrays`
+* Add `multiplyMatrixByVector` and `multiplyVectorByMatrix`
+
 # 0.5.6
 
 * Fix `(-.)` (it was incorrectly implemented as a flip of `(.-)`
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.5.6.0
+version:             0.5.7.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/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
@@ -18,8 +18,8 @@
   ( D(..)
   , Array(..)
   , delay
-  , eq
-  , ord
+  , eqArrays
+  , compareArrays
   , imap
   ) where
 
@@ -99,11 +99,11 @@
 
 
 instance (Eq e, Index ix) => Eq (Array D ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Ord e, Index ix) => Ord (Array D ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 instance Functor (Array D ix) where
@@ -218,6 +218,13 @@
 
 
 instance Num e => Numeric D e where
+  unsafeDotProduct a1 a2 = go 0 0
+    where
+      !len = totalElem (size a1)
+      go !acc i
+        | i < len = go (acc + unsafeLinearIndex a1 i * unsafeLinearIndex a2 i) (i + 1)
+        | otherwise = acc
+  {-# INLINE unsafeDotProduct #-}
   foldArray f !initAcc arr = go initAcc 0
     where
       !len = totalElem (dSize arr)
@@ -247,26 +254,30 @@
  #-}
 
 -- | /O(min (n1, n2))/ - Compute array equality by applying a comparing function to each element.
-eq :: (Source r1 ix e1, Source r2 ix e2) =>
+--
+-- @since 0.5.7
+eqArrays :: (Source r1 ix e1, Source r2 ix e2) =>
       (e1 -> e2 -> Bool) -> Array r1 ix e1 -> Array r2 ix e2 -> Bool
-eq f arr1 arr2 =
+eqArrays f arr1 arr2 =
   (size arr1 == size arr2) &&
   not (A.any not
        (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
            f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix)))
-{-# INLINE eq #-}
+{-# INLINE eqArrays #-}
 
 -- | /O(min (n1, n2))/ - Compute array ordering by applying a comparing function to each element.
 -- The exact ordering is unspecified so this is only intended for use in maps and the like where
 -- you need an ordering but do not care about which one is used.
-ord :: (Source r1 ix e1, Source r2 ix e2) =>
+--
+-- @since 0.5.7
+compareArrays :: (Source r1 ix e1, Source r2 ix e2) =>
        (e1 -> e2 -> Ordering) -> Array r1 ix e1 -> Array r2 ix e2 -> Ordering
-ord f arr1 arr2 =
+compareArrays f arr1 arr2 =
   compare (size arr1) (size arr2) <>
   A.fold
     (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
        f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
-{-# INLINE ord #-}
+{-# INLINE compareArrays #-}
 
 
 -- -- | The usual map.
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
@@ -204,8 +204,8 @@
     S.unstreamIntoM marr (stepsSize sts) (stepsStream sts)
   {-# INLINE unsafeLoadIntoS #-}
 
-  unsafeLoadInto marr arr = liftIO $ unsafeLoadIntoS marr arr
-  {-# INLINE unsafeLoadInto #-}
+  unsafeLoadIntoM marr arr = liftIO $ unsafeLoadIntoS marr arr
+  {-# INLINE unsafeLoadIntoM #-}
 
 
 -- cons :: e -> Array DS Ix1 e -> Array DS Ix1 e
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
@@ -49,7 +49,6 @@
 import Control.Monad.Primitive
 import Control.Monad.ST (runST)
 import qualified Data.Foldable as F (Foldable(..))
-import Data.Massiv.Array.Delayed.Pull (eq, ord)
 import Data.Massiv.Array.Delayed.Push (DL)
 import Data.Massiv.Array.Delayed.Stream (DS)
 import Data.Massiv.Array.Manifest.Internal (M, computeAs, toManifest)
@@ -112,11 +111,11 @@
   {-# INLINE rnf #-}
 
 instance (Index ix, Eq e) => Eq (Array B ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Index ix, Ord e) => Ord (Array B ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 instance Index ix => Construct B ix e where
@@ -284,11 +283,11 @@
   {-# INLINE rnf #-}
 
 instance (Index ix, NFData e, Eq e) => Eq (Array N ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Index ix, NFData e, Ord e) => Ord (Array N ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 
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
@@ -90,11 +90,11 @@
 
 
 instance (Eq e, Index ix) => Eq (Array M ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Ord e, Index ix) => Ord (Array M ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 
@@ -126,7 +126,6 @@
   {-# INLINE elem #-}
   toList arr = build (\ c n -> foldrFB c n arr)
   {-# INLINE toList #-}
-
 
 
 instance Index ix => Source M ix e where
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
@@ -87,11 +87,11 @@
   {-# INLINE rnf #-}
 
 instance (Prim e, Eq e, Index ix) => Eq (Array P ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Prim e, Ord e, Index ix) => Ord (Array P ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 instance (Prim e, Index ix) => Construct P ix e where
@@ -255,6 +255,13 @@
         | i < len = go (acc + unsafeLinearIndex a1 i * unsafeLinearIndex a2 i) (i + 1)
         | otherwise = acc
   {-# INLINE unsafeDotProduct #-}
+  powerSumArray arr p = go 0 0
+    where
+      !len = totalElem (size arr)
+      go !acc i
+        | i < len = go (acc + unsafeLinearIndex arr i ^ p) (i + 1)
+        | otherwise = acc
+  {-# INLINE powerSumArray #-}
   foldArray f !initAcc arr = go initAcc 0
     where
       !len = totalElem (size arr)
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
@@ -36,7 +36,7 @@
 import Control.DeepSeq (NFData(..), deepseq)
 import Control.Monad.IO.Unlift
 import Control.Monad.Primitive (unsafePrimToPrim)
-import Data.Massiv.Array.Delayed.Pull (eq, ord)
+import Data.Massiv.Array.Delayed.Pull (eqArrays, compareArrays)
 import Data.Massiv.Array.Manifest.Internal
 import Data.Massiv.Array.Manifest.Primitive (shrinkMutableByteArray)
 import Data.Primitive.ByteArray (MutableByteArray(..))
@@ -77,11 +77,11 @@
   {-# INLINE rnf #-}
 
 instance (Storable e, Eq e, Index ix) => Eq (Array S ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (Storable e, Ord e, Index ix) => Ord (Array S ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 instance (Storable e, Index ix) => Construct S ix e where
@@ -231,6 +231,13 @@
         | i < len = go (acc + unsafeLinearIndex a1 i * unsafeLinearIndex a2 i) (i + 1)
         | otherwise = acc
   {-# INLINE unsafeDotProduct #-}
+  powerSumArray arr p = go 0 0
+    where
+      !len = totalElem (size arr)
+      go !acc i
+        | i < len = go (acc + unsafeLinearIndex arr i ^ p) (i + 1)
+        | otherwise = acc
+  {-# INLINE powerSumArray #-}
   foldArray f !initAcc arr = go initAcc 0
     where
       !len = totalElem (size arr)
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
@@ -26,7 +26,7 @@
   ) where
 
 import Control.DeepSeq (NFData(..), deepseq)
-import Data.Massiv.Array.Delayed.Pull (eq, ord)
+import Data.Massiv.Array.Delayed.Pull (eqArrays, compareArrays)
 import Data.Massiv.Array.Manifest.Internal (M, toManifest)
 import Data.Massiv.Array.Manifest.List as A
 import Data.Massiv.Vector.Stream as S (steps, isteps)
@@ -68,11 +68,11 @@
 
 
 instance (VU.Unbox e, Eq e, Index ix) => Eq (Array U ix e) where
-  (==) = eq (==)
+  (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
 instance (VU.Unbox e, Ord e, Index ix) => Ord (Array U ix e) where
-  compare = ord compare
+  compare = compareArrays compare
   {-# INLINE compare #-}
 
 
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
@@ -280,7 +280,7 @@
 loadArray arr =
   liftIO $ do
     marr <- newMaybeInitialized arr
-    unsafeLoadInto marr arr
+    unsafeLoadIntoM marr arr
 {-# INLINE loadArray #-}
 
 
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
@@ -37,8 +37,10 @@
   -- ** Matrix multiplication
   , (.><)
   , (!><)
+  , multiplyMatrixByVector
   , (><.)
   , (><!)
+  , multiplyVectorByMatrix
   , (.><.)
   , (!><!)
   , multiplyMatrices
@@ -100,6 +102,7 @@
   , atan2A
   ) where
 
+import Data.Massiv.Array.Mutable
 import Data.Massiv.Array.Manifest
 import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Array.Delayed.Push
@@ -116,6 +119,7 @@
 import Control.Scheduler
 import Control.Monad (when)
 import qualified Data.Foldable as F
+import Data.Function
 
 infixr 8  .^, .^^
 infixl 7  !*!, .*., .*, *., !/!, ./., ./, /., `quotA`, `remA`, `divA`, `modA`
@@ -260,8 +264,10 @@
 {-# INLINE (.*.) #-}
 
 
--- | Multiplication of two arrays pointwise. Prefer to use monadic version of this
--- function `.*.` whenever possible, because it is better to avoid partial functions.
+-- | Multiplication of two arrays pointwise,
+-- i.e. <https://en.wikipedia.org/wiki/Hadamard_product_(matrices) Hadamard product>.
+-- Prefer to use monadic version of this function `.*.` whenever possible,
+-- because it is better to avoid partial functions.
 --
 -- [Partial] Mismatched array sizes will result in an impure exception being thrown.
 --
@@ -375,14 +381,6 @@
     comp = getComp v1 <> getComp v2
 {-# INLINE dotM #-}
 
--- | Compute L2 norm of an array.
---
--- @since 0.5.6
-normL2 :: (NumericFloat r e, Source r ix e) => Array r ix e -> e
-normL2 v
-  | getComp v == Seq = sqrt $! unsafeDotProduct v v
-  | otherwise = sqrt $! unsafePerformIO $ unsafeDotProductIO v v
-{-# INLINE normL2 #-}
 
 unsafeDotProductIO ::
      (MonadUnliftIO m, Numeric r b, Source r ix b)
@@ -409,16 +407,47 @@
 {-# INLINE unsafeDotProductIO #-}
 
 
+-- | Compute L2 norm of an array.
+--
+-- @since 0.5.6
+normL2 :: (Floating e, Numeric r e, Source r ix e) => Array r ix e -> e
+normL2 v
+  | getComp v == Seq = sqrt $! powerSumArray v 2
+  | otherwise = sqrt $! unsafePerformIO $ powerSumArrayIO v 2
+{-# INLINE normL2 #-}
+
+powerSumArrayIO ::
+     (MonadUnliftIO m, Numeric r b, Source r ix b)
+  => Array r ix b
+  -> Int
+  -> m b
+powerSumArrayIO v p = do
+  results <-
+    withScheduler (getComp v) $ \scheduler ->
+      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+        let n = SafeSz chunkLength
+        loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
+          scheduleWork scheduler $ pure $! powerSumArray (unsafeLinearSlice start n v) p
+        when (slackStart < totalLength) $ do
+          let k = SafeSz (totalLength - slackStart)
+          scheduleWork scheduler $ pure $! powerSumArray (unsafeLinearSlice slackStart k v) p
+  pure $! F.foldl' (+) 0 results
+  where
+    totalLength = totalElem (size v)
+{-# INLINE powerSumArrayIO #-}
+
+
 -- | Multiply a matrix by a column vector. Same as `!><` but produces monadic
 -- computation that allows for handling failure.
 --
 -- /__Throws Exception__/: `SizeMismatchException` when inner dimensions of arrays do not match.
 --
 -- @since 0.5.6
-(.><) :: (MonadThrow m, Numeric r e, Source r Ix1 e, Source r Ix2 e) =>
-         Matrix r e -- ^ Matrix
-      -> Vector r e -- ^ Column vector (Used many times, so make sure it is computed)
-      -> m (Vector D e)
+(.><) ::
+     (MonadThrow m, Numeric r e, Source r Ix1 e, Source r Ix2 e)
+  => Matrix r e -- ^ Matrix
+  -> Vector r e -- ^ Column vector (Used many times, so make sure it is computed)
+  -> m (Vector D e)
 (.><) mm v
   | mCols /= n = throwM $ SizeMismatchException (size mm) (Sz2 n 1)
   | otherwise = pure $ makeArray (getComp mm <> getComp v) (Sz1 mRows) $ \i ->
@@ -428,7 +457,20 @@
     sz@(Sz1 n) = size v
 {-# INLINE (.><) #-}
 
+-- | Multiply matrix by a column vector. Same as `.><` but returns computed version of a vector
+--
+-- /__Throws Exception__/: `SizeMismatchException` when inner dimensions of arrays do not match.
+--
+-- @since 0.5.7
+multiplyMatrixByVector ::
+     (MonadThrow m, Numeric r e, Mutable r Ix1 e, Mutable r Ix2 e)
+  => Matrix r e -- ^ Matrix
+  -> Vector r e -- ^ Column vector (Used many times, so make sure it is computed)
+  -> m (Vector r e)
+multiplyMatrixByVector mm v = compute <$> mm .>< v
+{-# INLINE multiplyMatrixByVector #-}
 
+
 -- | Multiply a matrix by a column vector
 --
 -- [Partial] Throws impure exception when inner dimensions do not agree
@@ -453,15 +495,44 @@
          Vector r e -- ^ Row vector
       -> Matrix r e -- ^ Matrix
       -> m (Vector D e)
-(><.) v mm
+(><.) v mm = delay <$> multiplyVectorByMatrix v mm
+{-# INLINE (><.) #-}
+
+-- | Multiply a row vector by a matrix. Same as `><.` but returns computed vector instead of
+-- a delayed one.
+--
+-- /__Throws Exception__/: `SizeMismatchException` when inner dimensions of arrays do not match.
+--
+-- @since 0.5.7
+multiplyVectorByMatrix ::
+     (MonadThrow m, Numeric r e, Mutable r Ix1 e, Mutable r Ix2 e)
+  => Vector r e -- ^ Row vector
+  -> Matrix r e -- ^ Matrix
+  -> m (Vector r e)
+multiplyVectorByMatrix v mm
   | mRows /= n = throwM $ SizeMismatchException (Sz2 1 n) (size mm)
-  | otherwise = pure $ makeArray (getComp mm <> getComp v) (Sz1 mCols) $ \i ->
-      unsafeDotProduct (unsafeLinearSlice (i * n) sz mm') v
+  | otherwise =
+    pure $!
+    unsafePerformIO $ do
+      mv <- new (Sz mCols)
+      withMassivScheduler_ comp $ \scheduler -> do
+        let loopCols x ivto =
+              fix $ \go im iv ->
+                when (iv < ivto) $ do
+                  _ <- unsafeLinearModify mv (\a -> pure $ a + unsafeLinearIndex mm im * x) iv
+                  go (im + 1) (iv + 1)
+            loopRows i0 from to =
+              flip fix i0 $ \go i ->
+                when (i < mRows) $ do
+                  loopCols (unsafeLinearIndex v i) to (i * mCols + from) from
+                  go (i + 1)
+        splitLinearlyM_ scheduler mCols (loopRows 0)
+      unsafeFreeze comp mv
   where
+    comp = getComp mm <> getComp v
     Sz2 mRows mCols = size mm
-    mm' = compute (transpose mm)
-    sz@(Sz1 n) = size v
-{-# INLINE (><.) #-}
+    Sz1 n = size v
+{-# INLINE multiplyVectorByMatrix #-}
 
 
 -- | Multiply a row vector by a matrix.
@@ -471,7 +542,7 @@
 -- @since 0.5.6
 (><!) ::
      (Numeric r e, Mutable r Ix1 e, Mutable r Ix2 e)
-  => Vector r e -- ^ Row vector
+  => Vector r e -- ^ Row vector (Used many times, so make sure it is computed)
   -> Matrix r e -- ^ Matrix
   -> Vector D e
 (><!) v mm = throwEither (v ><. mm)
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
@@ -36,6 +36,8 @@
   , all
   , any
   , elem
+  , eqArrays
+  , compareArrays
 
   -- ** Single dimension folds
   -- *** Safe inner most
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
@@ -36,6 +36,8 @@
   , unsafeThaw
   , unsafeFreeze
   , unsafeNew
+  , unsafeLoadIntoS
+  , unsafeLoadIntoM
   , unsafeCreateArray
   , unsafeCreateArray_
   , unsafeCreateArrayS
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
@@ -325,7 +325,7 @@
   -- | Load into a supplied mutable array sequentially. Returned array does not have to be
   -- the same
   --
-  -- @since 0.5.0
+  -- @since 0.5.7
   unsafeLoadIntoS ::
        (Mutable r' ix e, PrimMonad m)
     => MArray (PrimState m) r' ix e
@@ -337,17 +337,17 @@
 
   -- | Same as `unsafeLoadIntoS`, but respecting computation strategy.
   --
-  -- @since 0.5.0
-  unsafeLoadInto ::
+  -- @since 0.5.7
+  unsafeLoadIntoM ::
        (Mutable r' ix e, MonadIO m)
     => MArray RealWorld r' ix e
     -> Array r ix e
     -> m (MArray RealWorld r' ix e)
-  unsafeLoadInto marr arr = do
+  unsafeLoadIntoM marr arr = do
     liftIO $ withMassivScheduler_ (getComp arr) $ \scheduler ->
       loadArrayM scheduler arr (unsafeLinearWrite marr)
     pure marr
-  {-# INLINE unsafeLoadInto #-}
+  {-# INLINE unsafeLoadIntoM #-}
 
 -- | Selects an optimal scheduler for the supplied strategy, but it works only in `IO`
 withMassivScheduler_ :: Comp -> (Scheduler IO () -> IO ()) -> IO ()
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 @@
   , loopM_
   , loopDeepM
   , splitLinearly
+  , splitLinearlyM_
   , splitLinearlyWith_
   , splitLinearlyWithM_
   , splitLinearlyWithStartAtM_
@@ -21,6 +22,7 @@
   ) where
 
 import Control.Scheduler
+import Control.Monad
 
 -- | Efficient loop with an accumulator
 --
@@ -58,10 +60,22 @@
     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
@@ -108,6 +122,18 @@
     !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_ ::
+     Monad m => Scheduler m () -> 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_`.
@@ -126,11 +152,10 @@
 splitLinearlyWithM_ ::
      Monad m => Scheduler m () -> Int -> (Int -> m b) -> (Int -> b -> m c) -> m ()
 splitLinearlyWithM_ scheduler totalLength make write =
-  splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
-    loopM_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
-      scheduleWork_ scheduler $
-      loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> make k >>= write k
-    scheduleWork_ scheduler $ loopM_ slackStart (< totalLength) (+ 1) $ \ !k -> make k >>= write k
+  splitLinearlyM_ scheduler totalLength go
+  where
+    go start end = loopM_ start (< end) (+ 1) $ \ k -> make k >>= write k
+    {-# INLINE go #-}
 {-# INLINE splitLinearlyWithM_ #-}
 
 
@@ -144,8 +169,9 @@
     loopM_ startAt (< (slackStart + startAt)) (+ chunkLength) $ \ !start ->
       scheduleWork_ scheduler $
       loopM_ start (< (start + chunkLength)) (+ 1) $ \ !k -> make k >>= write k
-    scheduleWork_ scheduler $
-      loopM_ (slackStart + startAt) (< (totalLength + startAt)) (+ 1) $ \ !k -> make k >>= write k
+    when (slackStart < totalLength) $
+      scheduleWork_ scheduler $
+        loopM_ (slackStart + startAt) (< (totalLength + startAt)) (+ 1) $ \ !k -> make k >>= write k
 {-# INLINE splitLinearlyWithStartAtM_ #-}
 
 
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
@@ -1,4 +1,3 @@
-{-# LANGUAGE DefaultSignatures #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables #-}
@@ -16,7 +15,6 @@
   ) where
 
 import Data.Massiv.Core.Common
--- import Data.Massiv.Array.Ops.Fold.Internal
 
 
 class Num e => Numeric r e where
@@ -37,8 +35,12 @@
   productArray = foldArray (*) 1
   {-# INLINE productArray #-}
 
-  -- -- | Raise each element in the array to some non-negative power and sum the results
-  -- powerSumArray :: Array r Ix1 e -> Int -> e
+  -- | Raise each element in the array to some non-negative power and sum the results
+  --
+  -- @since 0.5.7
+  powerSumArray :: Index ix => Array r ix e -> Int -> e
+  powerSumArray arr p = sumArray $ unsafeLiftArray (^ p) arr
+  {-# INLINE powerSumArray #-}
 
   -- | Compute dot product without any extraneous checks
   --
@@ -80,6 +82,9 @@
   multiplicationPointwise = unsafeLiftArray2 (*)
   {-# INLINE multiplicationPointwise #-}
 
+  -- TODO:
+  --  - rename to powerScalar
+  --  - add? powerPointwise :: Array r ix e -> Array r ix Int -> Array r ix e
   -- | Raise each element of the array to the power
   powerPointwise :: Index ix => Array r ix e -> Int -> Array r ix e
   powerPointwise arr pow = unsafeLiftArray (^ pow) arr