diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,17 @@
+# 0.5.6
+
+* Fix `(-.)` (it was incorrectly implemented as a flip of `(.-)`
+* Addition of numeric functions:
+  * Partial: `!+!`, `!-!`, `!*!`, `!/!`
+  * Reciprocal division `/.`
+  * More efficient matrix-matrix multiplication: `.><.` and `!><!` (also helpers
+    `multiplyMatrices` and `multiplyMatricesTransposed`)
+  * More efficient matrix-vector multiplication: `.><` and `!><`
+  * New vector-matrix multiplication: `><.` and `><!`
+  * Dot product `dotM` and `!.!`
+  * Norm `normL2`
+* Deprecated `|*|` and `#>`
+
 # 0.5.5
 
 * Add `takeWhile`, `dropWhile` and `findIndex`
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.5.5.0
+version:             0.5.6.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
@@ -18,7 +18,8 @@
                     , GHC == 8.6.3
                     , GHC == 8.6.4
                     , GHC == 8.6.5
-                    , GHC == 8.8.1
+                    , GHC == 8.8.4
+                    , GHC == 8.10.2
 
 flag unsafe-checks
   description: Enable all the bounds checks for unsafe functions at the cost of
diff --git a/src/Data/Massiv/Array/Delayed/Pull.hs b/src/Data/Massiv/Array/Delayed/Pull.hs
--- a/src/Data/Massiv/Array/Delayed/Pull.hs
+++ b/src/Data/Massiv/Array/Delayed/Pull.hs
@@ -3,7 +3,6 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE NamedFieldPuns #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UndecidableInstances #-}
@@ -219,26 +218,13 @@
 
 
 instance Num e => Numeric D e where
-  -- plusScalar arr e = unsafeLiftArray (+ e) arr
-  -- {-# INLINE plusScalar #-}
-  -- minusScalar arr e = unsafeLiftArray (subtract e) arr
-  -- {-# INLINE minusScalar #-}
-  -- multiplyScalar arr e = unsafeLiftArray (* e) arr
-  -- {-# INLINE multiplyScalar #-}
-  -- absPointwise = unsafeLiftArray abs
-  -- {-# INLINE absPointwise #-}
-  -- additionPointwise = unsafeLiftArray2 (+)
-  -- {-# INLINE additionPointwise #-}
-  -- subtractionPointwise = unsafeLiftArray2 (-)
-  -- {-# INLINE subtractionPointwise #-}
-  -- multiplicationPointwise = unsafeLiftArray2 (*)
-  -- {-# INLINE multiplicationPointwise #-}
-  -- powerPointwise arr pow = unsafeLiftArray (^ pow) arr
-  -- {-# INLINE powerPointwise #-}
-  -- powerSumArray arr = sumArray . powerPointwise arr
-  -- {-# INLINE powerSumArray #-}
-  -- unsafeDotProduct a1 a2 = sumArray (multiplicationPointwise a1 a2)
-  -- {-# INLINE unsafeDotProduct #-}
+  foldArray f !initAcc arr = go initAcc 0
+    where
+      !len = totalElem (dSize arr)
+      go !acc i
+        | i < len = go (f acc (unsafeLinearIndex arr i)) (i + 1)
+        | otherwise = acc
+  {-# INLINE foldArray #-}
   unsafeLiftArray f arr = arr {dIndex = f . dIndex arr}
   {-# INLINE unsafeLiftArray #-}
   unsafeLiftArray2 f a1 a2 =
@@ -247,19 +233,7 @@
   {-# INLINE unsafeLiftArray2 #-}
 
 
-instance Floating e => NumericFloat D e where
-  -- recipPointwise = liftDArray recip
-  -- {-# INLINE recipPointwise #-}
-  -- sqrtPointwise = liftDArray sqrt
-  -- {-# INLINE sqrtPointwise #-}
-  -- floorPointwise = liftDArray floor
-  -- {-# INLINE floorPointwise #-}
-  -- ceilingPointwise = liftDArray ceiling
-  -- {-# INLINE ceilingPointwise #-}
-  -- divisionPointwise = liftDArray2 (/)
-  -- {-# INLINE divisionPointwise #-}
-  -- divideScalar arr e = liftDArray (/ e) arr
-  -- {-# INLINE divideScalar #-}
+instance Floating e => NumericFloat D e
 
 
 
@@ -272,15 +246,14 @@
 "delay" [~1] forall (arr :: Array D ix e) . delay arr = arr
  #-}
 
--- TODO: switch to zipWith
 -- | /O(min (n1, n2))/ - Compute array equality by applying a comparing function to each element.
 eq :: (Source r1 ix e1, Source r2 ix e2) =>
       (e1 -> e2 -> Bool) -> Array r1 ix e1 -> Array r2 ix e2 -> Bool
 eq f arr1 arr2 =
   (size arr1 == size arr2) &&
-  F.and
-    (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
-       f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
+  not (A.any not
+       (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
+           f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix)))
 {-# INLINE eq #-}
 
 -- | /O(min (n1, n2))/ - Compute array ordering by applying a comparing function to each element.
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
@@ -50,11 +50,12 @@
 
 import Control.DeepSeq (NFData(..), deepseq)
 import Control.Monad.Primitive (PrimMonad(..), primitive_)
-import Data.Massiv.Array.Delayed.Pull (eq, ord)
+import Data.Massiv.Array.Delayed.Pull -- (eq, ord)
 import Data.Massiv.Array.Manifest.Internal
 import Data.Massiv.Array.Manifest.List as A
 import Data.Massiv.Array.Mutable
 import Data.Massiv.Core.Common
+import Data.Massiv.Core.Operations
 import Data.Massiv.Core.List
 import Data.Massiv.Vector.Stream as S (steps, isteps)
 import Data.Maybe (fromMaybe)
@@ -97,8 +98,8 @@
   setComp c arr = arr { pComp = c }
   {-# INLINE setComp #-}
 
-  makeArray !comp !sz f = unsafePerformIO $ generateArray comp sz (return . f)
-  {-# INLINE makeArray #-}
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (return . f)
+  {-# INLINE makeArrayLinear #-}
 
 instance (Prim e, Index ix) => Source P ix e where
   unsafeLinearIndex _arr@(PArray _ _ o a) i =
@@ -244,6 +245,35 @@
   {-# INLINE toStream #-}
   toStreamIx = S.isteps
   {-# INLINE toStreamIx #-}
+
+
+instance (Prim e, Num e) => Numeric P 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 (size arr)
+      go !acc i
+        | i < len = go (f acc (unsafeLinearIndex arr i)) (i + 1)
+        | otherwise = acc
+  {-# INLINE foldArray #-}
+  unsafeLiftArray f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 f a1 a2 =
+    makeArrayLinear
+      (getComp a1 <> getComp a2)
+      (SafeSz (liftIndex2 min (unSz (size a1)) (unSz (size a2)))) $ \ !i ->
+      f (unsafeLinearIndex a1 i) (unsafeLinearIndex a2 i)
+  {-# INLINE unsafeLiftArray2 #-}
+
+
+instance (Prim e, Floating e) => NumericFloat P e
+
 
 instance ( Prim e
          , IsList (Array L ix e)
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
@@ -18,7 +18,7 @@
 module Data.Massiv.Array.Manifest.Storable
   ( S (..)
   , Array(..)
-  , VS.Storable
+  , Storable
   , toStorableVector
   , toStorableMVector
   , fromStorableVector
@@ -45,6 +45,7 @@
 import Data.Massiv.Array.Mutable
 import Data.Massiv.Core.Common
 import Data.Massiv.Core.List
+import Data.Massiv.Core.Operations
 import qualified Data.Vector.Generic.Mutable as VGM
 import qualified Data.Vector.Storable as VS
 import qualified Data.Vector.Storable.Mutable as MVS
@@ -67,7 +68,7 @@
                                     , sData :: !(VS.Vector e)
                                     }
 
-instance (Ragged L ix e, Show e, VS.Storable e) => Show (Array S ix e) where
+instance (Ragged L ix e, Show e, Storable e) => Show (Array S ix e) where
   showsPrec = showsArrayPrec id
   showList = showArrayList
 
@@ -75,15 +76,15 @@
   rnf (SArray c sz v) = c `deepseq` sz `deepseq` v `deepseq` ()
   {-# INLINE rnf #-}
 
-instance (VS.Storable e, Eq e, Index ix) => Eq (Array S ix e) where
+instance (Storable e, Eq e, Index ix) => Eq (Array S ix e) where
   (==) = eq (==)
   {-# INLINE (==) #-}
 
-instance (VS.Storable e, Ord e, Index ix) => Ord (Array S ix e) where
+instance (Storable e, Ord e, Index ix) => Ord (Array S ix e) where
   compare = ord compare
   {-# INLINE compare #-}
 
-instance (VS.Storable e, Index ix) => Construct S ix e where
+instance (Storable e, Index ix) => Construct S ix e where
   setComp c arr = arr { sComp = c }
   {-# INLINE setComp #-}
 
@@ -91,7 +92,7 @@
   {-# INLINE makeArray #-}
 
 
-instance (VS.Storable e, Index ix) => Source S ix e where
+instance (Storable e, Index ix) => Source S ix e where
   unsafeLinearIndex (SArray _ _ v) =
     INDEX_CHECK("(Source S ix e).unsafeLinearIndex", Sz . VS.length, VS.unsafeIndex) v
   {-# INLINE unsafeLinearIndex #-}
@@ -102,13 +103,13 @@
   unsafeResize !sz !arr = arr { sSize = sz }
   {-# INLINE unsafeResize #-}
 
-instance (VS.Storable e, Index ix) => Extract S ix e where
+instance (Storable e, Index ix) => Extract S ix e where
   unsafeExtract !sIx !newSz !arr = unsafeExtract sIx newSz (toManifest arr)
   {-# INLINE unsafeExtract #-}
 
 
 
-instance ( VS.Storable e
+instance ( Storable e
          , Index ix
          , Index (Lower ix)
          , Elt M ix e ~ Array M (Lower ix) e
@@ -118,7 +119,7 @@
   unsafeOuterSlice arr = unsafeOuterSlice (toManifest arr)
   {-# INLINE unsafeOuterSlice #-}
 
-instance ( VS.Storable e
+instance ( Storable e
          , Index ix
          , Index (Lower ix)
          , Elt M ix e ~ Array M (Lower ix) e
@@ -128,19 +129,19 @@
   unsafeInnerSlice arr = unsafeInnerSlice (toManifest arr)
   {-# INLINE unsafeInnerSlice #-}
 
-instance {-# OVERLAPPING #-} VS.Storable e => Slice S Ix1 e where
+instance {-# OVERLAPPING #-} Storable e => Slice S Ix1 e where
   unsafeSlice arr i _ _ = pure (unsafeLinearIndex arr i)
   {-# INLINE unsafeSlice #-}
 
 
-instance (Index ix, VS.Storable e) => Manifest S ix e where
+instance (Index ix, Storable e) => Manifest S ix e where
 
   unsafeLinearIndexM (SArray _ _ v) =
     INDEX_CHECK("(Manifest S ix e).unsafeLinearIndexM", Sz . VS.length, VS.unsafeIndex) v
   {-# INLINE unsafeLinearIndexM #-}
 
 
-instance (Index ix, VS.Storable e) => Mutable S ix e where
+instance (Index ix, Storable e) => Mutable S ix e where
   data MArray s S ix e = MSArray !(Sz ix) !(VS.MVector s e)
 
   msize (MSArray sz _) = sz
@@ -204,7 +205,7 @@
   {-# INLINE unsafeLinearGrow #-}
 
 
-instance (Index ix, VS.Storable e) => Load S ix e where
+instance (Index ix, Storable e) => Load S ix e where
   type R S = M
   size = sSize
   {-# INLINE size #-}
@@ -213,16 +214,44 @@
   loadArrayM !scheduler !arr = splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
   {-# INLINE loadArrayM #-}
 
-instance (Index ix, VS.Storable e) => StrideLoad S ix e
+instance (Index ix, Storable e) => StrideLoad S ix e
 
-instance (Index ix, VS.Storable e) => Stream S ix e where
+instance (Index ix, Storable e) => Stream S ix e where
   toStream = S.steps
   {-# INLINE toStream #-}
   toStreamIx = S.isteps
   {-# INLINE toStreamIx #-}
 
 
-instance ( VS.Storable e
+instance (Storable e, Num e) => Numeric S 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 (size arr)
+      go !acc i
+        | i < len = go (f acc (unsafeLinearIndex arr i)) (i + 1)
+        | otherwise = acc
+  {-# INLINE foldArray #-}
+  unsafeLiftArray f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 f a1 a2 =
+    makeArrayLinear
+      (getComp a1 <> getComp a2)
+      (SafeSz (liftIndex2 min (unSz (size a1)) (unSz (size a2)))) $ \ !i ->
+      f (unsafeLinearIndex a1 i) (unsafeLinearIndex a2 i)
+  {-# INLINE unsafeLiftArray2 #-}
+
+
+instance (Storable e, Floating e) => NumericFloat S e
+
+
+instance ( Storable e
          , IsList (Array L ix e)
          , Nested LN ix e
          , Nested L ix e
@@ -239,7 +268,7 @@
 -- referential transparency.
 --
 -- @since 0.1.3
-unsafeWithPtr :: (MonadUnliftIO m, VS.Storable a) => Array S ix a -> (Ptr a -> m b) -> m b
+unsafeWithPtr :: (MonadUnliftIO m, Storable a) => Array S ix a -> (Ptr a -> m b) -> m b
 unsafeWithPtr arr f = withRunInIO $ \run -> VS.unsafeWith (sData arr) (run . f)
 {-# INLINE unsafeWithPtr #-}
 
@@ -247,7 +276,7 @@
 -- | A pointer to the beginning of the mutable array.
 --
 -- @since 0.1.3
-withPtr :: (MonadUnliftIO m, VS.Storable a) => MArray RealWorld S ix a -> (Ptr a -> m b) -> m b
+withPtr :: (MonadUnliftIO m, Storable a) => MArray RealWorld S ix a -> (Ptr a -> m b) -> m b
 withPtr (MSArray _ mv) f = withRunInIO $ \run -> MVS.unsafeWith mv (run . f)
 {-# INLINE withPtr #-}
 
@@ -285,21 +314,21 @@
 -- | /O(1)/ - Yield the underlying `ForeignPtr` together with its length.
 --
 -- @since 0.3.0
-unsafeArrayToForeignPtr :: VS.Storable e => Array S ix e -> (ForeignPtr e, Int)
+unsafeArrayToForeignPtr :: Storable e => Array S ix e -> (ForeignPtr e, Int)
 unsafeArrayToForeignPtr = VS.unsafeToForeignPtr0 . toStorableVector
 {-# INLINE unsafeArrayToForeignPtr #-}
 
 -- | /O(1)/ - Yield the underlying `ForeignPtr` together with its length.
 --
 -- @since 0.3.0
-unsafeMArrayToForeignPtr :: VS.Storable e => MArray s S ix e -> (ForeignPtr e, Int)
+unsafeMArrayToForeignPtr :: Storable e => MArray s S ix e -> (ForeignPtr e, Int)
 unsafeMArrayToForeignPtr = MVS.unsafeToForeignPtr0 . toStorableMVector
 {-# INLINE unsafeMArrayToForeignPtr #-}
 
 -- | /O(1)/ - Wrap a `ForeignPtr` and it's size into a pure storable array.
 --
 -- @since 0.3.0
-unsafeArrayFromForeignPtr0 :: VS.Storable e => Comp -> ForeignPtr e -> Sz1 -> Array S Ix1 e
+unsafeArrayFromForeignPtr0 :: Storable e => Comp -> ForeignPtr e -> Sz1 -> Array S Ix1 e
 unsafeArrayFromForeignPtr0 comp ptr sz =
   SArray {sComp = comp, sSize = sz, sData = VS.unsafeFromForeignPtr0 ptr (unSz sz)}
 {-# INLINE unsafeArrayFromForeignPtr0 #-}
@@ -307,7 +336,7 @@
 -- | /O(1)/ - Wrap a `ForeignPtr`, an offset and it's size into a pure storable array.
 --
 -- @since 0.3.0
-unsafeArrayFromForeignPtr :: VS.Storable e => Comp -> ForeignPtr e -> Int -> Sz1 -> Array S Ix1 e
+unsafeArrayFromForeignPtr :: Storable e => Comp -> ForeignPtr e -> Int -> Sz1 -> Array S Ix1 e
 unsafeArrayFromForeignPtr comp ptr offset sz =
   SArray {sComp = comp, sSize = sz, sData = VS.unsafeFromForeignPtr ptr offset (unSz sz)}
 {-# INLINE unsafeArrayFromForeignPtr #-}
@@ -317,7 +346,7 @@
 -- modify the pointer, unless the array gets frozen prior to modification.
 --
 -- @since 0.3.0
-unsafeMArrayFromForeignPtr0 :: VS.Storable e => ForeignPtr e -> Sz1 -> MArray s S Ix1 e
+unsafeMArrayFromForeignPtr0 :: Storable e => ForeignPtr e -> Sz1 -> MArray s S Ix1 e
 unsafeMArrayFromForeignPtr0 fp sz =
   MSArray sz (MVS.unsafeFromForeignPtr0 fp (unSz sz))
 {-# INLINE unsafeMArrayFromForeignPtr0 #-}
@@ -327,7 +356,7 @@
 -- still safe to modify the pointer, unless the array gets frozen prior to modification.
 --
 -- @since 0.3.0
-unsafeMArrayFromForeignPtr :: VS.Storable e => ForeignPtr e -> Int -> Sz1 -> MArray s S Ix1 e
+unsafeMArrayFromForeignPtr :: Storable e => ForeignPtr e -> Int -> Sz1 -> MArray s S Ix1 e
 unsafeMArrayFromForeignPtr fp offset sz =
   MSArray sz (MVS.unsafeFromForeignPtr fp offset (unSz sz))
 {-# INLINE unsafeMArrayFromForeignPtr #-}
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
@@ -12,20 +12,44 @@
 -- Portability : non-portable
 --
 module Data.Massiv.Array.Numeric
-  ( -- * Num
-    (.+.)
+  ( -- * Numeric
+    Numeric
+  , NumericFloat
+    -- ** Pointwise addition
   , (.+)
   , (+.)
-  , (.-.)
+  , (.+.)
+  , (!+!)
+  -- ** Pointwise subtraction
   , (.-)
   , (-.)
-  , (.*.)
+  , (.-.)
+  , (!-!)
+  -- ** Pointwise multiplication
   , (.*)
   , (*.)
+  , (.*.)
+  , (!*!)
   , (.^)
+  -- ** Dot product
+  , (!.!)
+  , dotM
+  -- ** Matrix multiplication
+  , (.><)
+  , (!><)
+  , (><.)
+  , (><!)
+  , (.><.)
+  , (!><!)
+  , multiplyMatrices
+  , multiplyMatricesTransposed
+  -- Deprecated:
   , (#>)
   , (|*|)
   , multiplyTransposed
+  -- * Norms
+  , normL2
+  -- ** Simple matrices
   , identityMatrix
   , lowerTriangular
   , upperTriangular
@@ -41,8 +65,10 @@
   , quotRemA
   , divModA
   -- * Fractional
-  , (./.)
   , (./)
+  , (/.)
+  , (./.)
+  , (!/!)
   , (.^^)
   , recipA
   , fromRationalA
@@ -74,6 +100,7 @@
   , atan2A
   ) where
 
+import Data.Massiv.Array.Manifest
 import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Array.Delayed.Push
 import Data.Massiv.Array.Manifest.Internal
@@ -85,11 +112,14 @@
 import Data.Massiv.Core.Common
 import Data.Massiv.Core.Operations
 import Prelude as P
-
+import System.IO.Unsafe
+import Control.Scheduler
+import Control.Monad (when)
+import qualified Data.Foldable as F
 
 infixr 8  .^, .^^
-infixl 7  .*., .*, *., ./., ./, `quotA`, `remA`, `divA`, `modA`
-infixl 6  .+., .+, +., .-., .-, -.
+infixl 7  !*!, .*., .*, *., !/!, ./., ./, /., `quotA`, `remA`, `divA`, `modA`
+infixl 6  !+!, .+., .+, +., !-!, .-., .-, -.
 
 liftArray2Matching
   :: (Source r1 ix a, Source r2 ix b)
@@ -130,15 +160,35 @@
 {-# INLINE liftNumericArray2M #-}
 
 
--- | Add two arrays together pointwise. Throws `SizeMismatchException` if arrays sizes do
--- not match.
+-- | 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
 (.+.) ::
      (Load r ix e, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
 (.+.) = liftNumericArray2M 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
+(!+!) :: (Load r ix e, 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
@@ -153,9 +203,11 @@
 (+.) = flip plusScalar
 {-# INLINE (+.) #-}
 
--- | Subtract two arrays pointwise. Throws `SizeMismatchException` if arrays sizes do not
--- match.
+-- | 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
 (.-.) ::
      (Load r ix e, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
@@ -163,57 +215,123 @@
 {-# 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
+(!-!) :: (Load r ix e, 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.1.0
+-- @since 0.4.0
 (.-) :: (Index ix, Numeric r e) => Array r ix e -> e -> Array r ix e
 (.-) = minusScalar
 {-# INLINE (.-) #-}
 
--- | Subtract a scalar from each element of the array. Array is on the right.
+-- | Subtract each element of the array from a scalar. Array is on the right.
 --
--- @since 0.4.0
+-- @since 0.5.6
 (-.) :: (Index ix, Numeric r e) => e -> Array r ix e -> Array r ix e
-(-.) = flip minusScalar
+(-.) = scalarMinus
 {-# INLINE (-.) #-}
 
 
--- | Multiply two arrays together pointwise.
+-- | 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
 (.*.) ::
      (Load r ix e, Numeric r e, MonadThrow m) => Array r ix e -> Array r ix e -> m (Array r ix e)
 (.*.) = liftNumericArray2M multiplicationPointwise
 {-# INLINE (.*.) #-}
 
+
+-- | Multiplication of 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 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
+(!*!) :: (Load r ix e, 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 (.^) #-}
 
--- | Matrix multiplication
---
--- Inner dimensions must agree, otherwise `SizeMismatchException`.
-(|*|) ::
-     (Mutable r Ix2 e, Source r' Ix2 e, OuterSlice r Ix2 e, Source (R r) Ix1 e, Num e, MonadThrow m)
-  => Array r Ix2 e
-  -> Array r' Ix2 e
-  -> m (Array r Ix2 e)
-(|*|) a1 a2 = compute <$> multArrs a1 a2
-{-# INLINE [1] (|*|) #-}
-
-{-# RULES
-"multDoubleTranspose" [~1] forall arr1 arr2 . arr1 |*| transpose arr2 =
-    multiplyTransposedFused arr1 (convert arr2)
- #-}
-
 -- | Matrix-vector product
 --
 -- Inner dimensions must agree, otherwise `SizeMismatchException`
@@ -231,9 +349,244 @@
     Sz2 mRows mCols = size mm
     Sz1 n = size v
 {-# INLINE (#>) #-}
+{-# DEPRECATED (#>) "In favor of (`.><`)" #-}
 
 
+-- | 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 Ix1 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 :: (Numeric r e, Source r Ix1 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 #-}
+
+-- | 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)
+  => Array r ix b
+  -> Array r ix b
+  -> m b
+unsafeDotProductIO v1 v2 = do
+  results <-
+    withScheduler comp $ \scheduler ->
+      splitLinearly (numWorkers scheduler) totalLength $ \chunkLength slackStart -> do
+        let n = SafeSz chunkLength
+        loopM_ 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 #-}
+
+
+-- | 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)
+(.><) mm v
+  | mCols /= n = throwM $ SizeMismatchException (size mm) (Sz2 n 1)
+  | otherwise = pure $ makeArray (getComp mm <> getComp v) (Sz1 mRows) $ \i ->
+      unsafeDotProduct (unsafeLinearSlice (i * n) sz mm) v
+  where
+    Sz2 mRows mCols = size mm
+    sz@(Sz1 n) = size v
+{-# INLINE (.><) #-}
+
+
+-- | 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 Ix1 e, Source r Ix2 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, Mutable r Ix1 e, Mutable r Ix2 e) =>
+         Vector r e -- ^ Row vector
+      -> Matrix r e -- ^ Matrix
+      -> m (Vector D e)
+(><.) 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
+  where
+    Sz2 mRows mCols = size mm
+    mm' = compute (transpose mm)
+    sz@(Sz1 n) = size v
+{-# INLINE (><.) #-}
+
+
+-- | Multiply a row vector by a matrix.
+--
+-- [Partial] Throws impure exception when inner dimensions do not agree
+--
+-- @since 0.5.6
+(><!) ::
+     (Numeric r e, Mutable r Ix1 e, Mutable r Ix2 e)
+  => Vector r e -- ^ Row vector
+  -> Matrix r e -- ^ Matrix
+  -> Vector D e
+(><!) v mm = throwEither (v ><. mm)
+{-# INLINE (><!) #-}
+
+
+
+-- | Multiply two matrices together.
+--
+-- [Partial] Inner dimension must agree
+--
+-- ====__Examples__
+--
+-- >>> a1 = makeArrayR P Seq (Sz2 5 6) $ \(i :. j) -> i + j
+-- >>> a2 = makeArrayR D Seq (Sz2 6 5) $ \(i :. j) -> i - j
+-- >>> a1 !><! a2
+-- Array D 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, Mutable r Ix2 e, Source r' Ix2 e) => Matrix r e -> Matrix r' e -> Matrix D 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, Mutable r Ix2 e, Source r' Ix2 e, MonadThrow m)
+  => Matrix r e
+  -> Matrix r' e
+  -> m (Matrix D e)
+(.><.) = multiplyMatrices
+{-# INLINE (.><.) #-}
+
+
+-- | Synonym for `.><.`
+--
+-- @since 0.5.6
+multiplyMatrices ::
+     forall r r' e m. (Numeric r e, Mutable r Ix2 e, Source r' Ix2 e, MonadThrow m)
+  => Matrix r e
+  -> Matrix r' e
+  -> m (Matrix D e)
+multiplyMatrices arr1 arr2 = multiplyMatricesTransposed arr1 arr2'
+  where
+    arr2' :: Array r Ix2 e
+    arr2' = compute $ transpose arr2
+{-# INLINE [1] multiplyMatrices #-}
+
+{-# RULES
+"multiplyMatricesTransposed" [~1] forall arr1 arr2 . multiplyMatrices arr1 (transpose arr2) =
+    multiplyMatricesTransposedFused arr1 (convert arr2)
+ #-}
+
+-- | 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 Ix2 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)
+  | otherwise =
+    pure $
+    DArray (getComp arr1 <> getComp arr2) (SafeSz (m1 :. n2)) $ \(i :. j) ->
+      unsafeDotProduct (unsafeLinearSlice (i * n1) n arr1) (unsafeLinearSlice (j * n1) n arr2)
+  where
+    n = SafeSz n1
+    SafeSz (m1 :. n1) = size arr1
+    SafeSz (n2 :. m2) = size arr2
+{-# INLINE multiplyMatricesTransposed #-}
+
+multiplyMatricesTransposedFused ::
+     (Manifest r Ix2 e, Numeric r e, MonadThrow m) => Matrix r e -> Matrix r e -> m (Matrix D e)
+multiplyMatricesTransposedFused arr1 arr2 = multiplyMatricesTransposed arr1 arr2
+{-# INLINE multiplyMatricesTransposedFused #-}
+
+
+
+--- Below is outdated
+
+-- | Matrix multiplication
+--
+-- /__Throws Exception__/: `SizeMismatchException` when inner dimensions of arrays do not match.
+--
+-- @since 0.4.0
+(|*|) ::
+     (Mutable r Ix2 e, Source r' Ix2 e, OuterSlice r Ix2 e, Source (R r) Ix1 e, Num e, MonadThrow m)
+  => Array r Ix2 e
+  -> Array r' Ix2 e
+  -> m (Array r Ix2 e)
+(|*|) a1 a2 = compute <$> multArrs a1 a2
+{-# INLINE [1] (|*|) #-}
+{-# DEPRECATED (|*|) "In favor of `.><.`" #-}
+
+{-# RULES
+"multDoubleTranspose" [~1] forall arr1 arr2 . arr1 |*| transpose arr2 =
+    multiplyTransposedFused arr1 (convert arr2)
+ #-}
+
+
 multiplyTransposedFused ::
      ( Mutable r Ix2 e
      , OuterSlice r Ix2 e
@@ -357,23 +710,39 @@
              in wr (toLinearIndex sz ix) (f ix)
 {-# 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 #-}
 
+-- | Create a singleton array from an `Integer`
 fromIntegerA :: (Index ix, Num e) => Integer -> Array D ix e
 fromIntegerA = singleton . fromInteger
 {-# INLINE fromIntegerA #-}
+{-# DEPRECATED fromIntegerA "This almost never a desired behavior. Use `Data.Massiv.Array.Ops.Construct.replicate` instead" #-}
 
+-- | 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
 (./.) ::
      (Load r ix e, NumericFloat r e, MonadThrow m)
   => Array r ix e
@@ -382,6 +751,57 @@
 (./.) = liftNumericArray2M 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
+(!/!) :: (Load r ix e, 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 (./) #-}
@@ -392,6 +812,11 @@
 (.^^) 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 #-}
@@ -402,88 +827,201 @@
   => Rational -> Array D ix e
 fromRationalA = singleton . fromRational
 {-# INLINE fromRationalA #-}
+{-# DEPRECATED fromRationalA "This almost never a desired behavior. Use `Data.Massiv.Array.Ops.Construct.replicate` instead" #-}
 
 piA
   :: (Index ix, Floating e)
   => Array D ix e
 piA = singleton pi
 {-# INLINE piA #-}
+{-# DEPRECATED piA "This almost never a desired behavior. Use `Data.Massiv.Array.Ops.Construct.replicate` instead" #-}
 
+
+-- | 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
   :: (Source r1 ix e, Source r2 ix e, Floating e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
 logBaseA = liftArray2Matching 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
 (.**)
   :: (Source r1 ix e, Source r2 ix e, Floating e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
 (.**) = liftArray2Matching (**)
 {-# 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
@@ -491,18 +1029,44 @@
 {-# 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
 remA = liftArray2Matching 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
 divA = liftArray2Matching 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> Array D ix e
@@ -511,6 +1075,14 @@
 
 
 
+-- | 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> (Array D ix e, Array D ix e)
@@ -518,6 +1090,14 @@
 {-# 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
   :: (Source r1 ix e, Source r2 ix e, Integral e)
   => Array r1 ix e -> Array r2 ix e -> (Array D ix e, Array D ix e)
@@ -526,27 +1106,54 @@
 
 
 
+-- | Truncate each element of the array.
+--
+-- > truncateA arr == map truncate arr
+--
+-- @since 0.1.0
 truncateA
-  :: (Index ix, Numeric r e, RealFrac a, Integral e)
-  => Array r ix a -> Array r ix e
-truncateA = unsafeLiftArray truncate
+  :: (Source r ix a, RealFrac a, Integral e)
+  => Array r ix a -> Array D ix e
+truncateA = A.map truncate
 {-# INLINE truncateA #-}
 
 
-roundA :: (Index ix, Numeric r e, RealFrac a, Integral e) => Array r ix a -> Array r ix e
-roundA = unsafeLiftArray round
+-- | Round each element of the array.
+--
+-- > truncateA arr == map truncate arr
+--
+-- @since 0.1.0
+roundA :: (Source r ix a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
+roundA = A.map round
 {-# INLINE roundA #-}
 
 
-ceilingA :: (Index ix, Numeric r e, RealFrac a, Integral e) => Array r ix a -> Array r ix e
-ceilingA = unsafeLiftArray ceiling
+-- | Ceiling of each element of the array.
+--
+-- > truncateA arr == map truncate arr
+--
+-- @since 0.1.0
+ceilingA :: (Source r ix a, RealFrac a, Integral e) => Array r ix a -> Array D ix e
+ceilingA = A.map ceiling
 {-# INLINE ceilingA #-}
 
 
-floorA :: (Index ix, Numeric r e, RealFrac a, Integral e) => Array r ix a -> Array r ix e
-floorA = unsafeLiftArray floor
+-- | Floor each element of the array.
+--
+-- > truncateA arr == map truncate arr
+--
+-- @since 0.1.0
+floorA :: (Source r ix 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 ::
      (Load r ix e, Numeric r e, RealFloat e, MonadThrow m)
   => Array r ix e
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
@@ -153,10 +153,12 @@
 --
 -- @since 0.1.0
 lazyFoldlS :: Source r ix e => (a -> e -> a) -> a -> Array r ix e -> a
-lazyFoldlS f initAcc arr = go initAcc 0 where
+lazyFoldlS f initAcc arr = go initAcc 0
+  where
     len = totalElem (size arr)
-    go acc k | k < len = go (f acc (unsafeLinearIndex arr k)) (k + 1)
-             | otherwise = acc
+    go acc k
+      | k < len = go (f acc (unsafeLinearIndex arr k)) (k + 1)
+      | otherwise = acc
 {-# INLINE lazyFoldlS #-}
 
 
diff --git a/src/Data/Massiv/Core/Exception.hs b/src/Data/Massiv/Core/Exception.hs
--- a/src/Data/Massiv/Core/Exception.hs
+++ b/src/Data/Massiv/Core/Exception.hs
@@ -1,9 +1,11 @@
 {-# LANGUAGE CPP #-}
 {-# LANGUAGE GADTs #-}
+{-# LANGUAGE LambdaCase #-}
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 module Data.Massiv.Core.Exception
   ( ImpossibleException(..)
   , throwImpossible
+  , throwEither
   , Uninitialized(..)
   , guardNumberOfElements
   , Exception(..)
@@ -32,6 +34,13 @@
 throwImpossible :: Exception e => e -> a
 throwImpossible = throw . ImpossibleException . toException
 {-# NOINLINE throwImpossible #-}
+
+throwEither :: Either SomeException a -> a
+throwEither =
+  \case
+    Left exc -> throw exc
+    Right res -> res
+{-# INLINE throwEither #-}
 
 instance Exception ImpossibleException where
   displayException (ImpossibleException exc) =
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
@@ -21,23 +21,33 @@
 
 class Num e => Numeric r e where
 
-  {-# MINIMAL unsafeLiftArray, unsafeLiftArray2 #-}
+  {-# MINIMAL foldArray, unsafeLiftArray, unsafeLiftArray2 #-}
 
-  -- sumArray :: Array r Ix1 e -> e
-  -- default sumArray :: Source r Ix1 e => Array r Ix1 e -> e
-  -- sumArray = foldlS (+) 0
-  -- {-# INLINE sumArray #-}
+  -- | Compute sum of all elements in the array
+  --
+  -- @since 0.5.6
+  sumArray :: Index ix => Array r ix e -> e
+  sumArray = foldArray (+) 0
+  {-# INLINE sumArray #-}
 
-  -- productArray :: Array r Ix1 e -> e
-  -- default productArray :: Source r Ix1 e => Array r Ix1 e -> e
-  -- productArray = foldlS (*) 1
-  -- {-# INLINE productArray #-}
+  -- | Compute product of all elements in the array
+  --
+  -- @since 0.5.6
+  productArray :: Index ix => Array r ix e -> e
+  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
 
-  -- unsafeDotProduct :: Array r Ix1 e -> Array r Ix1 e -> e
+  -- | Compute dot product without any extraneous checks
+  --
+  -- @since 0.5.6
+  unsafeDotProduct :: Index ix => Array r ix e -> Array r ix e -> e
+  unsafeDotProduct v1 v2 = sumArray $ unsafeLiftArray2 (*) v1 v2
+  {-# INLINE unsafeDotProduct #-}
 
+
   plusScalar :: Index ix => Array r ix e -> e -> Array r ix e
   plusScalar arr e = unsafeLiftArray (+ e) arr
   {-# INLINE plusScalar #-}
@@ -46,6 +56,10 @@
   minusScalar arr e = unsafeLiftArray (subtract e) arr
   {-# INLINE minusScalar #-}
 
+  scalarMinus :: Index ix => e -> Array r ix e -> Array r ix e
+  scalarMinus e arr = unsafeLiftArray (e -) arr
+  {-# INLINE scalarMinus #-}
+
   multiplyScalar :: Index ix => Array r ix e -> e -> Array r ix e
   multiplyScalar arr e = unsafeLiftArray (* e) arr
   {-# INLINE multiplyScalar #-}
@@ -71,10 +85,14 @@
   powerPointwise arr pow = unsafeLiftArray (^ pow) arr
   {-# INLINE powerPointwise #-}
 
+  -- | Fold over an array
+  --
+  -- @since 0.5.6
+  foldArray :: Index ix => (e -> e -> e) -> e -> Array r ix e -> e
 
-  unsafeLiftArray :: Index ix => (a -> e) -> Array r ix a -> Array r ix e
+  unsafeLiftArray :: Index ix => (e -> e) -> Array r ix e -> Array r ix e
 
-  unsafeLiftArray2 :: Index ix => (a -> b -> e) -> Array r ix a -> Array r ix b -> Array r ix e
+  unsafeLiftArray2 :: Index ix => (e -> e -> e) -> Array r ix e -> Array r ix e -> Array r ix e
 
 
 
@@ -83,6 +101,10 @@
   divideScalar :: Index ix => Array r ix e -> e -> Array r ix e
   divideScalar arr e = unsafeLiftArray (/ e) arr
   {-# INLINE divideScalar #-}
+
+  scalarDivide :: Index ix => e -> Array r ix e -> Array r ix e
+  scalarDivide e arr = unsafeLiftArray (e /) arr
+  {-# INLINE scalarDivide #-}
 
   divisionPointwise :: Index ix => Array r ix e -> Array r ix e -> Array r ix e
   divisionPointwise = unsafeLiftArray2 (/)