Numeric.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}

-- |
-- Module      : Data.Massiv.Array.Numeric
-- Copyright   : (c) Alexey Kuleshevich 2018-2022
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
module Data.Massiv.Array.Numeric (
  -- * Numeric
  Numeric,
  NumericFloat,
  liftNumArray2M,

  -- ** Pointwise addition
  (.+),
  (+.),
  (.+.),
  (!+!),
  sumArraysM,
  sumArrays',

  -- ** Pointwise subtraction
  (.-),
  (-.),
  (.-.),
  (!-!),

  -- ** Pointwise multiplication
  (.*),
  (*.),
  (.*.),
  (!*!),
  (.^),
  productArraysM,
  productArrays',

  -- ** Dot product
  (!.!),
  dotM,

  -- ** Matrix multiplication
  (.><),
  (!><),
  multiplyMatrixByVector,
  (><.),
  (><!),
  multiplyVectorByMatrix,
  (.><.),
  (!><!),
  multiplyMatrices,
  multiplyMatricesTransposed,

  -- * Norms
  normL2,

  -- * Simple matrices
  identityMatrix,
  lowerTriangular,
  upperTriangular,
  negateA,
  absA,
  signumA,

  -- * Integral
  quotA,
  remA,
  divA,
  modA,
  quotRemA,
  divModA,

  -- * Fractional
  (./),
  (/.),
  (./.),
  (!/!),
  (.^^),
  recipA,

  -- * Floating
  expA,
  logA,
  sqrtA,
  (.**),
  logBaseA,
  sinA,
  cosA,
  tanA,
  asinA,
  acosA,
  atanA,
  sinhA,
  coshA,
  tanhA,
  asinhA,
  acoshA,
  atanhA,

  -- * RealFrac
  truncateA,
  roundA,
  ceilingA,
  floorA,

  -- * RealFloat
  atan2A,
) where

import Control.Monad (when)
import Control.Scheduler
import qualified Data.Foldable as F
import Data.Function
import Data.Massiv.Array.Delayed.Pull
import Data.Massiv.Array.Delayed.Push
import Data.Massiv.Array.Manifest.Internal
import Data.Massiv.Array.Ops.Construct
import Data.Massiv.Array.Ops.Map as A
import Data.Massiv.Core
import Data.Massiv.Core.Common as A
import Data.Massiv.Core.Operations
import System.IO.Unsafe
import Prelude as P

infixr 8 .^, .^^

infixl 7 !*!, .*., .*, *., !/!, ./., ./, /., `quotA`, `remA`, `divA`, `modA`

infixl 6 !+!, .+., .+, +., !-!, .-., .-, -.

-- | Similar to `liftArray2M`, except it can be applied only to representations
-- with `Numeric` instance and result representation stays the same.
--
-- @since 1.0.0
liftNumArray2M
  :: (Index ix, Numeric r e, MonadThrow m)
  => (e -> e -> e)
  -> Array r ix e
  -> Array r ix e
  -> m (Array r ix e)
