diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,11 @@
+# 0.5.2
+
+* Addition of `lowerTriangular` and `upperTriangular`
+* Relax `identityMatrix` type to return an array of any `Num` type, not just `Int`.
+* Addition of `unsafeMakeLoadArrayAdjusted`
+* Add matrix-vector product (`(#>)`)
+* Addition of `siterate`
+
 # 0.5.1
 
 * Fix `sfromListN` accepting a plain `Int` instead of `Sz1`, as well as switch to upper bound.
diff --git a/massiv.cabal b/massiv.cabal
--- a/massiv.cabal
+++ b/massiv.cabal
@@ -1,5 +1,5 @@
 name:                massiv
-version:             0.5.1.0
+version:             0.5.2.0
 synopsis:            Massiv (Массив) is an Array Library.
 description:         Multi-dimensional Arrays with fusion, stencils and parallel computation.
 homepage:            https://github.com/lehins/massiv
diff --git a/src/Data/Massiv/Array/Delayed/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
@@ -24,6 +24,7 @@
   , makeLoadArrayS
   , makeLoadArray
   , unsafeMakeLoadArray
+  , unsafeMakeLoadArrayAdjusted
   , fromStrideLoad
   , appendOuterM
   , concatOuterM
@@ -231,8 +232,7 @@
   -> Maybe e
   -- ^ An element to use for initialization of the mutable array that will be created in
   -- the future
-  -> (forall m. Monad m =>
-                  Scheduler m () -> Int -> (Int -> e -> m ()) -> m ())
+  -> (forall m. Monad m => Scheduler m () -> Int -> (Int -> e -> m ()) -> m ())
   -- ^ This function accepts:
   --
   -- * A scheduler that can be used for parallelization of loading
