diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,40 @@
+# 0.6.0
+
+* Fix semantics of `Applicative`, `Num` and `Fractional` instance for `D` arrays:
+  mismatched sizes will throw an error.
+* 20% speed improvement of matrix multiplication: `multiplyMatrices`, `.><.` and
+  `!><!`. Type signature has changed to `Mutable` for both arguments, thus it's a breaking
+  change.
+* Switch `><.` and `><!` from returning a delayed array to mutable, since that's what
+  `multiplyVectorByMatrix` returns.
+* Addition of synonym `HighIxN` and removing redundant `1 <= n` constraint.
+* Deprecating `makeStencilDef`, `unsafeMapStencil` and fix dangers of invalid stencils
+  reading out of bounds. Get rid of `Value`. Fix for
+  [#109](https://github.com/lehins/massiv/issues/109).
+* Addition of `appComp`
+* Addition of `mkSzM`
+* Addition of `SizeOverflowException` and `SizeNegativeException`
+* Fix setting computation for boxed vector when converted with `fromVectorM` and `fromVector'`
+* Add computation strategy argument to `fromUnboxedVector`, just so it matches other
+  vector conversion functions.
+* Removed `defaultElement`
+* Removed deprecated functions: `#>`, `|*|`, `multiplyTransposed`, `fromIntegerA`,
+  `fromRationalA`, `piA`
+* Addition of `BL` representation and related functionality, fix for [#111](https://github.com/lehins/massiv/issues/111).
+  * Addition of functions: `wrapLazyArray`, `unwrapLazyArray`, `toLazyArray`,
+    `evalLazyArray`, `forceLazyArray`, `unwrapMutableLazyArray`, `fromBoxedVector`,
+    `fromBoxedMVector`.
+  * Rename:
+    * `unsafeNormalBoxedArray` -> `coerceNormalBoxedArray`
+    * `unsafeBoxedArray` -> `coerceBoxedArray`
+  * Remove `unsafeFromBoxedVector`
+  * Conversion from vector with `castFromVector` will return `BL` representation for boxed vector
+  * Change type `B` -> `BL` for functions: `toBoxedVector` and `toBoxedMVector`
+* Rename `N` -> `BN` and add backwards compatibility shim.
+* Make `replicate` a function in `Construct` class
+* Add `newMArray`, `newMArray'` and deprecate `new`
+* Add custom implementation for `<$` in `Functor` instances for `BL` and `B`.
+
 # 0.5.9
 
 * Add `mallocCompute`, `mallocCopy` and `unsafeMallocMArray`
@@ -198,6 +235,8 @@
   * `generateArrayS`
 * Redefined most of the numeric operators with `Numeric` and `NumericFloat`. Will be
   required for SIMD operations.
+* `Num`, `Fractional` and `Applicative` for `D` changed behavior: instead of treating
+  singleton as a special array of any size it is treated as singleton.
 
 # 0.3.6
 
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright Alexey Kuleshevich (c) 2017-2020
+Copyright Alexey Kuleshevich (c) 2017-2021
 
 All rights reserved.
 
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.5.9.0
+version:             0.6.0.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
@@ -7,8 +7,8 @@
 license-file:        LICENSE
 author:              Alexey Kuleshevich
 maintainer:          alexey@kuleshevi.ch
-copyright:           2018-2020 Alexey Kuleshevich
-category:            Data, Data Structures, Parallelism
+copyright:           2018-2021 Alexey Kuleshevich
+category:            Array, Data, Data Structures, Parallelism
 build-type:          Simple
 extra-source-files:  README.md
                    , CHANGELOG.md
@@ -79,7 +79,6 @@
                      , Data.Massiv.Vector.Unsafe
   build-depends:       base >= 4.9 && < 5
                      , bytestring
-                     , data-default-class
                      , deepseq
                      , exceptions
                      , scheduler >= 1.5.0
diff --git a/src/Data/Massiv/Array.hs b/src/Data/Massiv/Array.hs
--- a/src/Data/Massiv/Array.hs
+++ b/src/Data/Massiv/Array.hs
@@ -1,7 +1,7 @@
 {-# OPTIONS_GHC -fno-warn-duplicate-exports #-}
 -- |
 -- Module      : Data.Massiv.Array
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -20,9 +20,11 @@
 --         element is a pointer to the actual value, therefore it is also the slowest
 --         representation. Elements are kept in a Weak Head Normal Form (WHNF).
 --
--- * `N` - Similar to `B`, is also a boxed type, except it's elements are always kept in a Normal
+-- * `BN` - Similar to `B`, is also a boxed type, except it's elements are always kept in a Normal
 --         Form (NF). This property is very useful for parallel processing, i.e. when calling
 --         `compute` you do want all of your elements to be fully evaluated.
+--
+-- * `BL` - Similar to `B`, is also a boxed type, but lazy. It's elements are not evaluated.
 --
 -- * `S` - Is a type of array that is backed by pinned memory, therefore pointers to those arrays
 --         can be passed to FFI calls, because Garbage Collector (GC) is guaranteed not to move
diff --git a/src/Data/Massiv/Array/Delayed.hs b/src/Data/Massiv/Array/Delayed.hs
--- a/src/Data/Massiv/Array/Delayed.hs
+++ b/src/Data/Massiv/Array/Delayed.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      : Data.Massiv.Array.Delayed
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Delayed/Interleaved.hs b/src/Data/Massiv/Array/Delayed/Interleaved.hs
--- a/src/Data/Massiv/Array/Delayed/Interleaved.hs
+++ b/src/Data/Massiv/Array/Delayed/Interleaved.hs
@@ -1,13 +1,13 @@
-{-# LANGUAGE BangPatterns               #-}
-{-# LANGUAGE FlexibleContexts           #-}
-{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE MultiParamTypeClasses      #-}
-{-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Interleaved
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -8,7 +8,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Pull
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -21,8 +21,10 @@
   , eqArrays
   , compareArrays
   , imap
+  , liftArray2Matching
   ) where
 
+import Control.Applicative
 import qualified Data.Foldable as F
 import Data.Massiv.Array.Ops.Fold.Internal as A
 import Data.Massiv.Vector.Stream as S (steps)
@@ -116,10 +118,8 @@
 instance Index ix => Applicative (Array D ix) where
   pure = singleton
   {-# INLINE pure #-}
-  (<*>) (DArray c1 (SafeSz sz1) uIndex1) (DArray c2 (SafeSz sz2) uIndex2) =
-    DArray (c1 <> c2) (SafeSz (liftIndex2 min sz1 sz2)) $ \ !ix ->
-      uIndex1 ix (uIndex2 ix)
-  {-# INLINE (<*>) #-}
+  liftA2 = liftArray2Matching
+  {-# INLINE liftA2 #-}
 
 
 -- | Row-major sequential folding over a Delayed array.
@@ -168,11 +168,11 @@
 {-# INLINE imap #-}
 
 instance (Index ix, Num e) => Num (Array D ix e) where
-  (+)         = unsafeLiftArray2 (+)
+  (+)         = liftArray2Matching (+)
   {-# INLINE (+) #-}
-  (-)         = unsafeLiftArray2 (-)
+  (-)         = liftArray2Matching (-)
   {-# INLINE (-) #-}
-  (*)         = unsafeLiftArray2 (*)
+  (*)         = liftArray2Matching (*)
   {-# INLINE (*) #-}
   abs         = unsafeLiftArray abs
   {-# INLINE abs #-}
@@ -182,7 +182,7 @@
   {-# INLINE fromInteger #-}
 
 instance (Index ix, Fractional e) => Fractional (Array D ix e) where
-  (/)          = unsafeLiftArray2 (/)
+  (/)          = liftArray2Matching (/)
   {-# INLINE (/) #-}
   fromRational = singleton . fromRational
   {-# INLINE fromRational #-}
@@ -217,24 +217,18 @@
   {-# INLINE acosh #-}
 
 
-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
+instance Num e => FoldNumeric D e where
+  unsafeDotProduct = defaultUnsafeDotProduct
   {-# 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
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
   {-# INLINE foldArray #-}
+
+instance Num e => Numeric D e where
   unsafeLiftArray f arr = arr {dIndex = f . dIndex arr}
   {-# INLINE unsafeLiftArray #-}
-  unsafeLiftArray2 f a1 a2 =
+  unsafeLiftArray2 f a1 a2 = -- TODO: possibly use the first size, it is unsafe anyways.
     DArray (dComp a1 <> dComp a2) (SafeSz (liftIndex2 min (unSz (dSize a1)) (unSz (dSize a2)))) $ \i ->
       f (dIndex a1 i) (dIndex a2 i)
   {-# INLINE unsafeLiftArray2 #-}
@@ -278,6 +272,22 @@
     (DArray (getComp arr1 <> getComp arr2) (size arr1) $ \ix ->
        f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
 {-# INLINE compareArrays #-}
+
+
+liftArray2Matching
+  :: (Source r1 ix a, Source r2 ix b)
+  => (a -> b -> e) -> Array r1 ix a -> Array r2 ix b -> Array D ix e
+liftArray2Matching f !arr1 !arr2
+  | sz1 == sz2 =
+    DArray
+      (getComp arr1 <> getComp arr2)
+      sz1
+      (\ !ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
+  | otherwise = throw $ SizeMismatchException (size arr1) (size arr2)
+  where
+    sz1 = size arr1
+    sz2 = size arr2
+{-# INLINE liftArray2Matching #-}
 
 
 -- -- | The usual map.
diff --git a/src/Data/Massiv/Array/Delayed/Push.hs b/src/Data/Massiv/Array/Delayed/Push.hs
--- a/src/Data/Massiv/Array/Delayed/Push.hs
+++ b/src/Data/Massiv/Array/Delayed/Push.hs
@@ -11,7 +11,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Push
--- Copyright   : (c) Alexey Kuleshevich 2019
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -64,6 +64,8 @@
         splitLinearlyWithStartAtM_ scheduler startAt (totalElem sz) (pure . f) dlWrite
       {-# INLINE load #-}
   {-# INLINE makeArrayLinear #-}
+  replicate comp !sz !e = makeLoadArray comp sz e $ \_ _ -> pure ()
+  {-# INLINE replicate #-}
 
 instance Index ix => Resize DL ix where
   unsafeResize !sz arr = arr { dlSize = sz }
@@ -315,8 +317,6 @@
   {-# INLINE getComp #-}
   loadArrayWithSetM scheduler DLArray {dlLoad} = dlLoad scheduler 0
   {-# INLINE loadArrayWithSetM #-}
-  defaultElement _ = Nothing
-  {-# INLINE defaultElement #-}
 
 instance Index ix => Functor (Array DL ix) where
   fmap f arr = arr {dlLoad = loadFunctor arr f}
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Stream
--- Copyright   : (c) Alexey Kuleshevich 2019
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Delayed/Windowed.hs b/src/Data/Massiv/Array/Delayed/Windowed.hs
--- a/src/Data/Massiv/Array/Delayed/Windowed.hs
+++ b/src/Data/Massiv/Array/Delayed/Windowed.hs
@@ -10,7 +10,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Delayed.Windowed
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Manifest.hs b/src/Data/Massiv/Array/Manifest.hs
--- a/src/Data/Massiv/Array/Manifest.hs
+++ b/src/Data/Massiv/Array/Manifest.hs
@@ -3,36 +3,45 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE PatternSynonyms #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
 -- Portability : non-portable
 --
 module Data.Massiv.Array.Manifest
-  (
-  -- * Manifest
+  ( -- * Manifest
     Manifest
   , toManifest
   , M
   -- * Boxed
   , B(..)
-  , N(..)
+  , BL(..)
+  , BN(..)
+  , N
+  , pattern N
   , Uninitialized(..)
   -- ** Access
   , findIndex
   -- ** Conversion
   -- $boxed_conversion_note
+  , toLazyArray
+  , evalLazyArray
+  , forceLazyArray
   , unwrapNormalForm
   , evalNormalForm
   -- *** Primitive Boxed Array
+  , unwrapLazyArray
+  , wrapLazyArray
   , unwrapArray
   , evalArray
   , unwrapMutableArray
+  , unwrapMutableLazyArray
   , evalMutableArray
   , unwrapNormalFormArray
   , evalNormalFormArray
@@ -41,6 +50,8 @@
   -- *** Boxed Vector
   , toBoxedVector
   , toBoxedMVector
+  , fromBoxedVector
+  , fromBoxedMVector
   , evalBoxedVector
   , evalBoxedMVector
   -- * Primitive
@@ -72,7 +83,7 @@
   , mallocCompute
   , mallocCopy
   -- ** Conversion
-  -- *** Primitive Vector
+  -- *** Storable Vector
   , toStorableVector
   , toStorableMVector
   , fromStorableVector
@@ -83,6 +94,7 @@
   , U(..)
   , Unbox
   -- ** Conversion
+  -- *** Unboxed Vector
   , toUnboxedVector
   , toUnboxedMVector
   , fromUnboxedVector
@@ -97,7 +109,6 @@
   ) where
 
 import Control.Monad
-import Data.Massiv.Array.Mutable
 import Data.ByteString as S hiding (findIndex)
 import Data.ByteString.Builder
 import Data.ByteString.Internal
@@ -107,6 +118,7 @@
 import Data.Massiv.Array.Manifest.Primitive
 import Data.Massiv.Array.Manifest.Storable
 import Data.Massiv.Array.Manifest.Unboxed
+import Data.Massiv.Array.Mutable
 import Data.Massiv.Array.Ops.Fold
 import Data.Massiv.Core.Common
 import Data.Word (Word8)
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
@@ -5,12 +5,13 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MagicHash #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE PatternSynonyms #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Boxed
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -18,13 +19,22 @@
 --
 module Data.Massiv.Array.Manifest.Boxed
   ( B(..)
-  , N(..)
+  , BL(..)
+  , BN(..)
+  , N
+  , pattern N
   , Array(..)
+  , wrapLazyArray
+  , unwrapLazyArray
   , unwrapNormalForm
   , evalNormalForm
   , unwrapArray
   , evalArray
+  , toLazyArray
+  , evalLazyArray
+  , forceLazyArray
   , unwrapMutableArray
+  , unwrapMutableLazyArray
   , evalMutableArray
   , unwrapNormalFormArray
   , evalNormalFormArray
@@ -32,13 +42,14 @@
   , evalNormalFormMutableArray
   , toBoxedVector
   , toBoxedMVector
+  , fromBoxedVector
+  , fromBoxedMVector
   , evalBoxedVector
   , evalBoxedMVector
   , evalNormalBoxedVector
   , evalNormalBoxedMVector
-  , unsafeBoxedArray
-  , unsafeNormalBoxedArray
-  , unsafeFromBoxedVector
+  , coerceBoxedArray
+  , coerceNormalBoxedArray
   , seqArray
   , deepseqArray
   ) where
@@ -47,24 +58,24 @@
 import Control.Exception
 import Control.Monad ((>=>))
 import Control.Monad.Primitive
-import Control.Monad.ST (runST)
 import qualified Data.Foldable as F (Foldable(..))
 import Data.Massiv.Array.Delayed.Push (DL)
 import Data.Massiv.Array.Delayed.Stream (DS)
 import Data.Massiv.Array.Manifest.Internal (M, computeAs, toManifest)
 import Data.Massiv.Array.Manifest.List as L
-import Data.Massiv.Vector.Stream as S (steps, isteps)
 import Data.Massiv.Array.Mutable
 import Data.Massiv.Array.Ops.Fold
 import Data.Massiv.Array.Ops.Fold.Internal
 import Data.Massiv.Array.Ops.Map (traverseA)
 import Data.Massiv.Core.Common
 import Data.Massiv.Core.List
+import Data.Massiv.Core.Operations
+import Data.Massiv.Vector.Stream as S (isteps, steps)
 import qualified Data.Primitive.Array as A
 import qualified Data.Vector as VB
 import qualified Data.Vector.Mutable as MVB
 import GHC.Exts as GHC
-import Prelude hiding (mapM)
+import Prelude hiding (mapM, replicate)
 import System.IO.Unsafe (unsafePerformIO)
 
 #include "massiv.h"
@@ -79,37 +90,232 @@
 sizeofMutableArray (A.MutableArray ma#) = I# (sizeofMutableArray# ma#)
 #endif
 
-------------------
--- Boxed Strict --
-------------------
+----------------
+-- Boxed Lazy --
+----------------
 
--- | Array representation for Boxed elements. This structure is element and
--- spine strict, but elements are strict to Weak Head Normal Form (WHNF) only.
-data B = B deriving Show
+-- | Array representation for Boxed elements. This data structure is lazy with respect to
+-- its elements, but is strict with respect to the spine.
+data BL = BL deriving Show
 
-data instance Array B ix e = BArray { bComp   :: !Comp
-                                    , bSize   :: !(Sz ix)
-                                    , bOffset :: {-# UNPACK #-} !Int
-                                    , bData   :: {-# UNPACK #-} !(A.Array e)
-                                    }
+data instance Array BL ix e = BLArray { blComp   :: !Comp
+                                      , blSize   :: !(Sz ix)
+                                      , blOffset :: {-# UNPACK #-} !Int
+                                      , blData   :: {-# UNPACK #-} !(A.Array e)
+                                      }
 
-instance (Ragged L ix e, Show e) => Show (Array B ix e) where
+instance (Ragged L ix e, Show e) => Show (Array BL ix e) where
   showsPrec = showsArrayPrec id
   showList = showArrayList
 
 instance (Ragged L ix e, Show e) => Show (Array DL ix e) where
-  showsPrec = showsArrayPrec (computeAs B)
+  showsPrec = showsArrayPrec (computeAs BL)
   showList = showArrayList
 
 instance Show e => Show (Array DS Ix1 e) where
-  showsPrec = showsArrayPrec (computeAs B)
+  showsPrec = showsArrayPrec (computeAs BL)
   showList = showArrayList
 
 
-instance (Index ix, NFData e) => NFData (Array B ix e) where
+instance (Index ix, NFData e) => NFData (Array BL ix e) where
   rnf = (`deepseqArray` ())
   {-# INLINE rnf #-}
 
+instance (Index ix, Eq e) => Eq (Array BL ix e) where
+  (==) = eqArrays (==)
+  {-# INLINE (==) #-}
+
+instance (Index ix, Ord e) => Ord (Array BL ix e) where
+  compare = compareArrays compare
+  {-# INLINE compare #-}
+
+instance Index ix => Construct BL ix e where
+  setComp c arr = arr { blComp = c }
+  {-# INLINE setComp #-}
+
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
+  {-# INLINE makeArrayLinear #-}
+
+  replicate comp sz e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
+
+instance Index ix => Source BL ix e where
+  unsafeLinearIndex (BLArray _ _sz o a) i =
+    INDEX_CHECK("(Source BL ix e).unsafeLinearIndex",
+                SafeSz . sizeofArray, A.indexArray) a (i + o)
+  {-# INLINE unsafeLinearIndex #-}
+
+  unsafeLinearSlice i k (BLArray c _ o a) = BLArray c k (o + i) a
+  {-# INLINE unsafeLinearSlice #-}
+
+
+instance Index ix => Resize BL ix where
+  unsafeResize !sz !arr = arr { blSize = sz }
+  {-# INLINE unsafeResize #-}
+
+instance Index ix => Extract BL ix e where
+  unsafeExtract !sIx !newSz !arr = unsafeExtract sIx newSz (toManifest arr)
+  {-# INLINE unsafeExtract #-}
+
+
+instance ( Index ix
+         , Index (Lower ix)
+         , Elt M ix e ~ Array M (Lower ix) e
+         , Elt BL ix e ~ Array M (Lower ix) e
+         ) =>
+         OuterSlice BL ix e where
+  unsafeOuterSlice arr = unsafeOuterSlice (toManifest arr)
+  {-# INLINE unsafeOuterSlice #-}
+
+instance ( Index ix
+         , Index (Lower ix)
+         , Elt M ix e ~ Array M (Lower ix) e
+         , Elt BL ix e ~ Array M (Lower ix) e
+         ) =>
+         InnerSlice BL ix e where
+  unsafeInnerSlice arr = unsafeInnerSlice (toManifest arr)
+  {-# INLINE unsafeInnerSlice #-}
+
+instance {-# OVERLAPPING #-} Slice BL Ix1 e where
+  unsafeSlice arr i _ _ = pure (unsafeLinearIndex arr i)
+  {-# INLINE unsafeSlice #-}
+
+
+instance Index ix => Manifest BL ix e where
+
+  unsafeLinearIndexM (BLArray _ _sz o a) i =
+    INDEX_CHECK("(Manifest BL ix e).unsafeLinearIndexM",
+                SafeSz . sizeofArray, A.indexArray) a (i + o)
+  {-# INLINE unsafeLinearIndexM #-}
+
+
+instance Index ix => Mutable BL ix e where
+  data MArray s BL ix e = MBLArray !(Sz ix) {-# UNPACK #-} !Int {-# UNPACK #-} !(A.MutableArray s e)
+
+  msize (MBLArray sz _ _) = sz
+  {-# INLINE msize #-}
+
+  unsafeThaw (BLArray _ sz o a) = MBLArray sz o <$> A.unsafeThawArray a
+  {-# INLINE unsafeThaw #-}
+
+  unsafeFreeze comp (MBLArray sz o ma) = BLArray comp sz o <$> A.unsafeFreezeArray ma
+  {-# INLINE unsafeFreeze #-}
+
+  unsafeNew sz = MBLArray sz 0 <$> A.newArray (totalElem sz) uninitialized
+  {-# INLINE unsafeNew #-}
+
+  initialize _ = return ()
+  {-# INLINE initialize #-}
+
+  newMArray sz e = MBLArray sz 0 <$> A.newArray (totalElem sz) e
+  {-# INLINE newMArray #-}
+
+  unsafeLinearRead (MBLArray _ o ma) i =
+    INDEX_CHECK("(Mutable BL ix e).unsafeLinearRead",
+                SafeSz . sizeofMutableArray, A.readArray) ma (i + o)
+  {-# INLINE unsafeLinearRead #-}
+
+  unsafeLinearWrite (MBLArray _sz o ma) i e = e `seq`
+    INDEX_CHECK("(Mutable BL ix e).unsafeLinearWrite",
+                SafeSz . sizeofMutableArray, A.writeArray) ma (i + o) e
+  {-# INLINE unsafeLinearWrite #-}
+
+instance Index ix => Load BL ix e where
+  type R BL = M
+  size = blSize
+  {-# INLINE size #-}
+  getComp = blComp
+  {-# INLINE getComp #-}
+  loadArrayM !scheduler !arr = splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
+  {-# INLINE loadArrayM #-}
+
+instance Index ix => StrideLoad BL ix e
+
+instance Index ix => Stream BL ix e where
+  toStream = S.steps
+  {-# INLINE toStream #-}
+  toStreamIx = S.isteps
+  {-# INLINE toStreamIx #-}
+
+
+-- | Row-major sequential folding over a Boxed array.
+instance Index ix => Foldable (Array BL ix) where
+  fold = fold
+  {-# INLINE fold #-}
+  foldMap = foldMono
+  {-# INLINE foldMap #-}
+  foldl = lazyFoldlS
+  {-# INLINE foldl #-}
+  foldl' = foldlS
+  {-# INLINE foldl' #-}
+  foldr = foldrFB
+  {-# INLINE foldr #-}
+  foldr' = foldrS
+  {-# INLINE foldr' #-}
+  null (BLArray _ sz _ _) = totalElem sz == 0
+  {-# INLINE null #-}
+  length = totalElem . size
+  {-# INLINE length #-}
+  toList arr = build (\ c n -> foldrFB c n arr)
+  {-# INLINE toList #-}
+
+
+instance Index ix => Functor (Array BL ix) where
+  fmap f arr = makeArrayLinear (blComp arr) (blSize arr) (f . unsafeLinearIndex arr)
+  {-# INLINE fmap #-}
+  (<$) e arr = replicate (getComp arr) (size arr) e
+  {-# INLINE (<$) #-}
+
+instance Index ix => Traversable (Array BL ix) where
+  traverse = traverseA
+  {-# INLINE traverse #-}
+
+instance ( IsList (Array L ix e)
+         , Nested LN ix e
+         , Nested L ix e
+         , Ragged L ix e
+         ) =>
+         IsList (Array BL ix e) where
+  type Item (Array BL ix e) = Item (Array L ix e)
+  fromList = L.fromLists' Seq
+  {-# INLINE fromList #-}
+  toList = GHC.toList . toListArray
+  {-# INLINE toList #-}
+
+instance Num e => FoldNumeric BL e where
+  unsafeDotProduct = defaultUnsafeDotProduct
+  {-# INLINE unsafeDotProduct #-}
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
+  {-# INLINE foldArray #-}
+
+instance Num e => Numeric BL e where
+  unsafeLiftArray = defaultUnsafeLiftArray
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
+  {-# INLINE unsafeLiftArray2 #-}
+
+
+
+------------------
+-- Boxed Strict --
+------------------
+
+-- | Array representation for Boxed elements. This structure is element and
+-- spine strict, but elements are strict to Weak Head Normal Form (WHNF) only.
+data B = B deriving Show
+
+newtype instance Array B ix e = BArray (Array BL ix e)
+
+instance (Ragged L ix e, Show e) => Show (Array B ix e) where
+  showsPrec = showsArrayPrec id
+  showList = showArrayList
+
+instance (Index ix, NFData e) => NFData (Array B ix e) where
+  rnf = (`deepseqArray` ()) . coerce
+  {-# INLINE rnf #-}
+
 instance (Index ix, Eq e) => Eq (Array B ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
@@ -119,24 +325,25 @@
   {-# INLINE compare #-}
 
 instance Index ix => Construct B ix e where
-  setComp c arr = arr { bComp = c }
+  setComp c = coerce (\arr -> arr { blComp = c })
   {-# INLINE setComp #-}
 
-  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (\ !i -> return $! f i)
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
+  replicate comp sz e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
+
 instance Index ix => Source B ix e where
-  unsafeLinearIndex (BArray _ _sz o a) i =
-    INDEX_CHECK("(Source B ix e).unsafeLinearIndex",
-                SafeSz . sizeofArray, A.indexArray) a (i + o)
+  unsafeLinearIndex arr = unsafeLinearIndex (coerce arr :: Array BL ix e)
   {-# INLINE unsafeLinearIndex #-}
 
-  unsafeLinearSlice i k (BArray c _ o a) = BArray c k (o + i) a
+  unsafeLinearSlice i k arr = coerce (unsafeLinearSlice i k (coerce arr :: Array BL ix e))
   {-# INLINE unsafeLinearSlice #-}
 
 
 instance Index ix => Resize B ix where
-  unsafeResize !sz !arr = arr { bSize = sz }
+  unsafeResize sz = coerce (\arr -> arr { blSize = sz })
   {-# INLINE unsafeResize #-}
 
 instance Index ix => Extract B ix e where
@@ -169,47 +376,44 @@
 
 instance Index ix => Manifest B ix e where
 
-  unsafeLinearIndexM (BArray _ _sz o a) i =
-    INDEX_CHECK("(Manifest B ix e).unsafeLinearIndexM",
-                SafeSz . sizeofArray, A.indexArray) a (i + o)
+  unsafeLinearIndexM = coerce unsafeLinearIndexM
   {-# INLINE unsafeLinearIndexM #-}
 
 
 instance Index ix => Mutable B ix e where
-  data MArray s B ix e = MBArray !(Sz ix) {-# UNPACK #-} !Int {-# UNPACK #-} !(A.MutableArray s e)
+  newtype MArray s B ix e = MBArray (MArray s BL ix e)
 
-  msize (MBArray sz _ _) = sz
+  msize = msize . coerce
   {-# INLINE msize #-}
 
-  unsafeThaw (BArray _ sz o a) = MBArray sz o <$> A.unsafeThawArray a
+  unsafeThaw arr = MBArray <$> unsafeThaw (coerce arr)
   {-# INLINE unsafeThaw #-}
 
-  unsafeFreeze comp (MBArray sz o ma) = BArray comp sz o <$> A.unsafeFreezeArray ma
+  unsafeFreeze comp marr = BArray <$> unsafeFreeze comp (coerce marr)
   {-# INLINE unsafeFreeze #-}
 
-  unsafeNew sz = MBArray sz 0 <$> A.newArray (totalElem sz) uninitialized
+  unsafeNew sz = MBArray <$> unsafeNew sz
   {-# INLINE unsafeNew #-}
 
   initialize _ = return ()
   {-# INLINE initialize #-}
 
-  unsafeLinearRead (MBArray _ o ma) i =
-    INDEX_CHECK("(Mutable B ix e).unsafeLinearRead",
-                SafeSz . sizeofMutableArray, A.readArray) ma (i + o)
+  newMArray sz !e = MBArray <$> newMArray sz e
+  {-# INLINE newMArray #-}
+
+  unsafeLinearRead ma = unsafeLinearRead (coerce ma)
   {-# INLINE unsafeLinearRead #-}
 
-  unsafeLinearWrite (MBArray _sz o ma) i e = e `seq`
-    INDEX_CHECK("(Mutable B ix e).unsafeLinearWrite",
-                SafeSz . sizeofMutableArray, A.writeArray) ma (i + o) e
+  unsafeLinearWrite ma i e = e `seq` unsafeLinearWrite (coerce ma) i e
   {-# INLINE unsafeLinearWrite #-}
 
 instance Index ix => Load B ix e where
   type R B = M
-  size = bSize
+  size = blSize . coerce
   {-# INLINE size #-}
-  getComp = bComp
+  getComp = blComp . coerce
   {-# INLINE getComp #-}
-  loadArrayM !scheduler !arr = splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
+  loadArrayM scheduler = coerce (loadArrayM scheduler)
   {-# INLINE loadArrayM #-}
 
 instance Index ix => StrideLoad B ix e
@@ -235,7 +439,7 @@
   {-# INLINE foldr #-}
   foldr' = foldrS
   {-# INLINE foldr' #-}
-  null (BArray _ sz _ _) = totalElem sz == 0
+  null arr = totalElem (size arr) == 0
   {-# INLINE null #-}
   length = totalElem . size
   {-# INLINE length #-}
@@ -244,8 +448,10 @@
 
 
 instance Index ix => Functor (Array B ix) where
-  fmap f arr = makeArrayLinear (bComp arr) (bSize arr) (f . unsafeLinearIndex arr)
+  fmap f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
   {-# INLINE fmap #-}
+  (<$) !e arr = replicate (getComp arr) (size arr) e
+  {-# INLINE (<$) #-}
 
 instance Index ix => Traversable (Array B ix) where
   traverse = traverseA
@@ -263,6 +469,20 @@
   toList = GHC.toList . toListArray
   {-# INLINE toList #-}
 
+instance Num e => FoldNumeric B e where
+  unsafeDotProduct = defaultUnsafeDotProduct
+  {-# INLINE unsafeDotProduct #-}
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
+  {-# INLINE foldArray #-}
+
+instance Num e => Numeric B e where
+  unsafeLiftArray = defaultUnsafeLiftArray
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
+  {-# INLINE unsafeLiftArray2 #-}
+
 -----------------------
 -- Boxed Normal Form --
 -----------------------
@@ -270,52 +490,55 @@
 -- | Array representation for Boxed elements. This structure is element and
 -- spine strict, and elements are always in Normal Form (NF), therefore `NFData`
 -- instance is required.
-data N = N deriving Show
+data BN = BN deriving Show
 
-newtype instance Array N ix e = NArray { bArray :: Array B ix e }
+-- | Type and pattern `N` have been added for backwards compatibility and will be replaced
+-- in the future in favor of `BN`
+type N = BN
+pattern N :: N
+pattern N = BN
+{-# COMPLETE N #-}
 
-instance (Ragged L ix e, Show e, NFData e) => Show (Array N ix e) where
+newtype instance Array N ix e = BNArray { bArray :: Array BL ix e }
+
+instance (Ragged L ix e, Show e, NFData e) => Show (Array BN ix e) where
   showsPrec = showsArrayPrec bArray
   showList = showArrayList
 
-instance (Index ix, NFData e) => NFData (Array N ix e) where
-  rnf (NArray barr) = barr `deepseqArray` ()
+-- | /O(1)/ - `BN` is already in normal form
+instance NFData (Array BN ix e) where
+  rnf = (`seq` ())
   {-# INLINE rnf #-}
 
-instance (Index ix, NFData e, Eq e) => Eq (Array N ix e) where
+instance (Index ix, NFData e, Eq e) => Eq (Array BN ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
 
-instance (Index ix, NFData e, Ord e) => Ord (Array N ix e) where
+instance (Index ix, NFData e, Ord e) => Ord (Array BN ix e) where
   compare = compareArrays compare
   {-# INLINE compare #-}
 
 
-instance (Index ix, NFData e) => Construct N ix e where
-  setComp c (NArray arr) = NArray (arr {bComp = c})
+instance (Index ix, NFData e) => Construct BN ix e where
+  setComp c (BNArray arr) = BNArray (arr {blComp = c})
   {-# INLINE setComp #-}
-  makeArray !comp !sz f =
-    unsafePerformIO $
-    generateArray
-      comp
-      sz
-      (\ !ix ->
-         let res = f ix
-          in res `deepseq` return res)
-  {-# INLINE makeArray #-}
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
+  {-# INLINE makeArrayLinear #-}
+  replicate comp sz e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
 
-instance (Index ix, NFData e) => Source N ix e where
-  unsafeLinearIndex (NArray arr) = unsafeLinearIndex arr
+instance (Index ix, NFData e) => Source BN ix e where
+  unsafeLinearIndex (BNArray arr) = unsafeLinearIndex arr
   {-# INLINE unsafeLinearIndex #-}
-  unsafeLinearSlice i k (NArray a) = NArray $ unsafeLinearSlice i k a
+  unsafeLinearSlice i k (BNArray a) = BNArray $ unsafeLinearSlice i k a
   {-# INLINE unsafeLinearSlice #-}
 
 
-instance Index ix => Resize N ix where
-  unsafeResize !sz = NArray . unsafeResize sz . bArray
+instance Index ix => Resize BN ix where
+  unsafeResize !sz = BNArray . unsafeResize sz . bArray
   {-# INLINE unsafeResize #-}
 
-instance (Index ix, NFData e) => Extract N ix e where
+instance (Index ix, NFData e) => Extract BN ix e where
   unsafeExtract !sIx !newSz !arr = unsafeExtract sIx newSz (toManifest arr)
   {-# INLINE unsafeExtract #-}
 
@@ -324,9 +547,9 @@
          , Index ix
          , Index (Lower ix)
          , Elt M ix e ~ Array M (Lower ix) e
-         , Elt N ix e ~ Array M (Lower ix) e
+         , Elt BN ix e ~ Array M (Lower ix) e
          ) =>
-         OuterSlice N ix e where
+         OuterSlice BN ix e where
   unsafeOuterSlice = unsafeOuterSlice . toManifest
   {-# INLINE unsafeOuterSlice #-}
 
@@ -334,59 +557,62 @@
          , Index ix
          , Index (Lower ix)
          , Elt M ix e ~ Array M (Lower ix) e
-         , Elt N ix e ~ Array M (Lower ix) e
+         , Elt BN ix e ~ Array M (Lower ix) e
          ) =>
-         InnerSlice N ix e where
+         InnerSlice BN ix e where
   unsafeInnerSlice = unsafeInnerSlice . toManifest
   {-# INLINE unsafeInnerSlice #-}
 
-instance {-# OVERLAPPING #-} NFData e => Slice N Ix1 e where
+instance {-# OVERLAPPING #-} NFData e => Slice BN Ix1 e where
   unsafeSlice arr i _ _ = pure (unsafeLinearIndex arr i)
   {-# INLINE unsafeSlice #-}
 
 
-instance (Index ix, NFData e) => Manifest N ix e where
+instance (Index ix, NFData e) => Manifest BN ix e where
 
-  unsafeLinearIndexM (NArray arr) = unsafeLinearIndexM arr
+  unsafeLinearIndexM arr = unsafeLinearIndexM (coerce arr)
   {-# INLINE unsafeLinearIndexM #-}
 
 
-instance (Index ix, NFData e) => Mutable N ix e where
-  newtype MArray s N ix e = MNArray { bmArray :: MArray s B ix e }
+instance (Index ix, NFData e) => Mutable BN ix e where
+  newtype MArray s BN ix e = MBNArray (MArray s BL ix e)
 
-  msize = msize . bmArray
+  msize = msize . coerce
   {-# INLINE msize #-}
 
-  unsafeThaw (NArray arr) = MNArray <$> unsafeThaw arr
+  unsafeThaw arr = MBNArray <$> unsafeThaw (coerce arr)
   {-# INLINE unsafeThaw #-}
 
-  unsafeFreeze comp (MNArray marr) = NArray <$> unsafeFreeze comp marr
+  unsafeFreeze comp marr = BNArray <$> unsafeFreeze comp (coerce marr)
   {-# INLINE unsafeFreeze #-}
 
-  unsafeNew sz = MNArray <$> unsafeNew sz
+  unsafeNew sz = MBNArray <$> unsafeNew sz
   {-# INLINE unsafeNew #-}
 
   initialize _ = return ()
   {-# INLINE initialize #-}
 
-  unsafeLinearRead (MNArray ma) = unsafeLinearRead ma
+  newMArray sz e = e `deepseq` (MBNArray <$> newMArray sz e)
+  {-# INLINE newMArray #-}
+
+  unsafeLinearRead ma = unsafeLinearRead (coerce ma)
   {-# INLINE unsafeLinearRead #-}
 
-  unsafeLinearWrite (MNArray ma) i e = e `deepseq` unsafeLinearWrite ma i e
+  unsafeLinearWrite ma i e = e `deepseq` unsafeLinearWrite (coerce ma) i e
   {-# INLINE unsafeLinearWrite #-}
 
-instance (Index ix, NFData e) => Load N ix e where
-  type R N = M
-  size = bSize . bArray
+instance (Index ix, NFData e) => Load BN ix e where
+  type R BN = M
+  size = blSize . coerce
   {-# INLINE size #-}
-  getComp = bComp . bArray
+  getComp = blComp . coerce
   {-# INLINE getComp #-}
   loadArrayM !scheduler !arr = splitLinearlyWith_ scheduler (elemsCount arr) (unsafeLinearIndex arr)
   {-# INLINE loadArrayM #-}
 
-instance (Index ix, NFData e) => StrideLoad N ix e
+instance (Index ix, NFData e) => StrideLoad BN ix e
 
-instance (Index ix, NFData e) => Stream N ix e where
+instance (Index ix, NFData e) => Stream BN ix e where
   toStream = toStream . coerce
   {-# INLINE toStream #-}
   toStreamIx = toStreamIx . coerce
@@ -399,14 +625,27 @@
          , Nested L ix e
          , Ragged L ix e
          ) =>
-         IsList (Array N ix e) where
-  type Item (Array N ix e) = Item (Array L ix e)
+         IsList (Array BN ix e) where
+  type Item (Array BN ix e) = Item (Array L ix e)
   fromList = L.fromLists' Seq
   {-# INLINE fromList #-}
   toList = GHC.toList . toListArray
   {-# INLINE toList #-}
 
+instance (NFData e, Num e) => FoldNumeric BN e where
+  unsafeDotProduct = defaultUnsafeDotProduct
+  {-# INLINE unsafeDotProduct #-}
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
+  {-# INLINE foldArray #-}
 
+instance (NFData e, Num e) => Numeric BN e where
+  unsafeLiftArray = defaultUnsafeLiftArray
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
+  {-# INLINE unsafeLiftArray2 #-}
+
 ----------------------
 -- Helper functions --
 ----------------------
@@ -423,7 +662,7 @@
 --
 -- @since 0.2.1
 unwrapArray :: Array B ix e -> A.Array e
-unwrapArray = bData
+unwrapArray = blData . coerce
 {-# INLINE unwrapArray #-}
 
 -- | /O(n)/ - Wrap a boxed array and evaluate all elements to a WHNF.
@@ -432,20 +671,66 @@
 evalArray ::
      Comp -- ^ Computation strategy
   -> A.Array e -- ^ Lazy boxed array from @primitive@ package.
-  -> Array B Ix1 e
-evalArray = fromArraySeq (\a -> a `seqArray` a)
+  -> Vector B e
+evalArray comp a = evalLazyArray $ setComp comp $ wrapLazyArray a
 {-# INLINE evalArray #-}
 
 
+-- | /O(1)/ - Unwrap boxed array. This will discard any possible slicing that has been
+-- applied to the array.
+--
+-- @since 0.6.0
+unwrapLazyArray :: Array BL ix e -> A.Array e
+unwrapLazyArray = blData
+{-# INLINE unwrapLazyArray #-}
+
+-- | /O(1)/ - Wrap a boxed array.
+--
+-- @since 0.6.0
+wrapLazyArray :: A.Array e -> Vector BL e
+wrapLazyArray a = BLArray Seq (SafeSz (sizeofArray a)) 0 a
+{-# INLINE wrapLazyArray #-}
+
+
+-- | /O(1)/ - Cast a strict boxed array into a lazy boxed array.
+--
+-- @since 0.6.0
+toLazyArray :: Array B ix e -> Array BL ix e
+toLazyArray = coerce
+{-# INLINE toLazyArray #-}
+
+-- | /O(n)/ - Evaluate all elements of a boxed lazy array to weak head normal form
+--
+-- @since 0.6.0
+evalLazyArray :: Index ix => Array BL ix e -> Array B ix e
+evalLazyArray arr = arr `seqArray` BArray arr
+{-# INLINE evalLazyArray #-}
+
+-- | /O(n)/ - Evaluate all elements of a boxed lazy array to normal form
+--
+-- @since 0.6.0
+forceLazyArray :: (NFData e, Index ix) => Array BL ix e -> Array N ix e
+forceLazyArray arr = arr `deepseqArray` BNArray arr
+{-# INLINE forceLazyArray #-}
+
 -- | /O(1)/ - Unwrap mutable boxed array. This will discard any possible slicing that has been
 -- applied to the array.
 --
 -- @since 0.2.1
 unwrapMutableArray :: MArray s B ix e -> A.MutableArray s e
-unwrapMutableArray (MBArray _ _ marr) = marr
+unwrapMutableArray (MBArray (MBLArray _ _ marr)) = marr
 {-# INLINE unwrapMutableArray #-}
 
 
+-- | /O(1)/ - Unwrap mutable boxed lazy array. This will discard any possible slicing that has been
+-- applied to the array.
+--
+-- @since 0.6.0
+unwrapMutableLazyArray :: MArray s BL ix e -> A.MutableArray s e
+unwrapMutableLazyArray (MBLArray _ _ marr) = marr
+{-# INLINE unwrapMutableLazyArray #-}
+
+
 -- | /O(n)/ - Wrap mutable boxed array and evaluate all elements to WHNF.
 --
 -- @since 0.2.1
@@ -453,7 +738,7 @@
      PrimMonad m
   => A.MutableArray (PrimState m) e -- ^ Mutable array that will get wrapped
   -> m (MArray (PrimState m) B Ix1 e)
-evalMutableArray = fromMutableArraySeq seq
+evalMutableArray = fmap MBArray . fromMutableArraySeq seq
 {-# INLINE evalMutableArray #-}
 
 -------------------
@@ -465,7 +750,7 @@
 --
 -- @since 0.2.1
 unwrapNormalFormArray :: Array N ix e -> A.Array e
-unwrapNormalFormArray = bData . bArray
+unwrapNormalFormArray = blData . coerce
 {-# INLINE unwrapNormalFormArray #-}
 
 -- | /O(n)/ - Wrap a boxed array and evaluate all elements to a Normal Form (NF).
@@ -476,7 +761,7 @@
   => Comp -- ^ Computation strategy
   -> A.Array e -- ^ Lazy boxed array
   -> Array N Ix1 e
-evalNormalFormArray = fromArraySeq (\a -> a `deepseqArray` NArray a)
+evalNormalFormArray comp = forceLazyArray . setComp comp . wrapLazyArray
 {-# INLINE evalNormalFormArray #-}
 
 
@@ -485,7 +770,7 @@
 --
 -- @since 0.2.1
 unwrapNormalFormMutableArray :: MArray s N ix e -> A.MutableArray s e
-unwrapNormalFormMutableArray (MNArray marr) = unwrapMutableArray marr
+unwrapNormalFormMutableArray = unwrapMutableLazyArray . coerce
 {-# INLINE unwrapNormalFormMutableArray #-}
 
 
@@ -496,7 +781,7 @@
      (PrimMonad m, NFData e)
   => A.MutableArray (PrimState m) e
   -> m (MArray (PrimState m) N Ix1 e)
-evalNormalFormMutableArray marr = MNArray <$> fromMutableArraySeq deepseq marr
+evalNormalFormMutableArray marr = MBNArray <$> fromMutableArraySeq deepseq marr
 {-# INLINE evalNormalFormMutableArray #-}
 
 
@@ -508,33 +793,25 @@
      PrimMonad m
   => (e -> m () -> m a)
   -> A.MutableArray (PrimState m) e
-  -> m (MArray (PrimState m) B Ix1 e)
+  -> m (MArray (PrimState m) BL Ix1 e)
 fromMutableArraySeq with ma = do
   let !sz = sizeofMutableArray ma
   loopM_ 0 (< sz) (+ 1) (A.readArray ma >=> (`with` return ()))
-  return $! MBArray (SafeSz sz) 0 ma
+  return $! MBLArray (SafeSz sz) 0 ma
 {-# INLINE fromMutableArraySeq #-}
 
-fromArraySeq ::
-     (Array B Ix1 e -> a)
-  -> Comp
-  -> A.Array e
-  -> a
-fromArraySeq with comp barr = with (BArray comp (SafeSz (sizeofArray barr)) 0 barr)
-{-# INLINE fromArraySeq #-}
 
-
-seqArray :: Index ix => Array B ix a -> t -> t
+seqArray :: Index ix => Array BL ix a -> t -> t
 seqArray !arr t = foldlInternal (flip seq) () (flip seq) () arr `seq` t
 {-# INLINE seqArray #-}
 
 
-deepseqArray :: (NFData a, Index ix) => Array B ix a -> t -> t
+deepseqArray :: (NFData a, Index ix) => Array BL ix a -> t -> t
 deepseqArray !arr t = foldlInternal (flip deepseq) () (flip seq) () arr `seq` t
 {-# INLINE deepseqArray #-}
 
 
--- | /O(n)/ - Compute all elements of a boxed array to NF (normal form)
+-- | /O(1)/ - Converts array from `N` to `B` representation.
 --
 -- @since 0.5.0
 unwrapNormalForm :: Array N ix e -> Array B ix e
@@ -545,21 +822,21 @@
 --
 -- @since 0.5.0
 evalNormalForm :: (Index ix, NFData e) => Array B ix e -> Array N ix e
-evalNormalForm arr = arr `deepseqArray` NArray arr
+evalNormalForm (BArray arr) = arr `deepseqArray` BNArray arr
 {-# INLINE evalNormalForm #-}
 
 -- | /O(1)/ - Converts a boxed `Array` into a `VB.Vector`.
 --
 -- @since 0.5.0
-toBoxedVector :: Index ix => Array B ix a -> VB.Vector a
+toBoxedVector :: Index ix => Array BL ix a -> VB.Vector a
 toBoxedVector arr = runST $ VB.unsafeFreeze . toBoxedMVector =<< unsafeThaw arr
 {-# INLINE toBoxedVector #-}
 
 -- | /O(1)/ - Converts a boxed `MArray` into a `VMB.MVector`.
 --
 -- @since 0.5.0
-toBoxedMVector :: Index ix => MArray s B ix a -> MVB.MVector s a
-toBoxedMVector (MBArray sz o marr) = MVB.MVector o (totalElem sz) marr
+toBoxedMVector :: Index ix => MArray s BL ix a -> MVB.MVector s a
+toBoxedMVector (MBLArray sz o marr) = MVB.MVector o (totalElem sz) marr
 {-# INLINE toBoxedMVector #-}
 
 -- | /O(n)/ - Convert a boxed vector and evaluate all elements to WHNF. Computation
@@ -567,9 +844,7 @@
 --
 -- @since 0.5.0
 evalBoxedVector :: Comp -> VB.Vector a -> Array B Ix1 a
-evalBoxedVector comp v = arr `seqArray` arr
-  where
-    arr = setComp comp $ unsafeFromBoxedVector v
+evalBoxedVector comp = evalLazyArray . setComp comp . fromBoxedVector
 {-# INLINE evalBoxedVector #-}
 
 
@@ -578,51 +853,58 @@
 --
 -- @since 0.5.0
 evalBoxedMVector :: PrimMonad m => MVB.MVector (PrimState m) a -> m (MArray (PrimState m) B Ix1 a)
-evalBoxedMVector (MVB.MVector o k ma) = do
-  let marr = MBArray (SafeSz k) o ma
-  loopM_ o (< k) (+ 1) (A.readArray ma >=> (`seq` pure ()))
-  pure marr
+evalBoxedMVector (MVB.MVector o k ma) =
+  let marr = MBArray (MBLArray (SafeSz k) o ma)
+   in marr <$ loopM_ o (< k) (+ 1) (A.readArray ma >=> (`seq` pure ()))
 {-# INLINE evalBoxedMVector #-}
 
 
--- | /O(n)/ - Cast a boxed vector without touching any elements. It is unsafe because it
--- violates the invariant that all elements of `B` array are in WHNF.
+-- | /O(1)/ - Cast a boxed vector without touching any elements.
 --
--- @since 0.5.0
-unsafeFromBoxedVector :: VB.Vector a -> Array B Ix1 a
-unsafeFromBoxedVector v =
+-- @since 0.6.0
+fromBoxedVector :: VB.Vector a -> Vector BL a
+fromBoxedVector v =
   runST $ do
     MVB.MVector o k ma <- VB.unsafeThaw v
-    unsafeFreeze Seq $ MBArray (SafeSz k) o ma
-{-# INLINE unsafeFromBoxedVector #-}
+    unsafeFreeze Seq $ MBLArray (SafeSz k) o ma
+{-# INLINE fromBoxedVector #-}
 
--- | /O(n)/ - Cast a boxed array. It is unsafe because it violates the invariant that all
--- elements of `N` array are in NF.
+
+-- | /O(1)/ - Convert mutable boxed vector to a lazy mutable boxed array. Both keep
+-- pointing to the same memory
 --
--- @since 0.5.0
-unsafeBoxedArray :: A.Array e -> Array B Ix1 e
-unsafeBoxedArray = fromArraySeq id Seq
-{-# INLINE unsafeBoxedArray #-}
+-- @since 0.6.0
+fromBoxedMVector :: MVB.MVector s a -> MArray s BL Ix1 a
+fromBoxedMVector (MVB.MVector o k ma) = MBLArray (SafeSz k) o ma
+{-# INLINE fromBoxedMVector #-}
 
 
--- | /O(n)/ - Cast a boxed array. It is unsafe because it violates the invariant that all
--- elements of `N` array are in NF.
+-- | /O(1)/ - Cast a boxed lazy array. It is unsafe because it can violate the invariant
+-- that all elements of `N` array are in NF.
 --
--- @since 0.5.0
-unsafeNormalBoxedArray :: Array B ix e -> Array N ix e
-unsafeNormalBoxedArray = coerce
-{-# INLINE unsafeNormalBoxedArray #-}
+-- @since 0.6.0
+coerceNormalBoxedArray :: Array BL ix e -> Array N ix e
+coerceNormalBoxedArray = coerce
+{-# INLINE coerceNormalBoxedArray #-}
 
+
+-- | /O(1)/ - Cast a boxed lazy array. It is unsafe because it can violate the invariant
+-- that all elements of `B` array are in WHNF.
+--
+-- @since 0.6.0
+coerceBoxedArray :: Array BL ix e -> Array B ix e
+coerceBoxedArray = coerce
+{-# INLINE coerceBoxedArray #-}
+
 -- | /O(n)/ - Convert mutable boxed vector and evaluate all elements to WHNF
 -- sequentially. Both keep pointing to the same memory
 --
 -- @since 0.5.0
 evalNormalBoxedMVector ::
      (NFData a, PrimMonad m) => MVB.MVector (PrimState m) a -> m (MArray (PrimState m) N Ix1 a)
-evalNormalBoxedMVector (MVB.MVector o k ma) = do
-  let marr = MNArray (MBArray (SafeSz k) o ma)
-  loopM_ o (< k) (+ 1) (A.readArray ma >=> (`deepseq` pure ()))
-  pure marr
+evalNormalBoxedMVector (MVB.MVector o k ma) =
+  let marr = MBNArray (MBLArray (SafeSz k) o ma)
+   in marr <$ loopM_ o (< k) (+ 1) (A.readArray ma >=> (`deepseq` pure ()))
 {-# INLINE evalNormalBoxedMVector #-}
 
 -- | /O(n)/ - Convert a boxed vector and evaluate all elements to WHNF. Computation
@@ -633,7 +915,6 @@
 evalNormalBoxedVector comp v =
   runST $ do
     MVB.MVector o k ma <- VB.unsafeThaw v
-    arr <- unsafeFreeze comp $ MBArray (SafeSz k) o ma
-    arr `deepseqArray` pure (NArray arr)
+    forceLazyArray <$> unsafeFreeze comp (MBLArray (SafeSz k) o ma)
 {-# INLINE evalNormalBoxedVector #-}
 
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
@@ -10,7 +10,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -55,6 +55,7 @@
 import Data.Massiv.Vector.Stream as S (steps, isteps)
 import Data.Massiv.Core.Common
 import Data.Massiv.Core.List
+import Data.Massiv.Core.Operations
 import Data.Maybe (fromMaybe)
 import Data.Typeable
 import GHC.Base hiding (ord)
@@ -207,6 +208,14 @@
   {-# INLINE toStreamIx #-}
 
 
+instance Num e => FoldNumeric M e where
+  unsafeDotProduct = defaultUnsafeDotProduct
+  {-# INLINE unsafeDotProduct #-}
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
+  {-# INLINE foldArray #-}
+
 -- | Ensure that Array is computed, i.e. represented with concrete elements in memory, hence is the
 -- `Mutable` type class restriction. Use `setComp` if you'd like to change computation strategy
 -- before calling @compute@
@@ -417,7 +426,7 @@
 -- ====__Example__
 --
 -- >>> import Data.Massiv.Array
--- >>> a = computeAs P $ makeLoadArrayS (Sz2 8 8) (0 :: Int) $ \ w -> w (0 :. 0) 1 >> pure ()
+-- >>> a = computeAs P $ makeLoadArrayS (Sz2 8 8) (0 :: Int) $ \ w -> () <$ w (0 :. 0) 1
 -- >>> a
 -- Array P Seq (Sz (8 :. 8))
 --   [ [ 1, 0, 0, 0, 0, 0, 0, 0 ]
@@ -430,7 +439,7 @@
 --   , [ 0, 0, 0, 0, 0, 0, 0, 0 ]
 --   ]
 -- >>> nextPascalRow cur above = if cur == 0 then above else cur
--- >>> pascal = makeStencil (Sz2 2 2) 1 $ \ get -> nextPascalRow <$> get (0 :. 0) <*> get (-1 :. -1) + get (-1 :. 0)
+-- >>> pascal = makeStencil (Sz2 2 2) 1 $ \ get -> nextPascalRow (get (0 :. 0)) (get (-1 :. -1) + get (-1 :. 0))
 -- >>> iterateUntil (\_ _ a -> (a ! (7 :. 7)) /= 0) (\ _ -> mapStencil (Fill 0) pascal) a
 -- Array P Seq (Sz (8 :. 8))
 --   [ [ 1, 0, 0, 0, 0, 0, 0, 0 ]
@@ -461,7 +470,7 @@
       let loadArr = iteration 1 initArr1
       marr <- unsafeNew (size loadArr)
       iterateLoop
-        (\n a a' _ -> pure $ convergence n a a')
+        (\n a comp marr' -> convergence n a <$> unsafeFreeze comp marr')
         iteration
         1
         initArr1
@@ -505,7 +514,7 @@
 
 iterateLoop ::
      (Load r' ix e, Mutable r ix e, PrimMonad m, MonadIO m, PrimState m ~ RealWorld)
-  => (Int -> Array r ix e -> Array r ix e -> MArray (PrimState m) r ix e -> m Bool)
+  => (Int -> Array r ix e -> Comp -> MArray (PrimState m) r ix e -> m Bool)
   -> (Int -> Array r ix e -> Array r' ix e)
   -> Int
   -> Array r ix e
@@ -518,6 +527,7 @@
       let !sz = size loadArr
           !k = totalElem sz
           !mk = totalElem (msize marr)
+          !comp = getComp loadArr
       marr' <-
         if k == mk
           then pure marr
@@ -525,8 +535,8 @@
                  then unsafeLinearShrink marr sz
                  else unsafeLinearGrow marr sz
       computeInto marr' loadArr
-      arr' <- unsafeFreeze (getComp loadArr) marr'
-      shouldStop <- convergence n arr arr' marr'
+      shouldStop <- convergence n arr comp marr'
+      arr' <- unsafeFreeze comp marr'
       if shouldStop
         then pure arr'
         else do
diff --git a/src/Data/Massiv/Array/Manifest/List.hs b/src/Data/Massiv/Array/Manifest/List.hs
--- a/src/Data/Massiv/Array/Manifest/List.hs
+++ b/src/Data/Massiv/Array/Manifest/List.hs
@@ -4,10 +4,9 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeOperators #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.List
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -11,7 +11,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Primitive
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -103,9 +103,12 @@
   setComp c arr = arr { pComp = c }
   {-# INLINE setComp #-}
 
-  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (return . f)
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
   {-# INLINE makeArrayLinear #-}
 
+  replicate comp !sz !e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
+
 instance (Prim e, Index ix) => Source P ix e where
   unsafeLinearIndex _arr@(PArray _ _ o a) i =
     INDEX_CHECK("(Source P ix e).unsafeLinearIndex",
@@ -253,35 +256,18 @@
   {-# 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
+instance (Prim e, Num e) => FoldNumeric P e where
+  unsafeDotProduct = defaultUnsafeDotProduct
   {-# 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
+  powerSumArray = defaultPowerSumArray
   {-# INLINE powerSumArray #-}
-  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
+  foldArray = defaultFoldArray
   {-# INLINE foldArray #-}
-  unsafeLiftArray f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
+
+instance (Prim e, Num e) => Numeric P e where
+  unsafeLiftArray = defaultUnsafeLiftArray
   {-# 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)
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
   {-# INLINE unsafeLiftArray2 #-}
 
 
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
@@ -9,7 +9,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Storable
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -35,6 +35,7 @@
   ) where
 
 import Control.DeepSeq (NFData(..), deepseq)
+import Control.Exception
 import Control.Monad.IO.Unlift
 import Control.Monad.Primitive (unsafePrimToPrim)
 import Data.Massiv.Array.Delayed.Pull (compareArrays, eqArrays)
@@ -50,16 +51,15 @@
 import qualified Data.Vector.Generic.Mutable as VGM
 import qualified Data.Vector.Storable as VS
 import qualified Data.Vector.Storable.Mutable as MVS
-import Foreign.ForeignPtr (withForeignPtr, newForeignPtr)
-import Foreign.Marshal.Array (advancePtr, copyArray)
+import Foreign.ForeignPtr (newForeignPtr, withForeignPtr)
 import Foreign.Marshal.Alloc
+import Foreign.Marshal.Array (advancePtr, copyArray)
 import Foreign.Ptr
 import Foreign.Storable
 import GHC.Exts as GHC (IsList(..))
 import GHC.ForeignPtr (ForeignPtr(..), ForeignPtrContents(..))
 import Prelude hiding (mapM)
 import System.IO.Unsafe (unsafePerformIO)
-import Control.Exception
 
 #include "massiv.h"
 
@@ -91,10 +91,13 @@
   setComp c arr = arr { sComp = c }
   {-# INLINE setComp #-}
 
-  makeArray !comp !sz f = unsafePerformIO $ generateArray comp sz (return . f)
-  {-# INLINE makeArray #-}
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
+  {-# INLINE makeArrayLinear #-}
 
+  replicate comp !sz !e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
 
+
 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
@@ -226,38 +229,20 @@
   {-# INLINE toStreamIx #-}
 
 
-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
+instance (Storable e, Num e) => FoldNumeric S e where
+  unsafeDotProduct = defaultUnsafeDotProduct
   {-# 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
+  powerSumArray = defaultPowerSumArray
   {-# INLINE powerSumArray #-}
-  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
+  foldArray = defaultFoldArray
   {-# INLINE foldArray #-}
-  unsafeLiftArray f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
+
+instance (Storable e, Num e) => Numeric S e where
+  unsafeLiftArray = defaultUnsafeLiftArray
   {-# 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)
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
   {-# INLINE unsafeLiftArray2 #-}
 
-
 instance (Storable e, Floating e) => NumericFloat S e
 
 
@@ -377,10 +362,10 @@
 --
 -- @since 0.5.9
 unsafeMallocMArray ::
-     forall ix e. (Index ix, Storable e)
+     forall ix e m. (Index ix, Storable e, MonadIO m)
   => Sz ix
-  -> IO (MArray RealWorld S ix e)
-unsafeMallocMArray sz = do
+  -> m (MArray RealWorld S ix e)
+unsafeMallocMArray sz = liftIO $ do
   let n = totalElem sz
   foreignPtr <- mask_ $ do
     ptr <- mallocBytes (sizeOf (undefined :: e) * n)
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
@@ -3,13 +3,12 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE NamedFieldPuns #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Unboxed
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -33,6 +32,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.Unboxed as VU
 import qualified Data.Vector.Unboxed.Mutable as MVU
@@ -63,10 +63,13 @@
   setComp c arr = arr { uComp = c }
   {-# INLINE setComp #-}
 
-  makeArray !comp !sz f = unsafePerformIO $ generateArray comp sz (return . f)
-  {-# INLINE makeArray #-}
+  makeArrayLinear !comp !sz f = unsafePerformIO $ generateArrayLinear comp sz (pure . f)
+  {-# INLINE makeArrayLinear #-}
 
+  replicate comp !sz !e = runST (newMArray sz e >>= unsafeFreeze comp)
+  {-# INLINE replicate #-}
 
+
 instance (VU.Unbox e, Eq e, Index ix) => Eq (Array U ix e) where
   (==) = eqArrays (==)
   {-# INLINE (==) #-}
@@ -210,6 +213,21 @@
   {-# INLINE toList #-}
 
 
+instance (VU.Unbox e, Num e) => FoldNumeric U e where
+  unsafeDotProduct = defaultUnsafeDotProduct
+  {-# INLINE unsafeDotProduct #-}
+  powerSumArray = defaultPowerSumArray
+  {-# INLINE powerSumArray #-}
+  foldArray = defaultFoldArray
+  {-# INLINE foldArray #-}
+
+instance (VU.Unbox e, Num e) => Numeric U e where
+  unsafeLiftArray = defaultUnsafeLiftArray
+  {-# INLINE unsafeLiftArray #-}
+  unsafeLiftArray2 = defaultUnsafeLiftArray2
+  {-# INLINE unsafeLiftArray2 #-}
+
+
 -- | /O(1)/ - Unwrap unboxed array and pull out the underlying unboxed vector.
 --
 -- @since 0.2.1
@@ -229,9 +247,9 @@
 
 -- | /O(1)/ - Wrap an unboxed vector and produce an unboxed flat array.
 --
--- @since 0.5.0
-fromUnboxedVector :: VU.Unbox e => VU.Vector e -> Array U Ix1 e
-fromUnboxedVector v = UArray Seq (SafeSz (VU.length v)) v
+-- @since 0.6.0
+fromUnboxedVector :: VU.Unbox e => Comp -> VU.Vector e -> Array U Ix1 e
+fromUnboxedVector comp v = UArray comp (SafeSz (VU.length v)) v
 {-# INLINE fromUnboxedVector #-}
 
 
diff --git a/src/Data/Massiv/Array/Manifest/Vector.hs b/src/Data/Massiv/Array/Manifest/Vector.hs
--- a/src/Data/Massiv/Array/Manifest/Vector.hs
+++ b/src/Data/Massiv/Array/Manifest/Vector.hs
@@ -6,7 +6,7 @@
 {-# LANGUAGE TypeOperators #-}
 -- |
 -- Module      : Data.Massiv.Array.Manifest.Vector
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -43,7 +43,7 @@
   ARepr VU.Vector = U
   ARepr VS.Vector = S
   ARepr VP.Vector = P
-  ARepr VB.Vector = B
+  ARepr VB.Vector = BL
 
 -- | Match array representation to a vector type
 type family VRepr r :: * -> * where
@@ -52,6 +52,7 @@
   VRepr P = VP.Vector
   VRepr B = VB.Vector
   VRepr N = VB.Vector
+  VRepr BL = VB.Vector
 
 
 -- | /O(1)/ - conversion from vector to an array with a corresponding representation. Will
@@ -77,8 +78,7 @@
          return $ PArray {pComp = comp, pSize = sz, pOffset = o, pData = ba}
     , do Refl <- eqT :: Maybe (v :~: VB.Vector)
          bVector <- join $ gcast1 (Just vector)
-         let ba = unsafeFromBoxedVector bVector
-         ba `seqArray` pure (unsafeResize sz ba)
+         pure $ unsafeResize sz $ setComp comp $ fromBoxedVector bVector
     ]
 {-# NOINLINE castFromVector #-}
 
@@ -141,10 +141,13 @@
          return $ VP.Vector (pOffset pArr) (totalElem (size arr)) $ pData pArr
     , do Refl <- eqT :: Maybe (r :~: B)
          bArr <- gcastArr arr
-         return $ toBoxedVector bArr
+         return $ toBoxedVector $ toLazyArray bArr
     , do Refl <- eqT :: Maybe (r :~: N)
          bArr <- gcastArr arr
-         return $ toBoxedVector $ bArray bArr
+         return $ toBoxedVector $ toLazyArray $ unwrapNormalForm bArr
+    , do Refl <- eqT :: Maybe (r :~: BL)
+         bArr <- gcastArr arr
+         return $ toBoxedVector bArr
     ]
 {-# NOINLINE castToVector #-}
 
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
@@ -6,7 +6,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -35,12 +35,14 @@
   , swap'
   -- ** Operations on @MArray@
   -- *** Immutable conversion
-  , new
   , thaw
   , thawS
   , freeze
   , freezeS
   -- *** Create mutable
+  , new
+  , newMArray
+  , newMArray'
   , makeMArray
   , makeMArrayLinear
   , makeMArrayS
@@ -109,10 +111,24 @@
 -- boxed arrays in will be a thunk with `Uninitialized` exception, while for others it will be
 -- simply zeros.
 --
+-- @since 0.1.0
+new ::
+     forall r ix e m. (Mutable r ix e, PrimMonad m)
+  => Sz ix
+  -> m (MArray (PrimState m) r ix e)
+new = initializeNew Nothing
+{-# INLINE new #-}
+{-# DEPRECATED new "In favor of a more robust and safer `newMArray` or a more consistently named `newMArray'`" #-}
+
+
+-- | /O(n)/ - Initialize a new mutable array. All elements will be set to some default value. For
+-- boxed arrays in will be a thunk with `Uninitialized` exception, while for others it will be
+-- simply zeros. This is a partial function.
+--
 -- ==== __Examples__
 --
 -- >>> import Data.Massiv.Array
--- >>> marr <- new (Sz2 2 6) :: IO (MArray RealWorld P Ix2 Int)
+-- >>> marr <- newMArray' (Sz2 2 6) :: IO (MArray RealWorld P Ix2 Int)
 -- >>> freeze Seq marr
 -- Array P Seq (Sz (2 :. 6))
 --   [ [ 0, 0, 0, 0, 0, 0 ]
@@ -122,22 +138,23 @@
 -- Or using @TypeApplications@:
 --
 -- >>> :set -XTypeApplications
--- >>> new @P @Ix2 @Int (Sz2 2 6) >>= freezeS
+-- >>> newMArray' @P @Ix2 @Int (Sz2 2 6) >>= freezeS
 -- Array P Seq (Sz (2 :. 6))
 --   [ [ 0, 0, 0, 0, 0, 0 ]
 --   , [ 0, 0, 0, 0, 0, 0 ]
 --   ]
--- >>> new @B @_ @Int (Sz2 2 6) >>= (`readM` 1)
+-- >>> newMArray' @B @_ @Int (Sz2 2 6) >>= (`readM` 1)
 -- *** Exception: Uninitialized
 --
--- @since 0.1.0
-new ::
+-- @since 0.6.0
+newMArray' ::
      forall r ix e m. (Mutable r ix e, PrimMonad m)
   => Sz ix
   -> m (MArray (PrimState m) r ix e)
-new = initializeNew Nothing
-{-# INLINE new #-}
+newMArray' sz = unsafeNew sz >>= \ma -> ma <$ initialize ma
+{-# INLINE newMArray' #-}
 
+
 -- | /O(n)/ - Make a mutable copy of a pure array. Keep in mind that both `freeze` and `thaw` trigger a
 -- copy of the full array.
 --
@@ -204,7 +221,7 @@
 -- ==== __Example__
 --
 -- >>> import Data.Massiv.Array
--- >>> marr <- new @P @_ @Int (Sz2 2 6)
+-- >>> marr <- newMArray @P @_ @Int (Sz2 2 6) 0
 -- >>> forM_ (range Seq 0 (Ix2 1 4)) $ \ix -> write marr ix 9
 -- >>> freeze Seq marr
 -- Array P Seq (Sz (2 :. 6))
@@ -1015,7 +1032,7 @@
 -- >>> :set -XTypeApplications
 -- >>> import Control.Monad.ST
 -- >>> import Data.Massiv.Array
--- >>> runST $ new @P @Ix1 @Int (Sz1 3) >>= (\ma -> modifyM_ ma (pure . (+10)) 1 >> freezeS ma)
+-- >>> runST $ newMArray' @P @Ix1 @Int (Sz1 3) >>= (\ma -> modifyM_ ma (pure . (+10)) 1 >> freezeS ma)
 -- Array P Seq (Sz1 3)
 --   [ 0, 10, 0 ]
 --
diff --git a/src/Data/Massiv/Array/Mutable/Algorithms.hs b/src/Data/Massiv/Array/Mutable/Algorithms.hs
--- a/src/Data/Massiv/Array/Mutable/Algorithms.hs
+++ b/src/Data/Massiv/Array/Mutable/Algorithms.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Algorithms
--- Copyright   : (c) Alexey Kuleshevich 2019
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Mutable/Atomic.hs b/src/Data/Massiv/Array/Mutable/Atomic.hs
--- a/src/Data/Massiv/Array/Mutable/Atomic.hs
+++ b/src/Data/Massiv/Array/Mutable/Atomic.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE RankNTypes #-}
@@ -6,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Atomic
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Mutable/Internal.hs b/src/Data/Massiv/Array/Mutable/Internal.hs
--- a/src/Data/Massiv/Array/Mutable/Internal.hs
+++ b/src/Data/Massiv/Array/Mutable/Internal.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE ExplicitForAll #-}
 -- |
 -- Module      : Data.Massiv.Array.Mutable.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Numeric
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -45,20 +45,15 @@
   , (!><!)
   , multiplyMatrices
   , multiplyMatricesTransposed
-  -- Deprecated:
-  , (#>)
-  , (|*|)
-  , multiplyTransposed
   -- * Norms
   , normL2
-  -- ** Simple matrices
+  -- * Simple matrices
   , identityMatrix
   , lowerTriangular
   , upperTriangular
   , negateA
   , absA
   , signumA
-  , fromIntegerA
   -- * Integral
   , quotA
   , remA
@@ -73,9 +68,7 @@
   , (!/!)
   , (.^^)
   , recipA
-  , fromRationalA
   -- * Floating
-  , piA
   , expA
   , logA
   , sqrtA
@@ -107,9 +100,7 @@
 import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Array.Delayed.Push
 import Data.Massiv.Array.Manifest.Internal
-import Data.Massiv.Array.Ops.Fold as A
 import Data.Massiv.Array.Ops.Map as A
-import Data.Massiv.Array.Ops.Transform as A
 import Data.Massiv.Array.Ops.Construct
 import Data.Massiv.Core
 import Data.Massiv.Core.Common
@@ -125,21 +116,6 @@
 infixl 7  !*!, .*., .*, *., !/!, ./., ./, /., `quotA`, `remA`, `divA`, `modA`
 infixl 6  !+!, .+., .+, +., !-!, .-., .-, -.
 
-liftArray2Matching
-  :: (Source r1 ix a, Source r2 ix b)
-  => (a -> b -> e) -> Array r1 ix a -> Array r2 ix b -> Array D ix e
-liftArray2Matching f !arr1 !arr2
-  | sz1 == sz2 =
-    makeArray
-      (getComp arr1 <> getComp arr2)
-      sz1
-      (\ !ix -> f (unsafeIndex arr1 ix) (unsafeIndex arr2 ix))
-  | otherwise = throw $ SizeMismatchException (size arr1) (size arr2)
-  where
-    sz1 = size arr1
-    sz2 = size arr2
-{-# INLINE liftArray2Matching #-}
-
 liftArray2M ::
      (Load r ix e, Numeric r e, MonadThrow m)
   => (e -> e -> e)
@@ -338,26 +314,7 @@
 (.^) = powerPointwise
 {-# INLINE (.^) #-}
 
--- | Matrix-vector product
---
--- Inner dimensions must agree, otherwise `SizeMismatchException`
---
--- @since 0.5.2
-(#>) :: (MonadThrow m, Num e, Source (R r) Ix1 e, Manifest r' Ix1 e, OuterSlice r Ix2 e) =>
-        Array r Ix2 e -- ^ Matrix
-     -> Array r' Ix1 e -- ^ Vector
-     -> m (Array D Ix1 e)
-mm #> v
-  | mCols /= n = throwM $ SizeMismatchException (size mm) (Sz2 n 1)
-  | otherwise = pure $ makeArray (getComp mm <> getComp v) (Sz1 mRows) $ \i ->
-      A.foldlS (+) 0 (A.zipWith (*) (unsafeOuterSlice mm i) v)
-  where
-    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
@@ -372,7 +329,7 @@
 -- /__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 :: (FoldNumeric 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
@@ -383,7 +340,7 @@
 
 
 unsafeDotProductIO ::
-     (MonadUnliftIO m, Numeric r b, Source r ix b)
+     (MonadUnliftIO m, FoldNumeric r b, Source r ix b)
   => Array r ix b
   -> Array r ix b
   -> m b
@@ -410,14 +367,14 @@
 -- | 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 :: (Floating e, FoldNumeric 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)
+     (MonadUnliftIO m, FoldNumeric r b, Source r ix b)
   => Array r ix b
   -> Int
   -> m b
@@ -444,7 +401,7 @@
 --
 -- @since 0.5.6
 (.><) ::
-     (MonadThrow m, Numeric r e, Source r Ix1 e, Source r Ix2 e)
+     (MonadThrow m, FoldNumeric 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)
@@ -496,8 +453,8 @@
 (><.) :: (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 = delay <$> multiplyVectorByMatrix v mm
+      -> m (Vector r e)
+(><.) = multiplyVectorByMatrix
 {-# INLINE (><.) #-}
 
 -- | Multiply a row vector by a matrix. Same as `><.` but returns computed vector instead of
@@ -517,7 +474,7 @@
   | otherwise =
     pure $!
     unsafePerformIO $ do
-      mv <- new (Sz mCols)
+      mv <- newMArray (Sz mCols) 0
       withMassivScheduler_ comp $ \scheduler -> do
         let loopCols x ivto =
               fix $ \go im iv ->
@@ -547,7 +504,7 @@
      (Numeric r e, Mutable r Ix1 e, Mutable r Ix2 e)
   => Vector r e -- ^ Row vector (Used many times, so make sure it is computed)
   -> Matrix r e -- ^ Matrix
-  -> Vector D e
+  -> Vector r e
 (><!) v mm = throwEither (v ><. mm)
 {-# INLINE (><!) #-}
 
@@ -560,9 +517,9 @@
 -- ====__Examples__
 --
 -- >>> a1 = makeArrayR P Seq (Sz2 5 6) $ \(i :. j) -> i + j
--- >>> a2 = makeArrayR D Seq (Sz2 6 5) $ \(i :. j) -> i - j
+-- >>> a2 = makeArrayR P Seq (Sz2 6 5) $ \(i :. j) -> i - j
 -- >>> a1 !><! a2
--- Array D Seq (Sz (5 :. 5))
+-- Array P Seq (Sz (5 :. 5))
 --   [ [ 55, 40, 25, 10, -5 ]
 --   , [ 70, 49, 28, 7, -14 ]
 --   , [ 85, 58, 31, 4, -23 ]
@@ -571,7 +528,7 @@
 --   ]
 --
 -- @since 0.5.6
-(!><!) :: (Numeric r e, Mutable r Ix2 e, Source r' Ix2 e) => Matrix r e -> Matrix r' e -> Matrix D e
+(!><!) :: (Numeric r e, Mutable r Ix2 e) => Matrix r e -> Matrix r e -> Matrix r e
 (!><!) a1 a2 = throwEither (a1 `multiplyMatrices` a2)
 {-# INLINE (!><!) #-}
 
@@ -581,11 +538,7 @@
 -- /__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)
+(.><.) :: (Numeric r e, Mutable r Ix2 e, MonadThrow m) => Matrix r e -> Matrix r e -> m (Matrix r e)
 (.><.) = multiplyMatrices
 {-# INLINE (.><.) #-}
 
@@ -594,21 +547,145 @@
 --
 -- @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 #-}
+     (Numeric r e, Mutable r Ix2 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 $ setComp comp empty
+  | 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)
 
-{-# RULES
-"multiplyMatricesTransposed" [~1] forall arr1 arr2 . multiplyMatrices arr1 (transpose arr2) =
-    multiplyMatricesTransposedFused arr1 (convert arr2)
- #-}
+            | 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')
 --
@@ -618,7 +695,7 @@
 multiplyMatricesTransposed ::
      (Numeric r e, Manifest r Ix2 e, MonadThrow m)
   => Matrix r e
-  -> Matrix r  e
+  -> Matrix r e
   -> m (Matrix D e)
 multiplyMatricesTransposed arr1 arr2
   | n1 /= m2 = throwM $ SizeMismatchException (size arr1) (Sz2 m2 n2)
@@ -634,86 +711,8 @@
     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
-     , Source (R r) Ix1 e
-     , Num e
-     , MonadThrow m
-     )
-  => Array r Ix2 e
-  -> Array r Ix2 e
-  -> m (Array r Ix2 e)
-multiplyTransposedFused arr1 arr2 = compute <$> multiplyTransposed arr1 arr2
-{-# INLINE multiplyTransposedFused #-}
-
-
-multArrs :: forall r r' e m.
-            ( 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 D Ix2 e)
-multArrs arr1 arr2 = multiplyTransposed arr1 arr2'
-  where
-    arr2' :: Array r Ix2 e
-    arr2' = compute $ transpose arr2
-{-# INLINE multArrs #-}
-
--- | Computes the matrix-matrix transposed product (i.e. A * A')
-multiplyTransposed ::
-     ( Manifest 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 D Ix2 e)
-multiplyTransposed arr1 arr2
-  | n1 /= m2 = throwM $ SizeMismatchException (size arr1) (size arr2)
-  | otherwise =
-    pure $
-    DArray (getComp arr1 <> getComp arr2) (SafeSz (m1 :. n2)) $ \(i :. j) ->
-      A.foldlS (+) 0 (A.zipWith (*) (unsafeOuterSlice arr1 i) (unsafeOuterSlice arr2 j))
-  where
-    SafeSz (m1 :. n1) = size arr1
-    SafeSz (n2 :. m2) = size arr2
-{-# INLINE multiplyTransposed #-}
-
 -- | Create an indentity matrix.
 --
 -- ==== __Example__
@@ -807,12 +806,6 @@
 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.
 --
@@ -898,21 +891,6 @@
 {-# INLINE recipA #-}
 
 
-fromRationalA
-  :: (Index ix, Fractional e)
-  => 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
@@ -1162,7 +1140,7 @@
 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)
-quotRemA arr1 = A.unzip . liftArray2Matching (quotRem) arr1
+quotRemA arr1 = A.unzip . liftArray2Matching quotRem arr1
 {-# INLINE quotRemA #-}
 
 
@@ -1177,7 +1155,7 @@
 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)
-divModA arr1 = A.unzip . liftArray2Matching (divMod) arr1
+divModA arr1 = A.unzip . liftArray2Matching divMod arr1
 {-# INLINE divModA #-}
 
 
diff --git a/src/Data/Massiv/Array/Numeric/Integral.hs b/src/Data/Massiv/Array/Numeric/Integral.hs
--- a/src/Data/Massiv/Array/Numeric/Integral.hs
+++ b/src/Data/Massiv/Array/Numeric/Integral.hs
@@ -1,8 +1,8 @@
-{-# LANGUAGE BangPatterns     #-}
+{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Numeric.Integral
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -32,14 +32,15 @@
   , fromFunctionMidpoint
   ) where
 
-import           Data.Coerce
-import           Data.Massiv.Array.Delayed.Pull      (D)
-import           Data.Massiv.Array.Delayed.Windowed  (DW)
-import           Data.Massiv.Array.Manifest.Internal
-import           Data.Massiv.Array.Ops.Construct     (rangeInclusive)
-import           Data.Massiv.Array.Ops.Transform     (extract')
-import           Data.Massiv.Array.Stencil
-import           Data.Massiv.Core.Common
+import Data.Coerce
+import Data.Massiv.Array.Delayed.Pull (D)
+import Data.Massiv.Array.Delayed.Windowed (DW)
+import Data.Massiv.Array.Manifest.Internal
+import Data.Massiv.Array.Ops.Construct (rangeInclusive)
+import Data.Massiv.Array.Ops.Transform (extract')
+import Data.Massiv.Array.Stencil
+import Data.Massiv.Array.Unsafe
+import Data.Massiv.Core.Common
 
 
 -- |
@@ -56,8 +57,8 @@
   -> Int -- ^ @n@ - number of sample points.
   -> Stencil ix e e
 midpointStencil dx dim k =
-  makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim k)) zeroIndex $ \g ->
-    pure dx * loop 0 (< k) (+ 1) 0 (\i -> (+ g (setDim' zeroIndex dim i)))
+  makeUnsafeStencil (Sz (setDim' (pureIndex 1) dim k)) zeroIndex $ \_ g ->
+    dx * loop 0 (< k) (+ 1) 0 (\i -> (+ g (setDim' zeroIndex dim i)))
 {-# INLINE midpointStencil #-}
 
 
@@ -75,8 +76,8 @@
   -> Int -- ^ @n@ - number of sample points.
   -> Stencil ix e e
 trapezoidStencil dx dim n =
-  makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \g ->
-    pure dx / 2 *
+  makeUnsafeStencil (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \_ g ->
+    dx / 2 *
     (loop 1 (< n) (+ 1) (g zeroIndex) (\i -> (+ 2 * g (setDim' zeroIndex dim i))) +
      g (setDim' zeroIndex dim n))
 {-# INLINE trapezoidStencil #-}
@@ -100,12 +101,12 @@
     error $
     "Number of sample points for Simpson's rule stencil should be even, but received: " ++ show n
   | otherwise =
-    makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \g ->
+    makeUnsafeStencil (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \_ g ->
       let simAcc i (prev, acc) =
             let !fx3 = g (setDim' zeroIndex dim (i + 2))
                 !newAcc = acc + prev + 4 * g (setDim' zeroIndex dim (i + 1)) + fx3
              in (fx3, newAcc)
-       in pure dx / 3 * snd (loop 2 (< n - 1) (+ 2) (simAcc 0 (g zeroIndex, 0)) simAcc)
+       in dx / 3 * snd (loop 2 (< n - 1) (+ 2) (simAcc 0 (g zeroIndex, 0)) simAcc)
 {-# INLINE simpsonsStencil #-}
 
 -- | Integrate with a stencil along a particular dimension.
diff --git a/src/Data/Massiv/Array/Ops/Construct.hs b/src/Data/Massiv/Array/Ops/Construct.hs
--- a/src/Data/Massiv/Array/Ops/Construct.hs
+++ b/src/Data/Massiv/Array/Ops/Construct.hs
@@ -1,15 +1,13 @@
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE DataKinds #-}
-{-# LANGUAGE ExplicitForAll #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeOperators #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Construct
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -113,14 +111,6 @@
 makeVectorR :: Construct r Ix1 e => r -> Comp -> Sz1 -> (Ix1 -> e) -> Array r Ix1 e
 makeVectorR _ = makeArray
 {-# INLINE makeVectorR #-}
-
-
--- | Replicate the same element
---
--- @since 0.3.0
-replicate :: forall ix e . Index ix => Comp -> Sz ix -> e -> Array DL ix e
-replicate comp sz e = makeLoadArray comp sz e $ \_ _ -> pure ()
-{-# INLINE replicate #-}
 
 
 newtype STA r ix a = STA {_runSTA :: forall s. MArray s r ix a -> ST s (Array r ix a)}
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Fold
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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
@@ -5,7 +5,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Fold.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Ops/Map.hs b/src/Data/Massiv/Array/Ops/Map.hs
--- a/src/Data/Massiv/Array/Ops/Map.hs
+++ b/src/Data/Massiv/Array/Ops/Map.hs
@@ -4,7 +4,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Map
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Ops/Slice.hs b/src/Data/Massiv/Array/Ops/Slice.hs
--- a/src/Data/Massiv/Array/Ops/Slice.hs
+++ b/src/Data/Massiv/Array/Ops/Slice.hs
@@ -3,7 +3,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Slice
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Ops/Sort.hs b/src/Data/Massiv/Array/Ops/Sort.hs
--- a/src/Data/Massiv/Array/Ops/Sort.hs
+++ b/src/Data/Massiv/Array/Ops/Sort.hs
@@ -3,7 +3,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Sort
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Ops/Transform.hs b/src/Data/Massiv/Array/Ops/Transform.hs
--- a/src/Data/Massiv/Array/Ops/Transform.hs
+++ b/src/Data/Massiv/Array/Ops/Transform.hs
@@ -6,7 +6,7 @@
 {-# OPTIONS_GHC -fno-warn-redundant-constraints #-}
 -- |
 -- Module      : Data.Massiv.Array.Ops.Transform
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Array/Stencil.hs b/src/Data/Massiv/Array/Stencil.hs
--- a/src/Data/Massiv/Array/Stencil.hs
+++ b/src/Data/Massiv/Array/Stencil.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -14,7 +14,6 @@
 module Data.Massiv.Array.Stencil
   ( -- * Stencil
     Stencil
-  , Value
   , makeStencil
   , makeStencilDef
   , getStencilSize
@@ -42,12 +41,9 @@
   , rmapStencil
   -- * Convolution
   , module Data.Massiv.Array.Stencil.Convolution
-  -- * Re-export
-  , Default(def)
   ) where
 
 import Data.Coerce
-import Data.Default.Class (Default(def))
 import Data.Massiv.Array.Delayed.Windowed
 import Data.Massiv.Array.Manifest
 import Data.Massiv.Array.Stencil.Convolution
@@ -175,6 +171,7 @@
     , paddingWithElement = border
     }
 
+
 -- | Apply a constructed stencil over an array. Resulting array must be
 -- `Data.Massiv.Array.compute`d in order to be useful. Unlike `mapStencil`, the size of
 -- the resulting array will not necesserally be the same as the source array, which will
@@ -197,7 +194,7 @@
       DArray
         (getComp arr)
         sz
-        (unValue . stencilF (Value . borderIndex border arr) . liftIndex2 (+) offset)
+        (stencilF (borderIndex border arr) (borderIndex border arr) . liftIndex2 (+) offset)
     -- Size by which the resulting array will shrink (not accounting for padding)
     !shrinkSz = Sz (liftIndex (subtract 1) (unSz sSz))
     !sz = liftSz2 (-) (SafeSz (liftIndex2 (+) po (liftIndex2 (+) pb (unSz (size arr))))) shrinkSz
@@ -206,7 +203,7 @@
       Window
         { windowStart = po
         , windowSize = wsz
-        , windowIndex = unValue . stencilF (Value . unsafeIndex arr) . liftIndex2 (+) offset
+        , windowIndex = stencilF (unsafeIndex arr) (index' arr) . liftIndex2 (+) offset
         , windowUnrollIx2 = unSz . fst <$> pullOutSzM sSz 2
         }
 {-# INLINE applyStencil #-}
@@ -218,6 +215,9 @@
 -- outside the stencil box will result in a runtime error upon stencil
 -- creation.
 --
+-- /Note/ - Once correctness of stencil is verified then switching to `makeUnsafeStencil`
+-- is recommended in order to get the most performance out of the `Stencil`
+--
 -- ==== __Example__
 --
 -- Below is an example of creating a `Stencil`, which, when mapped over a
@@ -226,7 +226,7 @@
 --
 -- /Note/ - Make sure to add an @INLINE@ pragma, otherwise performance will be terrible.
 --
--- > average3x3Stencil :: (Default a, Fractional a) => Stencil Ix2 a a
+-- > average3x3Stencil :: Fractional a => Stencil Ix2 a a
 -- > average3x3Stencil = makeStencil (Sz (3 :. 3)) (1 :. 1) $ \ get ->
 -- >   (  get (-1 :. -1) + get (-1 :. 0) + get (-1 :. 1) +
 -- >      get ( 0 :. -1) + get ( 0 :. 0) + get ( 0 :. 1) +
@@ -235,10 +235,10 @@
 --
 -- @since 0.1.0
 makeStencil
-  :: (Index ix, Default e)
+  :: Index ix
   => Sz ix -- ^ Size of the stencil
   -> ix -- ^ Center of the stencil
-  -> ((ix -> Value e) -> Value a)
+  -> ((ix -> e) -> a)
   -- ^ Stencil function that receives a "get" function as it's argument that can
   -- retrieve values of cells in the source array with respect to the center of
   -- the stencil. Stencil function must return a value that will be assigned to
@@ -246,7 +246,11 @@
   -- cannot go outside the boundaries of the stencil, otherwise an error will be
   -- raised during stencil creation.
   -> Stencil ix e a
-makeStencil = makeStencilDef def
+makeStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
+  where
+    stencil _ getVal !ix =
+      inline relStencil $ \ !ixD -> getVal (liftIndex2 (+) ix ixD)
+    {-# INLINE stencil #-}
 {-# INLINE makeStencil #-}
 
 -- | Same as `makeStencil`, but with ability to specify default value for stencil validation.
@@ -257,16 +261,17 @@
   => e -- ^ Default element that will be used for stencil validation only.
   -> Sz ix -- ^ Size of the stencil
   -> ix -- ^ Center of the stencil
-  -> ((ix -> Value e) -> Value a)
+  -> ((ix -> e) -> a)
   -- ^ Stencil function.
   -> Stencil ix e a
-makeStencilDef defVal !sSz !sCenter relStencil =
-  validateStencil defVal $ Stencil sSz sCenter stencil
+makeStencilDef _defVal !sSz !sCenter relStencil =
+  Stencil sSz sCenter stencil
   where
-    stencil getVal !ix =
+    stencil _ getVal !ix =
       inline relStencil $ \ !ixD -> getVal (liftIndex2 (+) ix ixD)
     {-# INLINE stencil #-}
 {-# INLINE makeStencilDef #-}
+{-# DEPRECATED makeStencilDef "In favor of `makeStencil`. Validation is no longer possible" #-}
 
 -- | Identity stencil that does not change the elements of the source array.
 --
diff --git a/src/Data/Massiv/Array/Stencil/Convolution.hs b/src/Data/Massiv/Array/Stencil/Convolution.hs
--- a/src/Data/Massiv/Array/Stencil/Convolution.hs
+++ b/src/Data/Massiv/Array/Stencil/Convolution.hs
@@ -2,7 +2,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Convolution
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -23,6 +23,10 @@
 -- | Create a convolution stencil by specifying border resolution technique and
 -- an accumulator function.
 --
+-- /Note/ - Using `Data.Massiv.Array.Stencil.Unsafe.makeUnsafeConvolutionStencil` will be
+-- much faster, therefore it is recommended to switch from this function, after manual
+-- verification that the created stencil behaves as expected.
+--
 -- ==== __Examples__
 --
 -- Here is how to create a 2D horizontal Sobel Stencil:
@@ -39,13 +43,13 @@
   :: (Index ix, Num e)
   => Sz ix
   -> ix
-  -> ((ix -> Value e -> Value e -> Value e) -> Value e -> Value e)
+  -> ((ix -> e -> e -> e) -> e -> e)
   -> Stencil ix e e
 makeConvolutionStencil !sz !sCenter relStencil =
-  validateStencil 0 $ Stencil sz sInvertCenter stencil
+  Stencil sz sInvertCenter stencil
   where
     !sInvertCenter = liftIndex2 (-) (liftIndex (subtract 1) (unSz sz)) sCenter
-    stencil getVal !ix =
+    stencil _ getVal !ix =
       (inline relStencil $ \ !ixD !kVal !acc -> getVal (liftIndex2 (-) ix ixD) * kVal + acc) 0
     {-# INLINE stencil #-}
 {-# INLINE makeConvolutionStencil #-}
@@ -66,10 +70,9 @@
     !szi1 = liftIndex (subtract 1) szi
     !sInvertCenter = liftIndex2 (-) szi1 sCenter
     !sCenter = liftIndex (`quot` 2) szi
-    stencil getVal !ix = Value (ifoldlS accum 0 kArr) where
+    stencil uget _ !ix = ifoldlS accum 0 kArr where
       !ixOff = liftIndex2 (+) ix sCenter
-      accum !acc !kIx !kVal =
-        unValue (getVal (liftIndex2 (-) ixOff kIx)) * kVal + acc
+      accum !acc !kIx !kVal = uget (liftIndex2 (-) ixOff kIx) * kVal + acc
       {-# INLINE accum #-}
     {-# INLINE stencil #-}
 {-# INLINE makeConvolutionStencilFromKernel #-}
@@ -77,16 +80,20 @@
 
 -- | Make a <https://en.wikipedia.org/wiki/Cross-correlation cross-correlation> stencil
 --
+-- /Note/ - Using `Data.Massiv.Array.Stencil.Unsafe.makeUnsafeCorrelationStencil` will be
+-- much faster, therefore it is recommended to switch from this function, after manual
+-- verification that the created stencil behaves as expected.
+--
 -- @since 0.1.5
 makeCorrelationStencil
   :: (Index ix, Num e)
   => Sz ix
   -> ix
-  -> ((ix -> Value e -> Value e -> Value e) -> Value e -> Value e)
+  -> ((ix -> e -> e -> e) -> e -> e)
   -> Stencil ix e e
-makeCorrelationStencil !sSz !sCenter relStencil = validateStencil 0 $ Stencil sSz sCenter stencil
+makeCorrelationStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
   where
-    stencil getVal !ix =
+    stencil _ getVal !ix =
       (inline relStencil $ \ !ixD !kVal !acc -> getVal (liftIndex2 (+) ix ixD) * kVal + acc) 0
     {-# INLINE stencil #-}
 {-# INLINE makeCorrelationStencil #-}
@@ -104,10 +111,9 @@
   where
     !sz = size kArr
     !sCenter = liftIndex (`div` 2) $ unSz sz
-    stencil getVal !ix = Value (ifoldlS accum 0 kArr) where
+    stencil uget _ !ix = ifoldlS accum 0 kArr where
       !ixOff = liftIndex2 (-) ix sCenter
-      accum !acc !kIx !kVal =
-        unValue (getVal (liftIndex2 (+) ixOff kIx)) * kVal + acc
+      accum !acc !kIx !kVal = uget (liftIndex2 (+) ixOff kIx) * kVal + acc
       {-# INLINE accum #-}
     {-# INLINE stencil #-}
 {-# INLINE makeCorrelationStencilFromKernel #-}
diff --git a/src/Data/Massiv/Array/Stencil/Internal.hs b/src/Data/Massiv/Array/Stencil/Internal.hs
--- a/src/Data/Massiv/Array/Stencil/Internal.hs
+++ b/src/Data/Massiv/Array/Stencil/Internal.hs
@@ -2,11 +2,10 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -14,123 +13,29 @@
 --
 module Data.Massiv.Array.Stencil.Internal
   ( Stencil(..)
-  , Value(..)
   , dimapStencil
   , lmapStencil
   , rmapStencil
-  , validateStencil
   ) where
 
 import Control.Applicative
 import Control.DeepSeq
-import Data.Massiv.Array.Delayed.Pull
 import Data.Massiv.Core.Common
 
--- | Stencil is abstract description of how to handle elements in the neighborhood of every array
--- cell in order to compute a value for the cells in the new array. Use `Data.Array.makeStencil` and
--- `Data.Array.makeConvolutionStencil` in order to create a stencil.
+-- | Stencil is abstract description of how to handle elements in the neighborhood of
+-- every array cell in order to compute a value for the cells in the new array. Use
+-- `Data.Massiv.Array.makeStencil` and `Data.Massiv.Array.makeConvolutionStencil` in order
+-- to create a stencil.
 data Stencil ix e a = Stencil
   { stencilSize   :: !(Sz ix)
   , stencilCenter :: !ix
-  , stencilFunc   :: (ix -> Value e) -> ix -> Value a
+  , stencilFunc   :: (ix -> e) -> (ix -> e) -> ix -> a
   }
 
 
 instance Index ix => NFData (Stencil ix e a) where
   rnf (Stencil sz ix f) = sz `deepseq` ix `deepseq` f `seq` ()
 
--- | This is a simple wrapper for value of an array cell. It is used in order to improve safety of
--- `Stencil` mapping. Using various class instances, such as `Num` and `Functor` for example, make
--- it possible to manipulate the value, without having direct access to it.
-newtype Value e = Value
-  { unValue :: e
-  } deriving (Bounded)
-
-instance Functor Value where
-  fmap f (Value e) = Value (f e)
-  {-# INLINE fmap #-}
-
-instance Applicative Value where
-  pure = Value
-  {-# INLINE pure #-}
-  (<*>) (Value f) (Value e) = Value (f e)
-  {-# INLINE (<*>) #-}
-
--- | @since 0.1.5
-instance Semigroup a => Semigroup (Value a) where
-  Value a <> Value b = Value (a <> b)
-  {-# INLINE (<>) #-}
-
--- | @since 0.1.5
-instance Monoid a => Monoid (Value a) where
-  mempty = Value mempty
-  {-# INLINE mempty #-}
-  Value a `mappend` Value b = Value (a `mappend` b)
-  {-# INLINE mappend #-}
-
-instance Num e => Num (Value e) where
-  (+) = liftA2 (+)
-  {-# INLINE (+) #-}
-  (*) = liftA2 (*)
-  {-# INLINE (*) #-}
-  negate = fmap negate
-  {-# INLINE negate #-}
-  abs = fmap abs
-  {-# INLINE abs #-}
-  signum = fmap signum
-  {-# INLINE signum #-}
-  fromInteger = Value . fromInteger
-  {-# INLINE fromInteger #-}
-
-instance Fractional e => Fractional (Value e) where
-  (/) = liftA2 (/)
-  {-# INLINE (/) #-}
-  recip = fmap recip
-  {-# INLINE recip #-}
-  fromRational = pure . fromRational
-  {-# INLINE fromRational #-}
-
-instance Floating e => Floating (Value e) where
-  pi = pure pi
-  {-# INLINE pi #-}
-  exp = fmap exp
-  {-# INLINE exp #-}
-  log = fmap log
-  {-# INLINE log #-}
-  sqrt = fmap sqrt
-  {-# INLINE sqrt #-}
-  (**) = liftA2 (**)
-  {-# INLINE (**) #-}
-  logBase = liftA2 logBase
-  {-# INLINE logBase #-}
-  sin = fmap sin
-  {-# INLINE sin #-}
-  cos = fmap cos
-  {-# INLINE cos #-}
-  tan = fmap tan
-  {-# INLINE tan #-}
-  asin = fmap asin
-  {-# INLINE asin #-}
-  acos = fmap acos
-  {-# INLINE acos #-}
-  atan = fmap atan
-  {-# INLINE atan #-}
-  sinh = fmap sinh
-  {-# INLINE sinh #-}
-  cosh = fmap cosh
-  {-# INLINE cosh #-}
-  tanh = fmap tanh
-  {-# INLINE tanh #-}
-  asinh = fmap asinh
-  {-# INLINE asinh #-}
-  acosh = fmap acosh
-  {-# INLINE acosh #-}
-  atanh = fmap atanh
-  {-# INLINE atanh #-}
-
-
-
-
 instance Functor (Stencil ix e) where
   fmap = rmapStencil
   {-# INLINE fmap #-}
@@ -144,7 +49,7 @@
 dimapStencil :: (c -> d) -> (a -> b) -> Stencil ix d a -> Stencil ix c b
 dimapStencil f g stencil@Stencil {stencilFunc = sf} = stencil {stencilFunc = sf'}
   where
-    sf' s = Value . g . unValue . sf (Value . f . unValue . s)
+    sf' us s = g . sf (f . us) (f . s)
     {-# INLINE sf' #-}
 {-# INLINE dimapStencil #-}
 
@@ -159,11 +64,12 @@
 lmapStencil :: (c -> d) -> Stencil ix d a -> Stencil ix c a
 lmapStencil f stencil@Stencil {stencilFunc = sf} = stencil {stencilFunc = sf'}
   where
-    sf' s = sf (Value . f . unValue . s)
+    sf' us s = sf (f . us) (f . s)
     {-# INLINE sf' #-}
 {-# INLINE lmapStencil #-}
 
--- | A covariant map over the right most type argument. In other words a usual Functor `fmap`:
+-- | A covariant map over the right most type argument. In other words the usual `fmap`
+-- from `Functor`:
 --
 -- > fmap == rmapStencil
 --
@@ -171,34 +77,39 @@
 rmapStencil :: (a -> b) -> Stencil ix e a -> Stencil ix e b
 rmapStencil f stencil@Stencil {stencilFunc = sf} = stencil {stencilFunc = sf'}
   where
-    sf' s = Value . f . unValue . sf s
+    sf' us s = f . sf us s
     {-# INLINE sf' #-}
 {-# INLINE rmapStencil #-}
 
+unionStencilCenters :: Index ix => Stencil ix e1 a1 -> Stencil ix e2 a2 -> ix
+unionStencilCenters (Stencil _ sC1 _) (Stencil _ sC2 _) = liftIndex2 max sC1 sC2
+{-# INLINE unionStencilCenters #-}
 
+unionStencilSizes :: Index ix => ix -> Stencil ix e1 a1 -> Stencil ix e2 a2 -> Sz ix
+unionStencilSizes maxCenter (Stencil (SafeSz sSz1) sC1 _) (Stencil (SafeSz sSz2) sC2 _) =
+  Sz $ liftIndex2 (+) maxCenter $ liftIndex2 max (liftIndex2 (-) sSz1 sC1) (liftIndex2 (-) sSz2 sC2)
+{-# INLINE unionStencilSizes #-}
 
--- TODO: Figure out interchange law (u <*> pure y = pure ($ y) <*> u) and issue
--- with discarding size and center. Best idea so far is to increase stencil size to
--- the maximum one and shift the center of the other stencil so that they both match
--- up. This approach would also remove requirement to validate the result
--- Stencil - both stencils are trusted, increasing the size will not affect the
--- safety.
+-- TODO: Test interchange law (u <*> pure y = pure ($ y) <*> u)
 instance Index ix => Applicative (Stencil ix e) where
-  pure a = Stencil oneSz zeroIndex (const (const (Value a)))
+  pure a = Stencil oneSz zeroIndex (\_ _ _ -> a)
   {-# INLINE pure #-}
-  (<*>) (Stencil (SafeSz sSz1) sC1 f1) (Stencil (SafeSz sSz2) sC2 f2) = Stencil newSz maxCenter stF
+  (<*>) s1@(Stencil _ _ f1) s2@(Stencil _ _ f2) = Stencil newSz maxCenter stF
     where
-      stF gV !ix = Value (unValue (f1 gV ix) (unValue (f2 gV ix)))
+      stF ug gV !ix = f1 ug gV ix (f2 ug gV ix)
       {-# INLINE stF #-}
-      !newSz =
-        Sz
-          (liftIndex2
-             (+)
-             maxCenter
-             (liftIndex2 max (liftIndex2 (-) sSz1 sC1) (liftIndex2 (-) sSz2 sC2)))
-      !maxCenter = liftIndex2 max sC1 sC2
+      !newSz = unionStencilSizes maxCenter s1 s2
+      !maxCenter = unionStencilCenters s1 s2
   {-# INLINE (<*>) #-}
 
+  liftA2 f s1@(Stencil _ _ f1) s2@(Stencil _ _ f2) = Stencil newSz maxCenter stF
+    where
+      stF ug gV !ix = f (f1 ug gV ix) (f2 ug gV ix)
+      {-# INLINE stF #-}
+      !newSz = unionStencilSizes maxCenter s1 s2
+      !maxCenter = unionStencilCenters s1 s2
+  {-# INLINE liftA2 #-}
+
 instance (Index ix, Num a) => Num (Stencil ix e a) where
   (+) = liftA2 (+)
   {-# INLINE (+) #-}
@@ -260,21 +171,3 @@
   {-# INLINE acosh #-}
   atanh = fmap atanh
   {-# INLINE atanh #-}
-
-
-safeStencilIndex :: Index ix => Array D ix e -> ix -> e
-safeStencilIndex DArray {..} ix
-  | isSafeIndex dSize ix = dIndex ix
-  | otherwise = throw $ IndexOutOfBoundsException dSize ix
-
-
--- | Make sure constructed stencil doesn't index outside the allowed stencil size boundary.
-validateStencil
-  :: Index ix
-  => e -> Stencil ix e a -> Stencil ix e a
-validateStencil d s@(Stencil sSz sCenter stencil)
-  | isSafeIndex sSz sCenter =
-    let valArr = DArray Seq sSz (const d)
-     in stencil (Value . safeStencilIndex valArr) sCenter `seq` s
-  | otherwise = throw $ IndexOutOfBoundsException sSz sCenter
-{-# INLINE validateStencil #-}
diff --git a/src/Data/Massiv/Array/Stencil/Unsafe.hs b/src/Data/Massiv/Array/Stencil/Unsafe.hs
--- a/src/Data/Massiv/Array/Stencil/Unsafe.hs
+++ b/src/Data/Massiv/Array/Stencil/Unsafe.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE RecordWildCards #-}
 -- |
 -- Module      : Data.Massiv.Array.Stencil.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -14,11 +14,11 @@
 module Data.Massiv.Array.Stencil.Unsafe
   ( -- * Stencil
     makeUnsafeStencil
+  , makeUnsafeConvolutionStencil
+  , makeUnsafeCorrelationStencil
   , unsafeTransformStencil
-  , unsafeMapStencil
   -- ** Deprecated
-  , mapStencilUnsafe
-  , forStencilUnsafe
+  , unsafeMapStencil
   ) where
 
 import Data.Massiv.Array.Delayed.Windowed (Array(..), DW, Window(..),
@@ -28,58 +28,8 @@
 import GHC.Exts (inline)
 
 
--- | Just as `mapStencilUnsafe` this is an unsafe version of the stencil
--- mapping. Arguments are in slightly different order and the indexing function returns
--- `Nothing` for elements outside the border.
---
--- @since 0.1.7
-forStencilUnsafe ::
-     (Source r ix e, Manifest r ix e)
-  => Array r ix e
-  -> Sz ix -- ^ Size of the stencil
-  -> ix -- ^ Center of the stencil
-  -> ((ix -> Maybe e) -> a)
-  -- ^ Stencil function that receives a "get" function as it's argument that can
-  -- retrieve values of cells in the source array with respect to the center of
-  -- the stencil. Stencil function must return a value that will be assigned to
-  -- the cell in the result array. Offset supplied to the "get" function
-  -- cannot go outside the boundaries of the stencil.
-  -> Array DW ix a
-forStencilUnsafe !arr !sSz !sCenter relStencil =
-  insertWindow (DArray (getComp arr) sz (stencil (index arr))) window
-  where
-    !window =
-      Window
-        { windowStart = sCenter
-        , windowSize = windowSz
-        , windowIndex = stencil (Just . unsafeIndex arr)
-        , windowUnrollIx2 = unSz . fst <$> pullOutSzM windowSz 2
-        }
-    !sz = size arr
-    !windowSz = Sz (liftIndex2 (-) (unSz sz) (liftIndex (subtract 1) (unSz sSz)))
-    stencil getVal !ix = inline relStencil $ \ !ixD -> getVal (liftIndex2 (+) ix ixD)
-    {-# INLINE stencil #-}
-{-# INLINE forStencilUnsafe #-}
-{-# DEPRECATED forStencilUnsafe "In favor of `unsafeMapStencil`" #-}
-
-
-mapStencilUnsafe ::
-     Manifest r ix e
-  => Border e
-  -> Sz ix
-  -> ix
-  -> ((ix -> e) -> a)
-  -> Array r ix e
-  -> Array DW ix a
-mapStencilUnsafe b sz ix f = unsafeMapStencil b sz ix (const f)
-{-# INLINE mapStencilUnsafe #-}
-{-# DEPRECATED mapStencilUnsafe "In favor of `unsafeMapStencil`" #-}
-
--- | This is an unsafe version of `Data.Massiv.Array.Stencil.mapStencil`, that does no
--- take `Stencil` as argument, as such it does no stencil validation. There is no
--- performance difference between the two, but the unsafe version has an advantage of not
--- requiring to deal with `Value` wrapper and has access to the actual index with the
--- array.
+-- | This is an unsafe version of `Data.Massiv.Array.Stencil.mapStencil`, which does not
+-- take a `Stencil`, but instead accepts all necessary information as separate arguments.
 --
 -- @since 0.5.0
 unsafeMapStencil ::
@@ -105,6 +55,7 @@
     stencil getVal !ix = inline (stencilF ix) $ \ !ixD -> getVal (liftIndex2 (+) ix ixD)
     {-# INLINE stencil #-}
 {-# INLINE unsafeMapStencil #-}
+{-# DEPRECATED unsafeMapStencil "In favor of `Data.Massiv.Array.mapStencil` that is applied to stencil created with `makeUnsafeStencil`" #-}
 
 
 -- | Similar to `Data.Massiv.Array.Stencil.makeStencil`, but there are no guarantees that the
@@ -121,12 +72,49 @@
   -> Stencil ix e a
 makeUnsafeStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
   where
-    stencil getVal !ix =
-      Value $ inline $ relStencil ix (unValue . getVal . liftIndex2 (+) ix)
+    stencil unsafeGetVal _getVal !ix =
+      inline $ relStencil ix (unsafeGetVal . liftIndex2 (+) ix)
     {-# INLINE stencil #-}
 {-# INLINE makeUnsafeStencil #-}
 
+-- | Same as `Data.Massiv.Array.Stencil.makeConvolutionStencil`, but will result in
+-- reading memory out of bounds and potential segfaults if supplied arguments are not valid.
+--
+-- @since 0.6.0
+makeUnsafeConvolutionStencil
+  :: (Index ix, Num e)
+  => Sz ix
+  -> ix
+  -> ((ix -> e -> e -> e) -> e -> e)
+  -> Stencil ix e e
+makeUnsafeConvolutionStencil !sz !sCenter relStencil =
+  Stencil sz sInvertCenter stencil
+  where
+    !sInvertCenter = liftIndex2 (-) (liftIndex (subtract 1) (unSz sz)) sCenter
+    stencil uget _ !ix =
+      (inline relStencil $ \ !ixD !kVal !acc -> uget (liftIndex2 (-) ix ixD) * kVal + acc) 0
+    {-# INLINE stencil #-}
+{-# INLINE makeUnsafeConvolutionStencil #-}
 
+-- | Same as `Data.Massiv.Array.Stencil.makeCorrelationStencil`, but will result in
+-- reading memory out of bounds and potential segfaults if supplied arguments are not
+-- valid.
+--
+-- @since 0.6.0
+makeUnsafeCorrelationStencil
+  :: (Index ix, Num e)
+  => Sz ix
+  -> ix
+  -> ((ix -> e -> e -> e) -> e -> e)
+  -> Stencil ix e e
+makeUnsafeCorrelationStencil !sSz !sCenter relStencil = Stencil sSz sCenter stencil
+  where
+    stencil _ getVal !ix =
+      (inline relStencil $ \ !ixD !kVal !acc -> getVal (liftIndex2 (+) ix ixD) * kVal + acc) 0
+    {-# INLINE stencil #-}
+{-# INLINE makeUnsafeCorrelationStencil #-}
+
+
 -- | Perform an arbitrary transformation of a stencil. This stencil modifier can be used for
 -- example to turn a vector stencil into a matrix stencil implement, or transpose a matrix
 -- stencil. It is really easy to get this wrong, so be extremely careful.
@@ -144,14 +132,14 @@
 --   , [ 4, 5, 6 ]
 --   , [ 7, 8, 9 ]
 --   ]
--- >>> let rowStencil = unsafeTransformStencil (\(Sz n) -> Sz (1 :. n)) (0 :.) $ \ f getVal (i :. j) -> f (getVal . (i :.)) j
+-- >>> let rowStencil = unsafeTransformStencil (\(Sz n) -> Sz (1 :. n)) (0 :.) $ \ f uget getVal (i :. j) -> f (uget  . (i :.)) (getVal . (i :.)) j
 -- >>> applyStencil noPadding (rowStencil (sumStencil (Sz1 3))) arr
 -- Array DW Seq (Sz (3 :. 1))
 --   [ [ 6 ]
 --   , [ 15 ]
 --   , [ 24 ]
 --   ]
--- >>> let columnStencil = unsafeTransformStencil (\(Sz n) -> Sz (n :. 1)) (:. 0) $ \ f getVal (i :. j) -> f (getVal . (:. j)) i
+-- >>> let columnStencil = unsafeTransformStencil (\(Sz n) -> Sz (n :. 1)) (:. 0) $ \ f uget getVal (i :. j) -> f (uget . (:. j)) (getVal . (:. j)) i
 -- >>> applyStencil noPadding (columnStencil (sumStencil (Sz1 3))) arr
 -- Array DW Seq (Sz (1 :. 3))
 --   [ [ 12, 15, 18 ]
@@ -163,7 +151,8 @@
   -- ^ Forward modifier for the size
   -> (ix' -> ix)
   -- ^ Forward index modifier
-  -> (((ix' -> Value e) -> ix' -> Value a) -> (ix -> Value e) -> ix -> Value a)
+  -> (((ix' -> e) -> (ix' -> e) -> ix' -> a)
+      -> (ix -> e) -> (ix -> e) -> ix -> a)
   -- ^ Inverse stencil function modifier
   -> Stencil ix' e a
   -- ^ Original stencil.
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
@@ -4,7 +4,7 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 -- |
 -- Module      : Data.Massiv.Array.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -80,9 +80,8 @@
   , unsafeAtomicXorIntArray
   , unsafeCasIntArray
     -- ** Other operations
-  , unsafeBoxedArray
-  , unsafeNormalBoxedArray
-  , unsafeFromBoxedVector
+  , coerceBoxedArray
+  , coerceNormalBoxedArray
   , unsafeUnstablePartitionRegionM
   , module Data.Massiv.Vector.Unsafe
   , module Data.Massiv.Array.Stencil.Unsafe
diff --git a/src/Data/Massiv/Core.hs b/src/Data/Massiv/Core.hs
--- a/src/Data/Massiv/Core.hs
+++ b/src/Data/Massiv/Core.hs
@@ -1,6 +1,6 @@
 -- |
 -- Module      : Data.Massiv.Core
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -14,7 +14,7 @@
   , MMatrix
   , Elt
   , Construct
-  , Load(R, loadArrayM, loadArrayWithSetM, defaultElement)
+  , Load(R, loadArrayM, loadArrayWithSetM)
   , Stream(..)
   , Source
   , Resize
@@ -34,6 +34,7 @@
   , Scheduler
   , SchedulerWS
   , Comp(Seq, Par, Par', ParOn, ParN)
+  , appComp
   , WorkerStates
   , initWorkerStates
   , module Data.Massiv.Core.Index
@@ -56,3 +57,10 @@
 import Data.Massiv.Core.List
 import Data.Massiv.Core.Exception
 
+
+-- | Append computation strategy using `Comp`'s `Monoid` instance.
+--
+-- @since 0.6.0
+appComp :: (Construct r ix e, Load r ix e) => Comp -> Array r ix e -> Array r ix e
+appComp comp arr = setComp (comp <> getComp arr) arr
+{-# INLINEABLE appComp #-}
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
@@ -8,7 +8,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Core.Common
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -83,6 +83,8 @@
   , Proxy(..)
   , Id(..)
   -- * Stateful Monads
+  , runST
+  , ST
   , MonadUnliftIO
   , MonadIO(liftIO)
   , PrimMonad(PrimState)
@@ -95,6 +97,7 @@
 import Control.Monad.Catch (MonadThrow(..))
 import Control.Monad.IO.Unlift (MonadIO(liftIO), MonadUnliftIO)
 import Control.Monad.Primitive
+import Control.Monad.ST
 import Control.Scheduler (Comp(..), Scheduler, WorkerStates, numWorkers,
                           scheduleWork, scheduleWork_, withScheduler_, trivialScheduler_)
 import Control.Scheduler.Global
@@ -223,7 +226,9 @@
   makeArrayLinear comp sz f = makeArray comp sz (f . toLinearIndex sz)
   {-# INLINE makeArrayLinear #-}
 
-
+  replicate :: Comp -> Sz ix -> e -> Array r ix e
+  replicate comp sz !e = makeArray comp sz (const e)
+  {-# INLINE replicate #-}
 
 class Index ix => Resize r ix where
   -- | /O(1)/ - Change the size of an array. Total number of elements should be the same, but it is
@@ -313,10 +318,6 @@
   loadArrayWithSetM scheduler arr uWrite _ = loadArrayM scheduler arr uWrite
   {-# INLINE loadArrayWithSetM #-}
 
-  defaultElement :: Array r ix e -> Maybe e
-  defaultElement _ = Nothing
-  {-# INLINE defaultElement #-}
-
   -- | /O(1)/ - Get the possible maximum size of an immutabe array. If the lookup of size
   -- in constant time is not possible, `Nothing` will be returned. This value will be used
   -- as the initial size of the mutable array into which the loading will happen.
@@ -370,9 +371,7 @@
     pure marr
   {-# INLINE unsafeLoadIntoM #-}
 
-{-# DEPRECATED defaultElement "This is no longer used by any of the loading functions and will be removed in massiv-0.6" #-}
 
-
 -- | Selects an optimal scheduler for the supplied strategy, but it works only in `IO`
 withMassivScheduler_ :: Comp -> (Scheduler IO () -> IO ()) -> IO ()
 withMassivScheduler_ comp f =
@@ -471,13 +470,18 @@
   --
   -- @since 0.3.0
   initializeNew :: PrimMonad m => Maybe e -> Sz ix -> m (MArray (PrimState m) r ix e)
-  initializeNew mdef sz = do
-    marr <- unsafeNew sz
-    case mdef of
-      Just val -> unsafeLinearSet marr 0 (SafeSz (totalElem sz)) val
-      Nothing  -> initialize marr
-    return marr
+  initializeNew Nothing sz = unsafeNew sz >>= \ma -> ma <$ initialize ma
+  initializeNew (Just e) sz = newMArray sz e
   {-# INLINE initializeNew #-}
+
+  -- | Create new mutable array while initializing all elements to the specified value.
+  --
+  -- @since 0.6.0
+  newMArray :: PrimMonad m => Sz ix -> e -> m (MArray (PrimState m) r ix e)
+  newMArray sz e = do
+    marr <- unsafeNew sz
+    marr <$ unsafeLinearSet marr 0 (SafeSz (totalElem sz)) e
+  {-# INLINE newMArray #-}
 
   -- | Set all cells in the mutable array within the range to a specified value.
   --
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
@@ -4,7 +4,7 @@
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Massiv.Core.Exception
--- Copyright   : (c) Alexey Kuleshevich 2019-2020
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Core/Index.hs b/src/Data/Massiv/Core/Index.hs
--- a/src/Data/Massiv/Core/Index.hs
+++ b/src/Data/Massiv/Core/Index.hs
@@ -5,7 +5,7 @@
 {-# LANGUAGE ExplicitNamespaces #-}
 -- |
 -- Module      : Data.Massiv.Core.Index
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
@@ -41,6 +41,7 @@
   , insertSzM
   , pullOutSzM
   , toLinearSz
+  , mkSzM
   -- ** Dimension
   , Dim(..)
   , Dimension(Dim1, Dim2, Dim3, Dim4, Dim5, DimN)
diff --git a/src/Data/Massiv/Core/Index/Internal.hs b/src/Data/Massiv/Core/Index/Internal.hs
--- a/src/Data/Massiv/Core/Index/Internal.hs
+++ b/src/Data/Massiv/Core/Index/Internal.hs
@@ -16,7 +16,7 @@
 {-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Internal
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
@@ -39,6 +39,7 @@
   , setSzM
   , insertSzM
   , pullOutSzM
+  , mkSzM
   , Dim(..)
   , Dimension(DimN)
   , pattern Dim1
@@ -62,6 +63,7 @@
 
 import Control.DeepSeq
 import Control.Exception (Exception(..), throw)
+import Control.Monad (when)
 import Control.Monad.Catch (MonadThrow(..))
 import Data.Coerce
 import Data.Massiv.Core.Iterator
@@ -129,7 +131,21 @@
   fromInteger = Sz . fromInteger
   {-# INLINE fromInteger #-}
 
+-- | Construct size from index while checking its correctness. Throws
+-- `SizeNegativeException` and `SizeOverflowException`.
+--
+-- @since 0.6.0
+mkSzM :: (Index ix, MonadThrow m) => ix -> m (Sz ix)
+mkSzM ix = do
+  let guardNegativeOverflow i !acc = do
+        when (i < 0) $ throwM $ SizeNegativeException (SafeSz ix)
+        let acc' = i * acc
+        when (acc' /= 0 && acc' < acc) $ throwM $ SizeOverflowException (SafeSz ix)
+        pure acc'
+  Sz ix <$ foldlIndex (\acc i -> acc >>= guardNegativeOverflow i) (pure 1) ix
 
+
+
 -- | Function for unwrapping `Sz`.
 --
 -- ==== __Example__
@@ -717,8 +733,8 @@
 instance NFData IndexException where
   rnf =
     \case
-      IndexZeroException i -> rnf i
-      IndexDimensionException i d -> i `deepseq` rnf d
+      IndexZeroException i           -> rnf i
+      IndexDimensionException i d    -> i `deepseq` rnf d
       IndexOutOfBoundsException sz i -> sz `deepseq` rnf i
 
 instance Exception IndexException
@@ -735,6 +751,14 @@
   SizeSubregionException :: Index ix => !(Sz ix) -> !ix -> !(Sz ix) -> SizeException
   -- | An array with the size cannot contain any elements.
   SizeEmptyException :: Index ix => !(Sz ix) -> SizeException
+  -- | Total number of elements is too large resulting in overflow.
+  --
+  -- @since 0.6.0
+  SizeOverflowException :: Index ix => !(Sz ix) -> SizeException
+  -- | At least one dimensions contain a negative value.
+  --
+  -- @since 0.6.0
+  SizeNegativeException :: Index ix => !(Sz ix) -> SizeException
 
 instance Eq SizeException where
   e1 == e2 =
@@ -746,15 +770,19 @@
       (SizeSubregionException sz1 i1 sz1', SizeSubregionException sz2 i2 sz2') ->
         show sz1 == show sz2 && show i1 == show i2 && show sz1' == show sz2'
       (SizeEmptyException sz1, SizeEmptyException sz2) -> show sz1 == show sz2
+      (SizeOverflowException sz1, SizeOverflowException sz2) -> show sz1 == show sz2
+      (SizeNegativeException sz1, SizeNegativeException sz2) -> show sz1 == show sz2
       _ -> False
 
 instance NFData SizeException where
   rnf =
     \case
-      SizeMismatchException sz sz' -> sz `deepseq` rnf sz'
+      SizeMismatchException sz sz'         -> sz `deepseq` rnf sz'
       SizeElementsMismatchException sz sz' -> sz `deepseq` rnf sz'
-      SizeSubregionException sz i sz' -> sz `deepseq` i `deepseq` rnf sz'
-      SizeEmptyException sz -> rnf sz
+      SizeSubregionException sz i sz'      -> sz `deepseq` i `deepseq` rnf sz'
+      SizeEmptyException sz                -> rnf sz
+      SizeOverflowException sz             -> rnf sz
+      SizeNegativeException sz             -> rnf sz
 
 instance Exception SizeException
 
@@ -769,6 +797,10 @@
     show sz' ++ ") is to small for " ++ show ix ++ " (" ++ show sz ++ ")"
   show (SizeEmptyException sz) =
     "SizeEmptyException: (" ++ show sz ++ ") corresponds to an empty array"
+  show (SizeOverflowException sz) =
+    "SizeOverflowException: (" ++ show sz ++ ") is too big"
+  show (SizeNegativeException sz) =
+    "SizeNegativeException: (" ++ show sz ++ ") contains negative value"
   showsPrec n exc = showsPrecWrapped n (show exc ++)
 
 -- | Exception that can happen upon conversion of a ragged type array into the rectangular kind. Which
diff --git a/src/Data/Massiv/Core/Index/Ix.hs b/src/Data/Massiv/Core/Index/Ix.hs
--- a/src/Data/Massiv/Core/Index/Ix.hs
+++ b/src/Data/Massiv/Core/Index/Ix.hs
@@ -12,7 +12,7 @@
 {-# LANGUAGE UndecidableInstances #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Ix
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -373,6 +373,8 @@
   {-# INLINE [1] repairIndex #-}
 
 -- | Constraint synonym that encapsulates all constraints needed for dimension 4 and higher.
+--
+-- @since 0.6.0
 type HighIxN n
    = (4 <= n, KnownNat n, KnownNat (n - 1), Index (Ix (n - 1)), IxN (n - 1) ~ Ix (n - 1))
 
diff --git a/src/Data/Massiv/Core/Index/Stride.hs b/src/Data/Massiv/Core/Index/Stride.hs
--- a/src/Data/Massiv/Core/Index/Stride.hs
+++ b/src/Data/Massiv/Core/Index/Stride.hs
@@ -3,7 +3,7 @@
 {-# LANGUAGE PatternSynonyms #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Stride
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Core/Index/Tuple.hs b/src/Data/Massiv/Core/Index/Tuple.hs
--- a/src/Data/Massiv/Core/Index/Tuple.hs
+++ b/src/Data/Massiv/Core/Index/Tuple.hs
@@ -6,7 +6,7 @@
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Core.Index.Tuple
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <alexey@kuleshevi.ch>
 -- Stability   : experimental
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
@@ -1,7 +1,7 @@
 {-# LANGUAGE BangPatterns #-}
 -- |
 -- Module      : Data.Massiv.Core.Iterator
--- Copyright   : (c) Alexey Kuleshevich 2018-2020
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Core/List.hs b/src/Data/Massiv/Core/List.hs
--- a/src/Data/Massiv/Core/List.hs
+++ b/src/Data/Massiv/Core/List.hs
@@ -10,7 +10,7 @@
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Massiv.Core.List
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2018-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
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,25 +1,30 @@
+{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 -- |
 -- Module      : Data.Massiv.Core.Operations
--- Copyright   : (c) Alexey Kuleshevich 2018-2019
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
 -- Portability : non-portable
 module Data.Massiv.Core.Operations
-  ( Numeric(..)
+  ( FoldNumeric(..)
+  , defaultPowerSumArray
+  , defaultUnsafeDotProduct
+  , defaultFoldArray
+  , Numeric(..)
+  , defaultUnsafeLiftArray
+  , defaultUnsafeLiftArray2
   , NumericFloat(..)
   ) where
 
 import Data.Massiv.Core.Common
 
 
-class Num e => Numeric r e where
-
-  {-# MINIMAL foldArray, unsafeLiftArray, unsafeLiftArray2 #-}
+class Num e => FoldNumeric r e where
 
   -- | Compute sum of all elements in the array
   --
@@ -39,17 +44,50 @@
   --
   -- @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
   --
   -- @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 #-}
 
+  -- | Fold over an array
+  --
+  -- @since 0.5.6
+  foldArray :: Index ix => (e -> e -> e) -> e -> Array r ix e -> e
 
+
+defaultUnsafeDotProduct ::
+     (Num e, Source r ix e) => Array r ix e -> Array r ix e -> e
+defaultUnsafeDotProduct 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 defaultUnsafeDotProduct #-}
+
+defaultPowerSumArray :: (Source r ix e, Num e) => Array r ix e -> Int -> e
+defaultPowerSumArray 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 defaultPowerSumArray #-}
+
+defaultFoldArray :: Source r ix e => (e -> e -> e) -> e -> Array r ix e -> e
+defaultFoldArray 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 defaultFoldArray #-}
+
+class FoldNumeric r e => Numeric r e where
+
+  {-# MINIMAL unsafeLiftArray, unsafeLiftArray2 #-}
+
   plusScalar :: Index ix => Array r ix e -> e -> Array r ix e
   plusScalar arr e = unsafeLiftArray (+ e) arr
   {-# INLINE plusScalar #-}
@@ -90,15 +128,29 @@
   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 => (e -> e) -> Array r ix e -> Array r ix e
 
   unsafeLiftArray2 :: Index ix => (e -> e -> e) -> Array r ix e -> Array r ix e -> Array r ix e
 
+
+defaultUnsafeLiftArray ::
+     (Construct r ix e, Source r ix e) => (e -> e) -> Array r ix e -> Array r ix e
+defaultUnsafeLiftArray f arr = makeArrayLinear (getComp arr) (size arr) (f . unsafeLinearIndex arr)
+{-# INLINE defaultUnsafeLiftArray #-}
+
+
+defaultUnsafeLiftArray2 ::
+     (Construct r ix e, Source r ix e)
+  => (e -> e -> e)
+  -> Array r ix e
+  -> Array r ix e
+  -> Array r ix e
+defaultUnsafeLiftArray2 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 defaultUnsafeLiftArray2 #-}
 
 
 class (Numeric r e, Floating e) => NumericFloat r e where
diff --git a/src/Data/Massiv/Vector.hs b/src/Data/Massiv/Vector.hs
--- a/src/Data/Massiv/Vector.hs
+++ b/src/Data/Massiv/Vector.hs
@@ -3,7 +3,7 @@
 {-# OPTIONS_GHC -fno-warn-duplicate-exports #-}
 -- |
 -- Module      : Data.Massiv.Vector
--- Copyright   : (c) Alexey Kuleshevich 2020
+-- Copyright   : (c) Alexey Kuleshevich 2020-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
diff --git a/src/Data/Massiv/Vector/Stream.hs b/src/Data/Massiv/Vector/Stream.hs
--- a/src/Data/Massiv/Vector/Stream.hs
+++ b/src/Data/Massiv/Vector/Stream.hs
@@ -4,12 +4,11 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TupleSections #-}
 {-# OPTIONS_HADDOCK hide, not-home #-}
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Massiv.Vector.Stream
--- Copyright   : (c) Alexey Kuleshevich 2019-2020
+-- Copyright   : (c) Alexey Kuleshevich 2019-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental
@@ -132,7 +131,7 @@
 import qualified Control.Monad as M
 import Control.Monad.ST
 import qualified Data.Foldable as F
-import Data.Massiv.Core.Common hiding (empty, singleton)
+import Data.Massiv.Core.Common hiding (empty, singleton, replicate)
 import Data.Maybe (catMaybes)
 import qualified Data.Traversable as Traversable (traverse)
 import qualified Data.Vector.Fusion.Bundle.Monadic as B
diff --git a/src/Data/Massiv/Vector/Unsafe.hs b/src/Data/Massiv/Vector/Unsafe.hs
--- a/src/Data/Massiv/Vector/Unsafe.hs
+++ b/src/Data/Massiv/Vector/Unsafe.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 -- |
 -- Module      : Data.Massiv.Vector.Unsafe
--- Copyright   : (c) Alexey Kuleshevich 2020
+-- Copyright   : (c) Alexey Kuleshevich 2020-2021
 -- License     : BSD3
 -- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
 -- Stability   : experimental