liftNumArray2M f a1 a2
  | size a1 == size a2 = pure $ unsafeLiftArray2 f a1 a2
  | isZeroSz sz1 && isZeroSz sz2 = pure $ unsafeResize zeroSz a1
  | otherwise = throwM $ SizeMismatchException sz1 sz2
  where
    !sz1 = size a1
    !sz2 = size a2
{-# INLINE liftNumArray2M #-}

applyExactSize2M
  :: (Index ix, Size r, MonadThrow m)
  => (Array r ix e -> Array r ix e -> Array r ix e)
  -> Array r ix e
  -> Array r ix e
  -> m (Array r ix e)
applyExactSize2M f a1 a2
  | size a1 == size a2 = pure $! f a1 a2
  | isZeroSz sz1 && isZeroSz sz2 = pure $! unsafeResize zeroSz a1
  | otherwise = throwM $! SizeMismatchException sz1 sz2
  where
    !sz1 = size a1
    !sz2 = size a2
{-# INLINE applyExactSize2M #-}

-- | Add two arrays together pointwise. Same as `!+!` but produces monadic computation
-- that allows for handling failure.
--
-- /__Throws Exception__/: `SizeMismatchException` when array sizes do not match.
--
-- @since 0.4.0
(.+.) :: (Index ix, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
(.+.) = applyExactSize2M additionPointwise
{-# INLINE (.+.) #-}

-- | Add two arrays together pointwise. 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.
--
-- ====__Example__
--
-- >>> let a1 = Ix1 0 ... 10
-- >>> let a2 = Ix1 20 ... 30
-- >>> a1 !+! a2
-- Array D Seq (Sz1 11)
--   [ 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40 ]
--
-- @since 0.5.6
(!+!) :: (HasCallStack, Index ix, Numeric r e) => Array r ix e -> Array r ix e -> Array r ix e
(!+!) a1 a2 = throwEither (a1 .+. a2)
{-# INLINE (!+!) #-}

-- | Add a scalar to each element of the array. Array is on the left.
--
-- @since 0.1.0
(.+) :: (Index ix, Numeric r e) => Array r ix e -> e -> Array r ix e
(.+) = plusScalar
{-# INLINE (.+) #-}

-- | Add a scalar to each element of the array. Array is on the right.
--
-- @since 0.4.0
(+.) :: (Index ix, Numeric r e) => e -> Array r ix e -> Array r ix e
(+.) = flip plusScalar
{-# INLINE (+.) #-}

-- | Subtract two arrays pointwise. Same as `!-!` but produces monadic computation that
-- allows for handling failure.
--
-- /__Throws Exception__/: `SizeMismatchException` when array sizes do not match.
--
-- @since 0.4.0
(.-.)
  :: (Index ix, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
(.-.) = applyExactSize2M subtractionPointwise
{-# INLINE (.-.) #-}

-- | Subtract one array from another pointwise. 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.
--
-- ====__Example__
--
-- >>> let a1 = Ix1 0 ... 10
-- >>> let a2 = Ix1 20 ... 30
-- >>> a1 !-! a2
-- Array D Seq (Sz1 11)
--   [ -20, -20, -20, -20, -20, -20, -20, -20, -20, -20, -20 ]
--
-- @since 0.5.6
(!-!) :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e -> Array r ix e
(!-!) a1 a2 = throwEither (a1 .-. a2)
{-# INLINE (!-!) #-}

-- | Subtract a scalar from each element of the array. Array is on the left.
--
-- @since 0.4.0
(.-) :: (Index ix, Numeric r e) => Array r ix e -> e -> Array r ix e
(.-) = minusScalar
{-# INLINE (.-) #-}

-- | Subtract each element of the array from a scalar. Array is on the right.
--
-- @since 0.5.6
(-.) :: (Index ix, Numeric r e) => e -> Array r ix e -> Array r ix e
(-.) = scalarMinus
{-# INLINE (-.) #-}

-- | Multiply two arrays together pointwise. Same as `!*!` but produces monadic
-- computation that allows for handling failure.
--
-- /__Throws Exception__/: `SizeMismatchException` when array sizes do not match.
--
-- @since 0.4.0
(.*.)
  :: (Index ix, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
(.*.) = applyExactSize2M multiplicationPointwise
{-# INLINE (.*.) #-}

-- | 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.
--
-- ====__Example__
--
-- >>> let a1 = Ix1 0 ... 10
-- >>> let a2 = Ix1 20 ... 30
-- >>> a1 !*! a2
-- Array D Seq (Sz1 11)
--   [ 0, 21, 44, 69, 96, 125, 156, 189, 224, 261, 300 ]
--
-- @since 0.5.6
(!*!) :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e -> Array r ix e
(!*!) a1 a2 = throwEither (a1 .*. a2)
{-# INLINE (!*!) #-}

-- | Multiply each element of the array by a scalar value. Scalar is on the right.
--
-- ====__Example__
--
-- >>> let arr = Ix1 20 ..: 25
-- >>> arr
-- Array D Seq (Sz1 5)
--   [ 20, 21, 22, 23, 24 ]
-- >>> arr .* 10
-- Array D Seq (Sz1 5)
--   [ 200, 210, 220, 230, 240 ]
--
-- @since 0.4.0
(.*) :: (Index ix, Numeric r e) => Array r ix e -> e -> Array r ix e
(.*) = multiplyScalar
{-# INLINE (.*) #-}

-- | Multiply each element of the array by a scalar value. Scalar is on the left.
--
-- ====__Example__
--
-- >>> let arr = Ix1 20 ..: 25
-- >>> arr
-- Array D Seq (Sz1 5)
--   [ 20, 21, 22, 23, 24 ]
-- >>> 10 *. arr
-- Array D Seq (Sz1 5)
--   [ 200, 210, 220, 230, 240 ]
--
-- @since 0.4.0
(*.) :: (Index ix, Numeric r e) => e -> Array r ix e -> Array r ix e
(*.) = flip multiplyScalar
{-# INLINE (*.) #-}

-- | Raise each element of the array to a power.
--
-- ====__Example__
--
-- >>> let arr = Ix1 20 ..: 25
-- >>> arr
-- Array D Seq (Sz1 5)
--   [ 20, 21, 22, 23, 24 ]
-- >>> arr .^ 3
-- Array D Seq (Sz1 5)
--   [ 8000, 9261, 10648, 12167, 13824 ]
--
-- @since 0.4.0
(.^) :: (Index ix, Numeric r e) => Array r ix e -> Int -> Array r ix e
(.^) = powerPointwise
{-# INLINE (.^) #-}

-- | Dot product of two vectors.
--
-- [Partial] Throws an impure exception when lengths of vectors do not match
--
-- @since 0.5.6
(!.!) :: (Numeric r e, Source r e) => Vector r e -> Vector r e -> e
(!.!) v1 v2 = throwEither $ dotM v1 v2
{-# INLINE (!.!) #-}

-- | Dot product of two vectors.
--
-- /__Throws Exception__/: `SizeMismatchException` when lengths of vectors do not match
--
-- @since 0.5.6
dotM :: (FoldNumeric r e, Source r e, MonadThrow m) => Vector r e -> Vector r e -> m e
dotM v1 v2
  | size v1 /= size v2 = throwM $ SizeMismatchException (size v1) (size v2)
  | comp == Seq = pure $! unsafeDotProduct v1 v2
  | otherwise = pure $! unsafePerformIO $ unsafeDotProductIO v1 v2
  where
    comp = getComp v1 <> getComp v2
{-# INLINE dotM #-}

unsafeDotProductIO
  :: (MonadUnliftIO m, Index ix, FoldNumeric r b, Source r b)
  => Array r ix b
  -> Array r ix b
  -> m b
unsafeDotProductIO v1 v2 = do
  results <-
    withScheduler comp $ \scheduler ->
      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> liftIO $ do
        let n = SafeSz chunkLength
        loopA_ 0 (< slackStart) (+ chunkLength) $ \ !start ->
          scheduleWork scheduler $
            pure $!
              unsafeDotProduct (unsafeLinearSlice start n v1) (unsafeLinearSlice start n v2)
        when (slackStart < totalLength) $ do
          let k = SafeSz (totalLength - slackStart)
          scheduleWork scheduler $
            pure $!
              unsafeDotProduct (unsafeLinearSlice slackStart k v1) (unsafeLinearSlice slackStart k v2)
  pure $! F.foldl' (+) 0 results
  where
    totalLength = totalElem (size v1)
    comp = getComp v1 <> getComp v2
{-# INLINE unsafeDotProductIO #-}

-- | Compute L2 norm of an array.
--
-- @since 0.5.6
normL2 :: (FoldNumeric r e, Source r e, Index ix, Floating 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, Index ix, FoldNumeric r b, Source r b)
  => Array r ix b
  -> Int
  -> m b
powerSumArrayIO v p = do
  results <-
    withScheduler (getComp v) $ \scheduler ->
      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> liftIO $ do
        let n = SafeSz chunkLength
        loopA_ 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, FoldNumeric r e, Source r 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)
  | mRows == 0 || mCols == 0 = pure $ setComp comp empty
  | otherwise = pure $ makeArray comp (Sz1 mRows) $ \i ->
      unsafeDotProduct (unsafeLinearSlice (i * n) sz mm) v
  where
    comp = getComp mm <> getComp v
    Sz2 mRows mCols = size mm
    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, Manifest r 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
--
-- @since 0.5.6
(!><)
  :: (Numeric r e, Source r e)
  => Matrix r e
  -- ^ Matrix
  -> Vector r e
  -- ^ Column vector (Used many times, so make sure it is computed)
  -> Vector D e
(!><) mm v = throwEither (mm .>< v)
{-# INLINE (!><) #-}

-- | Multiply a row vector by a matrix. 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, Manifest r e)
  => Vector r e
  -- ^ Row vector
  -> Matrix r e
  -- ^ Matrix
  -> m (Vector r e)
(><.) = multiplyVectorByMatrix
{-# 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, Manifest r 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)
  | mRows == 0 || mCols == 0 = pure $ runST (unsafeFreeze comp =<< unsafeNew zeroSz)
  | otherwise =
      pure $!
        unsafePerformIO $ do
          mv <- newMArray (Sz mCols) 0
          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
    Sz1 n = size v
{-# INLINE multiplyVectorByMatrix #-}

-- | Multiply a row vector by a matrix.
--
-- [Partial] Throws impure exception when inner dimensions do not agree
--
-- @since 0.5.6
(><!)
  :: (Numeric r e, Manifest r e)
  => Vector r e
  -- ^ Row vector (Used many times, so make sure it is computed)
  -> Matrix r e
  -- ^ Matrix
  -> Vector r e
(><!) v mm = throwEither (v ><. mm)
{-# INLINE (><!) #-}

-- | Multiply two matrices together.
--
-- [Partial] Inner dimension must agree
--
-- ====__Examples__
--
-- >>> import Data.Massiv.Array
-- >>> a1 = makeArrayR P Seq (Sz2 5 6) $ \(i :. j) -> i + j
-- >>> a2 = makeArrayR P Seq (Sz2 6 5) $ \(i :. j) -> i - j
-- >>> a1 !><! a2
-- Array P Seq (Sz (5 :. 5))
--   [ [ 55, 40, 25, 10, -5 ]
--   , [ 70, 49, 28, 7, -14 ]
--   , [ 85, 58, 31, 4, -23 ]
--   , [ 100, 67, 34, 1, -32 ]
--   , [ 115, 76, 37, -2, -41 ]
--   ]
--
-- @since 0.5.6
(!><!) :: (Numeric r e, Manifest r e) => Matrix r e -> Matrix r e -> Matrix r e
(!><!) a1 a2 = throwEither (a1 `multiplyMatrices` a2)
{-# INLINE (!><!) #-}

-- | Matrix multiplication. 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
(.><.) :: (Numeric r e, Manifest r e, MonadThrow m) => Matrix r e -> Matrix r e -> m (Matrix r e)
(.><.) = multiplyMatrices
{-# INLINE (.><.) #-}

-- | Synonym for `.><.`
--
-- @since 0.5.6
multiplyMatrices
  :: (Numeric r e, Manifest r e, MonadThrow m) => Matrix r e -> Matrix r e -> m (Matrix r e)
multiplyMatrices arrA arrB
  -- mA == 1 = -- TODO: call multiplyVectorByMatrix
  -- nA == 1 = -- TODO: call multiplyMatrixByVector
  | nA /= mB = throwM $ SizeMismatchException (size arrA) (size arrB)
  | isEmpty arrA || isEmpty arrB = pure $ runST (unsafeFreeze comp =<< unsafeNew zeroSz)
  | otherwise = pure $! unsafePerformIO $ do
      marrC <- newMArray (SafeSz (mA :. nB)) 0
      withScheduler_ comp $ \scheduler -> do
        let withC00 iA jB f =
              let !ixC00 = iA * nB + jB
               in f ixC00 =<< unsafeLinearRead marrC ixC00
            withC01 ixC00 f =
              let !ixC01 = ixC00 + 1
               in f ixC01 =<< unsafeLinearRead marrC ixC01
            withC10 ixC00 f =
              let !ixC10 = ixC00 + nB
               in f ixC10 =<< unsafeLinearRead marrC ixC10
            withC11 ixC01 f =
              let !ixC11 = ixC01 + nB
               in f ixC11 =<< unsafeLinearRead marrC ixC11
            withB00 iB jB f =
              let !ixB00 = iB * nB + jB
               in f ixB00 $! unsafeLinearIndex arrB ixB00
            withB00B10 iB jB f =
              withB00 iB jB $ \ixB00 b00 ->
                let !ixB10 = ixB00 + nB
                 in f ixB00 b00 ixB10 $! unsafeLinearIndex arrB ixB10
            withA00 iA jA f =
              let !ixA00 = iA * nA + jA
               in f ixA00 $! unsafeLinearIndex arrA ixA00
            withA00A10 iA jA f =
              withA00 iA jA $ \ixA00 a00 ->
                let !ixA10 = ixA00 + nA
                 in f ixA00 a00 ixA10 $! unsafeLinearIndex arrA ixA10
        let loopColsB_UnRowBColA_UnRowA a00 a01 a10 a11 iA iB jB
              | jB < n2B = do
                  withB00B10 iB jB $ \ixB00 b00 ixB10 b10 -> do
                    let !b01 = unsafeLinearIndex arrB (ixB00 + 1)
                        !b11 = unsafeLinearIndex arrB (ixB10 + 1)
                    withC00 iA jB $ \ixC00 c00 -> do
                      unsafeLinearWrite marrC ixC00 (c00 + a00 * b00 + a01 * b10)
                      withC01 ixC00 $ \ixC01 c01 -> do
                        unsafeLinearWrite marrC ixC01 (c01 + a00 * b01 + a01 * b11)
                        withC10 ixC00 $ \ixC10 c10 ->
                          unsafeLinearWrite marrC ixC10 (c10 + a10 * b00 + a11 * b10)
                        withC11 ixC01 $ \ixC11 c11 ->
                          unsafeLinearWrite marrC ixC11 (c11 + a10 * b01 + a11 * b11)
                  loopColsB_UnRowBColA_UnRowA a00 a01 a10 a11 iA iB (jB + 2)
              | jB < nB = withB00B10 iB jB $ \_ b00 _ b10 ->
                  withC00 iA jB $ \ixC00 c00 -> do
                    unsafeLinearWrite marrC ixC00 (c00 + a00 * b00 + a01 * b10)
                    withC10 ixC00 $ \ixC10 c10 ->
                      unsafeLinearWrite marrC ixC10 (c10 + a10 * b00 + a11 * b10)
              | otherwise = pure ()

            loopColsB_UnRowBColA_RowA a00 a01 iA iB jB
              | jB < n2B = do
                  withB00B10 iB jB $ \ixB00 b00 ixB10 b10 -> do
                    let !b01 = unsafeLinearIndex arrB (ixB00 + 1)
                        !b11 = unsafeLinearIndex arrB (ixB10 + 1)
                    withC00 iA jB $ \ixC00 c00 -> do
                      unsafeLinearWrite marrC ixC00 (c00 + a00 * b00 + a01 * b10)
                      withC01 ixC00 $ \ixC01 c01 ->
                        unsafeLinearWrite marrC ixC01 (c01 + a00 * b01 + a01 * b11)
                  loopColsB_UnRowBColA_RowA a00 a01 iA iB (jB + 2)
              | jB < nB = withB00B10 iB jB $ \_ b00 _ b10 ->
                  withC00 iA jB $ \ixC00 c00 ->
                    unsafeLinearWrite marrC ixC00 (c00 + a00 * b00 + a01 * b10)
              | otherwise = pure ()

            loopColsB_RowBColA_UnRowA a00 a10 iA iB jB
              | jB < n2B = do
                  withB00 iB jB $ \ixB00 b00 -> do
                    let !b01 = unsafeLinearIndex arrB (ixB00 + 1)
                    withC00 iA jB $ \ixC00 c00 -> do
                      unsafeLinearWrite marrC ixC00 (c00 + a00 * b00)
                      withC01 ixC00 $ \ixC01 c01 -> do
                        unsafeLinearWrite marrC ixC01 (c01 + a00 * b01)
                        withC10 ixC00 $ \ixC10 c10 ->
                          unsafeLinearWrite marrC ixC10 (c10 + a10 * b00)
                        withC11 ixC01 $ \ixC11 c11 ->
                          unsafeLinearWrite marrC ixC11 (c11 + a10 * b01)
                  loopColsB_RowBColA_UnRowA a00 a10 iA iB (jB + 2)
              | jB < nB = withB00 iB jB $ \_ b00 ->
                  withC00 iA jB $ \ixC00 c00 -> do
                    unsafeLinearWrite marrC ixC00 (c00 + a00 * b00)
                    withC10 ixC00 $ \ixC10 c10 ->
                      unsafeLinearWrite marrC ixC10 (c10 + a10 * b00)
              | otherwise = pure ()

            loopColsB_RowBColA_RowA a00 iA iB jB
              | jB < n2B = do
                  withB00 iB jB $ \ixB00 b00 -> do
                    let !b01 = unsafeLinearIndex arrB (ixB00 + 1)
                    withC00 iA jB $ \ixC00 c00 -> do
                      unsafeLinearWrite marrC ixC00 (c00 + a00 * b00)
                      withC01 ixC00 $ \ixC01 c01 -> do
                        unsafeLinearWrite marrC ixC01 (c01 + a00 * b01)
                  loopColsB_RowBColA_RowA a00 iA iB (jB + 2)
              | jB < nB = withB00 iB jB $ \_ b00 ->
                  withC00 iA jB $ \ixC00 c00 ->
                    unsafeLinearWrite marrC ixC00 (c00 + a00 * b00)
              | otherwise = pure ()

            loopRowsB_UnRowA iA iB
              | iB < m2B = do
                  withA00A10 iA iB $ \ixA00 a00 ixA10 a10 -> do
                    let !a01 = unsafeLinearIndex arrA (ixA00 + 1)
                        !a11 = unsafeLinearIndex arrA (ixA10 + 1)
                    loopColsB_UnRowBColA_UnRowA a00 a01 a10 a11 iA iB 0
                  loopRowsB_UnRowA iA (iB + 2)
              | iB < mB =
                  withA00A10 iA iB $ \_ a00 _ a10 -> loopColsB_RowBColA_UnRowA a00 a10 iA iB 0
              | otherwise = pure ()

            loopRowsB_RowA iA iB
              | iB < m2B = do
                  withA00 iA iB $ \ixA00 a00 -> do
                    let !a01 = unsafeLinearIndex arrA (ixA00 + 1)
                    loopColsB_UnRowBColA_RowA a00 a01 iA iB 0
                  loopRowsB_RowA iA (iB + 2)
              | iB < mB = withA00 iA iB $ \_ a00 -> loopColsB_RowBColA_RowA a00 iA iB 0
              | otherwise = pure ()

            loopRowsA iA
              | iA < m2A = do
                  scheduleWork_ scheduler $ loopRowsB_UnRowA iA 0
                  loopRowsA (iA + 2)
              | iA < mA = scheduleWork_ scheduler $ loopRowsB_RowA iA 0
              | otherwise = pure ()
        loopRowsA 0

      unsafeFreeze comp marrC
  where
    comp = getComp arrA <> getComp arrB
    m2A = mA - mA `rem` 2
    m2B = mB - mB `rem` 2
    n2B = nB - nB `rem` 2
    Sz (mA :. nA) = size arrA
    Sz (mB :. nB) = size arrB
{-# INLINEABLE multiplyMatrices #-}

-- | Computes the matrix-matrix multiplication where second matrix is transposed (i.e. M
-- x N')
--
-- > m1 .><. transpose m2 == multiplyMatricesTransposed m1 m2
--
-- @since 0.5.6
multiplyMatricesTransposed
  :: (Numeric r e, Manifest r e, MonadThrow m)
  => Matrix r e
  -> Matrix r e
  -> m (Matrix D e)
multiplyMatricesTransposed arr1 arr2
  | n1 /= m2 = throwM $ SizeMismatchException (size arr1) (Sz2 m2 n2)
  | isEmpty arr1 || isEmpty arr2 = pure $ setComp comp empty
  | otherwise =
      pure $
        makeArray comp (SafeSz (m1 :. n2)) $ \(i :. j) ->
          unsafeDotProduct (unsafeLinearSlice (i * n1) n arr1) (unsafeLinearSlice (j * n1) n arr2)
  where
    comp = getComp arr1 <> getComp arr2
    n = SafeSz n1
    SafeSz (m1 :. n1) = size arr1
    SafeSz (n2 :. m2) = size arr2
{-# INLINE multiplyMatricesTransposed #-}

-- | Create an indentity matrix.
--
-- ==== __Example__
--
-- >>> import Data.Massiv.Array
-- >>> identityMatrix 5
-- Array DL Seq (Sz (5 :. 5))
--   [ [ 1, 0, 0, 0, 0 ]
--   , [ 0, 1, 0, 0, 0 ]
--   , [ 0, 0, 1, 0, 0 ]
--   , [ 0, 0, 0, 1, 0 ]
--   , [ 0, 0, 0, 0, 1 ]
--   ]
--
-- @since 0.3.6
identityMatrix :: Num e => Sz1 -> Matrix DL e
identityMatrix (Sz n) =
  makeLoadArrayS (Sz2 n n) 0 $ \w -> loopA_ 0 (< n) (+ 1) $ \i -> w (i :. i) 1
{-# INLINE identityMatrix #-}

-- | Create a lower triangular (L in LU decomposition) matrix of size @NxN@
--
-- ==== __Example__
--
-- >>> import Data.Massiv.Array
-- >>> lowerTriangular Seq 5 (\(i :. j) -> i + j)
-- Array DL Seq (Sz (5 :. 5))
--   [ [ 0, 0, 0, 0, 0 ]
--   , [ 1, 2, 0, 0, 0 ]
--   , [ 2, 3, 4, 0, 0 ]
--   , [ 3, 4, 5, 6, 0 ]
--   , [ 4, 5, 6, 7, 8 ]
--   ]
--
-- @since 0.5.2
lowerTriangular :: forall e. Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
lowerTriangular comp (Sz1 n) f = DLArray comp (SafeSz (n :. n)) load
  where
    load :: Loader e
    load scheduler startAt uWrite uSet = do
      forM_ (0 ..: n) $ \i -> do
        let !k = startAt + i * n
        scheduleWork_ scheduler $ do
          forM_ (0 ... i) $ \j -> uWrite (k + j) (f (i :. j))
          uSet (k + i + 1) (Sz (n - i - 1)) 0
{-# INLINE lowerTriangular #-}

-- | Create an upper triangular (U in LU decomposition) matrix of size @NxN@
--
-- ==== __Example__
--
-- >>> import Data.Massiv.Array
-- >>> upperTriangular Par 5 (\(i :. j) -> i + j)
-- Array DL Par (Sz (5 :. 5))
--   [ [ 0, 1, 2, 3, 4 ]
--   , [ 0, 2, 3, 4, 5 ]
--   , [ 0, 0, 4, 5, 6 ]
--   , [ 0, 0, 0, 6, 7 ]
--   , [ 0, 0, 0, 0, 8 ]
--   ]
--
-- @since 0.5.2
upperTriangular :: forall e. Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
upperTriangular comp (Sz1 n) f = DLArray comp (SafeSz (n :. n)) load
  where
    load :: Loader e
    load scheduler startAt uWrite uSet = do
      forM_ (0 ..: n) $ \i -> do
        let !k = startAt + i * n
        scheduleWork_ scheduler $ do
          uSet k (SafeSz i) 0
          forM_ (i ..: n) $ \j -> uWrite (k + j) (f (i :. j))
{-# INLINE upperTriangular #-}

-- | Negate each element of the array
--
-- @since 0.4.0
negateA :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e
negateA = unsafeLiftArray negate
{-# INLINE negateA #-}

-- | Apply `abs` to each element of the array
--
-- @since 0.4.0
absA :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e
absA = absPointwise
{-# INLINE absA #-}

-- | Apply `signum` to each element of the array
--
-- @since 0.4.0
signumA :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e
signumA = unsafeLiftArray signum
{-# INLINE signumA #-}

-- | Divide each element of one array by another pointwise. Same as `!/!` but produces
-- monadic computation that allows for handling failure.
--
-- /__Throws Exception__/: `SizeMismatchException` when array sizes do not match.
--
-- @since 0.4.0
(./.)
  :: (Index ix, NumericFloat r e, MonadThrow m)
  => Array r ix e
  -> Array r ix e
  -> m (Array r ix e)
(./.) = applyExactSize2M divisionPointwise
{-# INLINE (./.) #-}

-- | Divide two arrays pointwise. 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.
--
-- ====__Example__
--
-- >>> let arr1 = fromIntegral <$> (Ix1 20 ..: 25) :: Array D Ix1 Float
-- >>> let arr2 = fromIntegral <$> (Ix1 100 ..: 105) :: Array D Ix1 Float
-- >>> arr1 !/! arr2
-- Array D Seq (Sz1 5)
--   [ 0.2, 0.20792079, 0.21568628, 0.22330096, 0.23076923 ]
--
-- @since 0.5.6
(!/!) :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e -> Array r ix e
(!/!) a1 a2 = throwEither (a1 ./. a2)
{-# INLINE (!/!) #-}

-- | Divide a scalar value by each element of the array.
--
-- > e /. arr == e *. recipA arr
--
-- ====__Example__
--
-- >>> let arr = fromIntegral <$> (Ix1 20 ..: 25) :: Array D Ix1 Float
-- >>> arr
-- Array D Seq (Sz1 5)
--   [ 20.0, 21.0, 22.0, 23.0, 24.0 ]
-- >>> 100 /. arr
-- Array D Seq (Sz1 5)
--   [ 5.0, 4.7619047, 4.5454545, 4.347826, 4.1666665 ]
--
-- @since 0.5.6
(/.) :: (Index ix, NumericFloat r e) => e -> Array r ix e -> Array r ix e
(/.) = scalarDivide
{-# INLINE (/.) #-}

-- | Divide each element of the array by a scalar value.
--
-- ====__Example__
--
-- >>> let arr = fromIntegral <$> (Ix1 20 ..: 25) :: Array D Ix1 Float
-- >>> arr
-- Array D Seq (Sz1 5)
--   [ 20.0, 21.0, 22.0, 23.0, 24.0 ]
-- >>> arr ./ 100
-- Array D Seq (Sz1 5)
--   [ 0.2, 0.21, 0.22, 0.23, 0.24 ]
--
-- @since 0.4.0
(./) :: (Index ix, NumericFloat r e) => Array r ix e -> e -> Array r ix e
(./) = divideScalar
{-# INLINE (./) #-}

(.^^)
  :: (Index ix, Numeric r e, Fractional e, Integral b)
  => Array r ix e
  -> b
  -> Array r ix e
(.^^) arr n = unsafeLiftArray (^^ n) arr
{-# INLINE (.^^) #-}

-- | Apply reciprocal to each element of the array.
--
-- > recipA arr == 1 /. arr
--
-- @since 0.4.0
recipA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
recipA = recipPointwise
{-# INLINE recipA #-}

-- | Apply exponent to each element of the array.
--
-- > expA arr == map exp arr
--
-- @since 0.4.0
expA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
expA = unsafeLiftArray exp
{-# INLINE expA #-}

-- | Apply square root to each element of the array.
--
-- > sqrtA arr == map sqrt arr
--
-- @since 0.4.0
sqrtA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
sqrtA = unsafeLiftArray sqrt
{-# INLINE sqrtA #-}

-- | Apply logarithm to each element of the array.
--
-- > logA arr == map log arr
--
-- @since 0.4.0
logA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
logA = unsafeLiftArray log
{-# INLINE logA #-}

-- | Apply logarithm to each element of the array where the base is in the same cell in
-- the second array.
--
-- > logBaseA arr1 arr2 == zipWith logBase arr1 arr2
--
-- [Partial] Throws an error when arrays do not have matching sizes
--
-- @since 0.4.0
logBaseA
  :: (Index ix, Source r1 e, Source r2 e, Floating e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
logBaseA = liftArray2' logBase
{-# INLINE logBaseA #-}

-- TODO: siwtch to
-- (breaking) logBaseA :: Array r ix e -> e -> Array D ix e
-- logBasesM :: Array r ix e -> Array r ix e -> m (Array D ix e)

-- | Apply power to each element of the array where the power value is in the same cell
-- in the second array.
--
-- > arr1 .** arr2 == zipWith (**) arr1 arr2
--
-- [Partial] Throws an error when arrays do not have matching sizes
--
-- @since 0.4.0
(.**)
  :: (Index ix, Source r1 e, Source r2 e, Floating e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
(.**) = liftArray2' (**)
{-# INLINE (.**) #-}

-- TODO:
-- !**! :: Array r1 ix e -> Array r2 ix e -> Array D ix e
-- .**. :: Array r1 ix e -> Array r2 ix e -> m (Array D ix e)
-- (breaking) .** :: Array r1 ix e -> e -> Array D ix e

-- | Apply sine function to each element of the array.
--
-- > sinA arr == map sin arr
--
-- @since 0.4.0
sinA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
sinA = unsafeLiftArray sin
{-# INLINE sinA #-}

-- | Apply cosine function to each element of the array.
--
-- > cosA arr == map cos arr
--
-- @since 0.4.0
cosA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
cosA = unsafeLiftArray cos
{-# INLINE cosA #-}

-- | Apply tangent function to each element of the array.
--
-- > tanA arr == map tan arr
--
-- @since 0.4.0
tanA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
tanA = unsafeLiftArray tan
{-# INLINE tanA #-}

-- | Apply arcsine function to each element of the array.
--
-- > asinA arr == map asin arr
--
-- @since 0.4.0
asinA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
asinA = unsafeLiftArray asin
{-# INLINE asinA #-}

-- | Apply arctangent function to each element of the array.
--
-- > atanA arr == map atan arr
--
-- @since 0.4.0
atanA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
atanA = unsafeLiftArray atan
{-# INLINE atanA #-}

-- | Apply arccosine function to each element of the array.
--
-- > acosA arr == map acos arr
--
-- @since 0.4.0
acosA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
acosA = unsafeLiftArray acos
{-# INLINE acosA #-}

-- | Apply hyperbolic sine function to each element of the array.
--
-- > sinhA arr == map sinh arr
--
-- @since 0.4.0
sinhA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
sinhA = unsafeLiftArray sinh
{-# INLINE sinhA #-}

-- | Apply hyperbolic tangent function to each element of the array.
--
-- > tanhA arr == map tanh arr
--
-- @since 0.4.0
tanhA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
tanhA = unsafeLiftArray tanh
{-# INLINE tanhA #-}

-- | Apply hyperbolic cosine function to each element of the array.
--
-- > coshA arr == map cosh arr
--
-- @since 0.4.0
coshA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
coshA = unsafeLiftArray cosh
{-# INLINE coshA #-}

-- | Apply inverse hyperbolic sine function to each element of the array.
--
-- > asinhA arr == map asinh arr
--
-- @since 0.4.0
asinhA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
asinhA = unsafeLiftArray asinh
{-# INLINE asinhA #-}

-- | Apply inverse hyperbolic cosine function to each element of the array.
--
-- > acoshA arr == map acosh arr
--
-- @since 0.4.0
acoshA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
acoshA = unsafeLiftArray acosh
{-# INLINE acoshA #-}

-- | Apply inverse hyperbolic tangent function to each element of the array.
--
-- > atanhA arr == map atanh arr
--
-- @since 0.4.0
atanhA :: (Index ix, NumericFloat r e) => Array r ix e -> Array r ix e
atanhA = unsafeLiftArray atanh
{-# INLINE atanhA #-}

-- | Perform a pointwise quotient where first array contains numerators and the second
-- one denominators
--
-- > quotA arr1 arr2 == zipWith quot arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
quotA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
quotA = liftArray2' quot
{-# INLINE quotA #-}

-- | Perform a pointwise remainder computation
--
-- > remA arr1 arr2 == zipWith rem arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
remA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
remA = liftArray2' rem
{-# INLINE remA #-}

-- | Perform a pointwise integer division where first array contains numerators and the
-- second one denominators
--
-- > divA arr1 arr2 == zipWith div arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
divA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
divA = liftArray2' div
{-# INLINE divA #-}

-- TODO:
--  * Array r ix e -> Array r ix e -> m (Array r ix e)
--  * Array r ix e -> e -> Array r ix e
--  * e -> Array r ix e -> Array r ix e

-- | Perform a pointwise modulo computation
--
-- > modA arr1 arr2 == zipWith mod arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
modA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> Array D ix e
modA = liftArray2' mod
{-# INLINE modA #-}

-- | Perform a pointwise quotient with remainder where first array contains numerators
-- and the second one denominators
--
-- > quotRemA arr1 arr2 == zipWith quotRem arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
quotRemA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> (Array D ix e, Array D ix e)
quotRemA arr1 = A.unzip . liftArray2' quotRem arr1
{-# INLINE quotRemA #-}

-- | Perform a pointwise integer division with modulo where first array contains
-- numerators and the second one denominators
--
-- > divModA arr1 arr2 == zipWith divMod arr1 arr2
--
-- [Partial] Mismatched array sizes will result in an impure exception being thrown.
--
-- @since 0.1.0
divModA
  :: (HasCallStack, Index ix, Source r1 e, Source r2 e, Integral e)
  => Array r1 ix e
  -> Array r2 ix e
  -> (Array D ix e, Array D ix e)
divModA arr1 = A.unzip . liftArray2' divMod arr1
{-# INLINE divModA #-}

-- | Truncate each element of the array.
--
-- > truncateA arr == map truncate arr
--
-- @since 0.1.0
truncateA :: (Index ix, Source r a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
truncateA = A.map truncate
{-# INLINE truncateA #-}

-- | Round each element of the array.
--
-- > truncateA arr == map truncate arr
--
-- @since 0.1.0
roundA :: (Index ix, Source r a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
roundA = A.map round
{-# INLINE roundA #-}

-- | Ceiling of each element of the array.
--
-- > truncateA arr == map truncate arr
--
-- @since 0.1.0
ceilingA :: (Index ix, Source r a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
ceilingA = A.map ceiling
{-# INLINE ceilingA #-}

-- | Floor each element of the array.
--
-- > truncateA arr == map truncate arr
--
-- @since 0.1.0
floorA :: (Index ix, Source r a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
floorA = A.map floor
{-# INLINE floorA #-}

-- | Perform atan2 pointwise
--
-- > atan2A arr1 arr2 == zipWith atan2 arr1 arr2
--
-- /__Throws Exception__/: `SizeMismatchException` when array sizes do not match.
--
-- @since 0.1.0
atan2A
  :: (Index ix, Numeric r e, RealFloat e, MonadThrow m)
  => Array r ix e
  -> Array r ix e
  -> m (Array r ix e)
atan2A = liftNumArray2M atan2
{-# INLINE atan2A #-}

-- | Same as `sumArraysM`, compute sum of arrays pointwise. All arrays must have the same
-- size, otherwise it will result in an error.
--
-- @since 1.0.0
sumArrays'
  :: (HasCallStack, Foldable t, Load r ix e, Numeric r e) => t (Array r ix e) -> Array r ix e
sumArrays' = throwEither . sumArraysM
{-# INLINE sumArrays' #-}

-- | Compute sum of arrays pointwise. All arrays must have the same size.
--
-- ====__Examples__
--
-- >>> import Data.Massiv.Array as A
-- >>> sumArraysM [] :: IO (Array P Ix3 Int)
-- Array P Seq (Sz (0 :> 0 :. 0))
--   [  ]
-- >>> arr = A.makeArrayR P Seq (Sz3 4 5 6) $ \(i :> j :. k) -> i + j * k
-- >>> arr
-- Array P Seq (Sz (4 :> 5 :. 6))
--   [ [ [ 0, 0, 0, 0, 0, 0 ]
--     , [ 0, 1, 2, 3, 4, 5 ]
--     , [ 0, 2, 4, 6, 8, 10 ]
--     , [ 0, 3, 6, 9, 12, 15 ]
--     , [ 0, 4, 8, 12, 16, 20 ]
--     ]
--   , [ [ 1, 1, 1, 1, 1, 1 ]
--     , [ 1, 2, 3, 4, 5, 6 ]
--     , [ 1, 3, 5, 7, 9, 11 ]
--     , [ 1, 4, 7, 10, 13, 16 ]
--     , [ 1, 5, 9, 13, 17, 21 ]
--     ]
--   , [ [ 2, 2, 2, 2, 2, 2 ]
--     , [ 2, 3, 4, 5, 6, 7 ]
--     , [ 2, 4, 6, 8, 10, 12 ]
--     , [ 2, 5, 8, 11, 14, 17 ]
--     , [ 2, 6, 10, 14, 18, 22 ]
--     ]
--   , [ [ 3, 3, 3, 3, 3, 3 ]
--     , [ 3, 4, 5, 6, 7, 8 ]
--     , [ 3, 5, 7, 9, 11, 13 ]
--     , [ 3, 6, 9, 12, 15, 18 ]
--     , [ 3, 7, 11, 15, 19, 23 ]
--     ]
--   ]
-- >>> sumArraysM $ outerSlices arr
-- Array P Seq (Sz (5 :. 6))
--   [ [ 6, 6, 6, 6, 6, 6 ]
--   , [ 6, 10, 14, 18, 22, 26 ]
--   , [ 6, 14, 22, 30, 38, 46 ]
--   , [ 6, 18, 30, 42, 54, 66 ]
--   , [ 6, 22, 38, 54, 70, 86 ]
--   ]
-- >>> sumArraysM $ innerSlices arr
-- Array D Seq (Sz (4 :. 5))
--   [ [ 0, 15, 30, 45, 60 ]
--   , [ 6, 21, 36, 51, 66 ]
--   , [ 12, 27, 42, 57, 72 ]
--   , [ 18, 33, 48, 63, 78 ]
--   ]
--
-- @since 1.0.0
sumArraysM
  :: (Foldable t, Load r ix e, Numeric r e, MonadThrow m) => t (Array r ix e) -> m (Array r ix e)
sumArraysM as =
  case F.toList as of
    [] -> pure empty
    (x : xs) -> F.foldlM (.+.) x xs
{-# INLINE sumArraysM #-}

-- OPTIMIZE: Allocate a single result array and write sums into it incrementally.

-- | Same as `productArraysM`. Compute product of arrays pointwise. All arrays must have
-- the same size, otherwise it
-- will result in an error.
--
-- @since 1.0.0
productArrays'
  :: (HasCallStack, Foldable t, Load r ix e, Numeric r e) => t (Array r ix e) -> Array r ix e
productArrays' = throwEither . productArraysM
{-# INLINE productArrays' #-}

-- | Compute product of arrays pointwise. All arrays must have the same size.
--
-- ====__Examples__
--
-- >>> import Data.Massiv.Array as A
-- >>> productArraysM [] :: IO (Array P Ix3 Int)
-- Array P Seq (Sz (0 :> 0 :. 0))
--   [  ]
-- >>> arr = A.makeArrayR P Seq (Sz3 4 5 6) $ \(i :> j :. k) -> i + j * k
-- >>> arr
-- Array P Seq (Sz (4 :> 5 :. 6))
--   [ [ [ 0, 0, 0, 0, 0, 0 ]
--     , [ 0, 1, 2, 3, 4, 5 ]
--     , [ 0, 2, 4, 6, 8, 10 ]
--     , [ 0, 3, 6, 9, 12, 15 ]
--     , [ 0, 4, 8, 12, 16, 20 ]
--     ]
--   , [ [ 1, 1, 1, 1, 1, 1 ]
--     , [ 1, 2, 3, 4, 5, 6 ]
--     , [ 1, 3, 5, 7, 9, 11 ]
--     , [ 1, 4, 7, 10, 13, 16 ]
--     , [ 1, 5, 9, 13, 17, 21 ]
--     ]
--   , [ [ 2, 2, 2, 2, 2, 2 ]
--     , [ 2, 3, 4, 5, 6, 7 ]
--     , [ 2, 4, 6, 8, 10, 12 ]
--     , [ 2, 5, 8, 11, 14, 17 ]
--     , [ 2, 6, 10, 14, 18, 22 ]
--     ]
--   , [ [ 3, 3, 3, 3, 3, 3 ]
--     , [ 3, 4, 5, 6, 7, 8 ]
--     , [ 3, 5, 7, 9, 11, 13 ]
--     , [ 3, 6, 9, 12, 15, 18 ]
--     , [ 3, 7, 11, 15, 19, 23 ]
--     ]
--   ]
-- >>> productArraysM $ outerSlices arr
-- Array P Seq (Sz (5 :. 6))
--   [ [ 0, 0, 0, 0, 0, 0 ]
--   , [ 0, 24, 120, 360, 840, 1680 ]
--   , [ 0, 120, 840, 3024, 7920, 17160 ]
--   , [ 0, 360, 3024, 11880, 32760, 73440 ]
--   , [ 0, 840, 7920, 32760, 93024, 212520 ]
--   ]
-- >>> productArraysM $ innerSlices arr
-- Array D Seq (Sz (4 :. 5))
--   [ [ 0, 0, 0, 0, 0 ]
--   , [ 1, 720, 10395, 58240, 208845 ]
--   , [ 64, 5040, 46080, 209440, 665280 ]
--   , [ 729, 20160, 135135, 524880, 1514205 ]
--   ]
--
-- @since 1.0.0
productArraysM
  :: (Foldable t, Load r ix e, Numeric r e, MonadThrow m) => t (Array r ix e) -> m (Array r ix e)
productArraysM as =
  case F.toList as of
    [] -> pure empty
    (x : xs) -> F.foldlM (.*.) x xs
{-# INLINE productArraysM #-}