@@ -244,6 +244,21 @@
   -> Array DL ix e
 unsafeMakeLoadArray = DLArray
 {-# INLINE unsafeMakeLoadArray #-}
+
+-- | Same as `unsafeMakeLoadArray`, except will ensure that starting index is correctly
+-- adjusted. Which means the writing function gets one less argument.
+--
+-- @since 0.5.2
+unsafeMakeLoadArrayAdjusted ::
+     Comp
+  -> Sz ix
+  -> Maybe e
+  -> (forall m. Monad m => Scheduler m () -> (Int -> e -> m ()) -> m ())
+  -> Array DL ix e
+unsafeMakeLoadArrayAdjusted comp sz mDefVal writer =
+  DLArray comp sz mDefVal $ \scheduler !startAt uWrite ->
+    writer scheduler (\i -> uWrite (startAt + i))
+{-# INLINE unsafeMakeLoadArrayAdjusted #-}
 
 -- | Convert any `Load`able array into `DL` representation.
 --
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
@@ -23,9 +23,12 @@
   , (.*)
   , (*.)
   , (.^)
+  , (#>)
   , (|*|)
   , multiplyTransposed
   , identityMatrix
+  , lowerTriangular
+  , upperTriangular
   , negateA
   , absA
   , signumA
@@ -77,6 +80,7 @@
 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
 import Data.Massiv.Core.Operations
@@ -195,7 +199,9 @@
 (.^) = powerPointwise
 {-# INLINE (.^) #-}
 
--- | Perform matrix multiplication. Inner dimensions must agree, otherwise `SizeMismatchException`.
+-- | Matrix multiplication
+--
+-- Inner dimensions must agree, otherwise `SizeMismatchException`.
 (|*|) ::
      (Mutable r Ix2 e, Source r' Ix2 e, OuterSlice r Ix2 e, Source (R r) Ix1 e, Num e, MonadThrow m)
   => Array r Ix2 e
@@ -209,6 +215,26 @@
     multiplyTransposedFused arr1 (convert arr2)
  #-}
 
+-- | 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 (#>) #-}
+
+
+
 multiplyTransposedFused ::
      ( Mutable r Ix2 e
      , OuterSlice r Ix2 e
@@ -238,8 +264,7 @@
     arr2' = compute $ transpose arr2
 {-# INLINE multArrs #-}
 
--- | It is quite often that second matrix gets transposed before multiplication (eg. A * A'), but
--- due to layout of data in memory it is more efficient to transpose the second array again.
+-- | Computes the matrix-matrix transposed product (i.e. A * A')
 multiplyTransposed ::
      ( Manifest r Ix2 e
      , OuterSlice r Ix2 e
@@ -276,9 +301,62 @@
 --   ]
 --
 -- @since 0.3.6
-identityMatrix :: Sz1 -> Array DL Ix2 Int
-identityMatrix (Sz n) = makeLoadArrayS (Sz2 n n) 0 $ \ w -> loopM_ 0 (< n) (+1) $ \ i -> w (i :. i) 1
+identityMatrix :: Num e => Sz1 -> Matrix DL e
+identityMatrix (Sz n) =
+  makeLoadArrayS (Sz2 n n) 0 $ \ w -> loopM_ 0 (< n) (+1) $ \ i -> w (i :. i) 1
 {-# INLINE identityMatrix #-}
+
+-- | Create a lower triangular (L in LU decomposition) matrix of size @NxN@
+--
+-- ==== __Example__
+--
+-- >>> import Data.Massiv.Array
+-- >>> lowerTriangular Seq 5 (\(i :. j) -> i + j)
+-- Array DL Seq (Sz (5 :. 5))
+--   [ [ 0, 0, 0, 0, 0 ]
+--   , [ 1, 2, 0, 0, 0 ]
+--   , [ 2, 3, 4, 0, 0 ]
+--   , [ 3, 4, 5, 6, 0 ]
+--   , [ 4, 5, 6, 7, 8 ]
+--   ]
+--
+-- @since 0.5.2
+lowerTriangular :: Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
+lowerTriangular comp (Sz1 n) f =
+  let sz = Sz2 n n
+   in unsafeMakeLoadArrayAdjusted comp sz (Just 0) $ \scheduler wr ->
+        forM_ (0 ..: n) $ \i ->
+          scheduleWork scheduler $
+          forM_ (0 ... i) $ \j ->
+            let ix = i :. j
+             in wr (toLinearIndex sz ix) (f ix)
+{-# INLINE lowerTriangular #-}
+
+-- | Create an upper triangular (U in LU decomposition) matrix of size @NxN@
+--
+-- ==== __Example__
+--
+-- >>> import Data.Massiv.Array
+-- >>> upperTriangular Par 5 (\(i :. j) -> i + j)
+-- Array DL Par (Sz (5 :. 5))
+--   [ [ 0, 1, 2, 3, 4 ]
+--   , [ 0, 2, 3, 4, 5 ]
+--   , [ 0, 0, 4, 5, 6 ]
+--   , [ 0, 0, 0, 6, 7 ]
+--   , [ 0, 0, 0, 0, 8 ]
+--   ]
+--
+-- @since 0.5.2
+upperTriangular :: Num e => Comp -> Sz1 -> (Ix2 -> e) -> Matrix DL e
+upperTriangular comp (Sz1 n) f =
+  let sz = Sz2 n n
+   in unsafeMakeLoadArrayAdjusted comp sz (Just 0) $ \scheduler wr ->
+        forM_ (0 ..: n) $ \i ->
+          scheduleWork scheduler $
+          forM_ (i ..: n) $ \j ->
+            let ix = i :. j
+             in wr (toLinearIndex sz ix) (f ix)
+{-# INLINE upperTriangular #-}
 
 
 negateA :: (Index ix, Numeric r e) => Array r ix e -> Array r ix e
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
@@ -14,6 +14,7 @@
 module Data.Massiv.Array.Unsafe
   ( -- * Creation
     unsafeMakeLoadArray
+  , unsafeMakeLoadArrayAdjusted
     -- * Indexing
   , Sz(SafeSz)
   , Stride(SafeStride)
@@ -86,7 +87,7 @@
   ) where
 
 import Data.Massiv.Array.Delayed.Pull (D)
-import Data.Massiv.Array.Delayed.Push (unsafeMakeLoadArray)
+import Data.Massiv.Array.Delayed.Push (unsafeMakeLoadArray, unsafeMakeLoadArrayAdjusted)
 import Data.Massiv.Array.Manifest.Boxed
 import Data.Massiv.Array.Manifest.Primitive
 import Data.Massiv.Array.Manifest.Storable
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
@@ -73,6 +73,7 @@
   , sgenerate
   -- , iterateN
   -- , iiterateN
+  , siterate
   , siterateN
   -- ** Monadic initialization
   , sreplicateM
@@ -1016,10 +1017,28 @@
 {-# INLINE sgenerate #-}
 
 
--- | Create a delayed stream vector of length @n@ by repeatedly apply a function to the
+-- | Create a delayed stream vector of infinite length by repeatedly applying a function to the
 -- initial value.
 --
 -- ==== __Examples__
+--
+-- >>> stake 10 $ siterate succ 'a'
+-- Array DS Seq (Sz1 10)
+--   [ 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j' ]
+--
+-- @since 0.5.2
+siterate :: (e -> e) -> e -> Vector DS e
+siterate f = fromSteps . S.unfoldr (\a -> Just (a, f a))
+{-# INLINE siterate #-}
+
+-- | Create a delayed stream vector of length @n@ by repeatedly applying a function to the
+-- initial value.
+--
+-- ==== __Examples__
+--
+-- >>> siterateN 10 succ 'a'
+-- Array DS Seq (Sz1 10)
+--   [ 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j' ]
 --
 -- @since 0.5.0
 siterateN :: Sz1 -> (e -> e) -> e -> Vector DS e