From ac3ddc20543208ad0550a7a50e3f306d1bf5825d Mon Sep 17 00:00:00 2001
From: Nick Robinson <nicholas.robinson@invenialabs.co.uk>
Date: Tue, 2 Feb 2021 14:54:57 +0000
Subject: [PATCH 1/5] Allows blocks to be <:AbstractVector or <:Tuple

- This is a breaking change, as the type param `V` is now different,
  being the type of the `blocks` collection, rather than the `eltype`
  of that collection.
- This allows allows the collection to be something other than a
  `Vector` e.g.
  - a `CuArray` of matrices (allowing use on GPU hopefully)
  - a `Tuple` with different concrete subtypes of `AbstractMatrix{T}`
    (avoiding heterogeneous block to being abstractly typed
    i.e. a way to avoid `Vector{AbstractMatrix{T}}`)
---
 Project.toml          |  2 +-
 src/blockdiagonal.jl  | 23 ++++++++++++++---------
 src/chainrules.jl     |  4 ++--
 src/linalg.jl         | 13 ++++++++++++-
 test/base_maths.jl    | 16 +++++++++++-----
 test/blockdiagonal.jl | 16 ++++++++++++----
 test/linalg.jl        | 26 ++++++++++++++++----------
 7 files changed, 68 insertions(+), 32 deletions(-)

diff --git a/Project.toml b/Project.toml
index a3b029d..ee4dbff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockDiagonals"
 uuid = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 authors = ["Invenia Technical Computing Corporation"]
-version = "0.1.14"
+version = "0.2.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
diff --git a/src/blockdiagonal.jl b/src/blockdiagonal.jl
index aefd77d..f974700 100644
--- a/src/blockdiagonal.jl
+++ b/src/blockdiagonal.jl
@@ -1,19 +1,22 @@
 # Core functionality for the `BlockDiagonal` type
 
 """
-    BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
+    BlockDiagonal{T, V} <: AbstractMatrix{T}
+    BlockDiagonal(blocks::V) -> BlockDiagonal{T,V}
 
 A matrix with matrices on the diagonal, and zeros off the diagonal.
-"""
-struct BlockDiagonal{T, V<:AbstractMatrix{T}} <: AbstractMatrix{T}
-    blocks::Vector{V}
 
-    function BlockDiagonal{T, V}(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
-        return new{T, V}(blocks)
-    end
+!!! info "`V` type"
+    `blocks::V` should be a `Tuple` or `AbstractVector` where each component (each block) is
+    `<:AbstractMatrix{T}` for some common element type `T`.
+"""
+struct BlockDiagonal{T, V} <: AbstractMatrix{T}
+    blocks::V
 end
 
-function BlockDiagonal(blocks::Vector{V}) where {T, V<:AbstractMatrix{T}}
+function BlockDiagonal(blocks::V) where {
+    T, V<:Union{Tuple{Vararg{<:AbstractMatrix{T}}}, AbstractVector{<:AbstractMatrix{T}}}
+}
     return BlockDiagonal{T, V}(blocks)
 end
 
@@ -151,7 +154,9 @@ function _block_indices(B::BlockDiagonal, i::Integer, j::Integer)
         p += 1
         j -= ncols[p]
     end
-    i -= sum(nrows[1:(p-1)])
+    if !isempty(nrows[1:(p-1)])
+        i -= sum(nrows[1:(p-1)])
+    end
     # if row `i` outside of block `p`, set `p` to place-holder value `-1`
     if i <= 0 || i > nrows[p]
         p = -1
diff --git a/src/chainrules.jl b/src/chainrules.jl
index c640488..db0e4ef 100644
--- a/src/chainrules.jl
+++ b/src/chainrules.jl
@@ -1,5 +1,5 @@
 # constructor
-function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::Vector{V}) where {V}
+function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::V) where {V<:AbstractVector}
     BlockDiagonal_pullback(Δ::Composite) = (NO_FIELDS, Δ.blocks)
     return BlockDiagonal(blocks), BlockDiagonal_pullback
 end
@@ -27,7 +27,7 @@ function ChainRulesCore.rrule(
         ::typeof(*),
         bm::BlockDiagonal{T, V},
         v::StridedVector{T}
-    ) where {T<:Union{Real, Complex}, V<:Matrix{T}}
+    ) where {T<:Union{Real, Complex}, V<:Vector{Matrix{T}}}
 
     y = bm * v
 
diff --git a/src/linalg.jl b/src/linalg.jl
index 2ea9b42..183f38a 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -76,7 +76,7 @@ svdvals_blockwise(B::BlockDiagonal) = mapreduce(svdvals, vcat, blocks(B))
 LinearAlgebra.svdvals(B::BlockDiagonal) = sort!(svdvals_blockwise(B); rev=true)
 
 # `B = U * Diagonal(S) * Vt` with `U` and `Vt` `BlockDiagonal` (`S` only sorted block-wise).
-function svd_blockwise(B::BlockDiagonal{T}; full::Bool=false) where T
+function svd_blockwise(B::BlockDiagonal{T, <:AbstractVector}; full::Bool=false) where T
     U = Matrix{float(T)}[]
     S = Vector{float(T)}()
     Vt = Matrix{float(T)}[]
@@ -88,6 +88,17 @@ function svd_blockwise(B::BlockDiagonal{T}; full::Bool=false) where T
     end
     return BlockDiagonal(U), S, BlockDiagonal(Vt)
 end
+function svd_blockwise(B::BlockDiagonal{T, <:Tuple}; full::Bool=false) where T
+    S = Vector{float(T)}()
+    U_Vt = ntuple(length(blocks(B))) do i
+        F = svd(getblock(B, i), full=full)
+        append!(S, F.S)
+        (F.U, F.Vt)
+    end
+    U = first.(U_Vt)
+    Vt = last.(U_Vt)
+    return BlockDiagonal(U), S, BlockDiagonal(Vt)
+end
 
 function LinearAlgebra.svd(B::BlockDiagonal; full::Bool=false)::SVD
     U, S, Vt = svd_blockwise(B, full=full)
diff --git a/test/base_maths.jl b/test/base_maths.jl
index fb5be5f..caaedf6 100644
--- a/test/base_maths.jl
+++ b/test/base_maths.jl
@@ -7,9 +7,14 @@ using Test
     rng = MersenneTwister(123456)
     N1, N2, N3 = 3, 4, 5
     N = N1 + N2 + N3
-    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
-    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
-    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
+    blocks1 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
+    blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
+    blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
+
+    @testset "$T" for (T, (b1, b2, b3)) in (
+        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
+        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
+    )
     A = rand(rng, N, N + N1)
     B = rand(rng, N + N1, N + N2)
     A′, B′ = A', B'
@@ -127,8 +132,8 @@ using Test
         end
 
         @testset "Non-Square BlockDiagonal * Non-Square BlockDiagonal" begin
-    	    b4 = BlockDiagonal([ones(2, 4), 2 * ones(3, 2)])
-            b5 = BlockDiagonal([3 * ones(2, 2), 2 * ones(4, 1)])
+            b4 = BlockDiagonal(T([ones(2, 4), 2 * ones(3, 2)]))
+            b5 = BlockDiagonal(T([3 * ones(2, 2), 2 * ones(4, 1)]))
 
             @test b4 * b5 isa Array
             @test b4 * b5 == [6 * ones(2, 2) 4 * ones(2, 1); zeros(3, 2) 8 * ones(3, 1)]
@@ -138,3 +143,4 @@ using Test
         end
     end  # Multiplication
 end
+end
diff --git a/test/blockdiagonal.jl b/test/blockdiagonal.jl
index 464e69b..16963ac 100644
--- a/test/blockdiagonal.jl
+++ b/test/blockdiagonal.jl
@@ -7,9 +7,14 @@ using Test
     rng = MersenneTwister(123456)
     N1, N2, N3 = 3, 4, 5
     N = N1 + N2 + N3
-    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
-    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
-    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
+    blocks1 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
+    blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
+    blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
+
+    @testset "$T" for (T, (b1, b2, b3)) in (
+        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
+        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
+    )
     A = rand(rng, N, N + N1)
     B = rand(rng, N + N1, N + N2)
     A′, B′ = A', B'
@@ -42,8 +47,10 @@ using Test
         end
 
         @testset "parent" begin
-            @test parent(b1) isa Vector{<:AbstractMatrix}
+            @test parent(b1) isa Union{Tuple,AbstractVector}
+            @test eltype(parent(b1)) <: AbstractMatrix
             @test parent(BlockDiagonal([X, Y])) == [X, Y]
+            @test parent(BlockDiagonal((X, Y))) == (X, Y)
         end
 
         @testset "similar" begin
@@ -117,3 +124,4 @@ using Test
         @test_throws DimensionMismatch copy!(b2, b1)
     end
 end
+end
diff --git a/test/linalg.jl b/test/linalg.jl
index 86634a1..03c1583 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -8,9 +8,14 @@ using Test
     rng = MersenneTwister(123456)
     N1, N2, N3 = 3, 4, 5
     N = N1 + N2 + N3
-    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
-    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
-    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
+    blocks1 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
+    blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
+    blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
+
+    @testset "$T" for (T, (b1, b2, b3)) in (
+        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
+        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
+    )
     A = rand(rng, N, N + N1)
     B = rand(rng, N + N1, N + N2)
     A′, B′ = A', B'
@@ -143,7 +148,7 @@ using Test
 
         @testset "eigvals on LinearAlgebra types" begin
             # `eigvals` has different methods for different types, e.g. Hermitian
-            b_herm = BlockDiagonal([Hermitian(rand(rng, 3, 3) + I) for _ in 1:3])
+            b_herm = BlockDiagonal(T(Hermitian(rand(rng, 3, 3) + I) for _ in 1:3))
             @test eigvals(b_herm) ≈ eigvals(Matrix(b_herm))
             @test eigvals(b_herm, 1.0, 2.0) ≈ eigvals(Hermitian(Matrix(b_herm)), 1.0, 2.0)
         end
@@ -163,12 +168,12 @@ using Test
               0.0 1.0  5.0
               0.0 0.0  3.0]
 
-        B = BlockDiagonal([X, X])
+        B = BlockDiagonal(T([X, X]))
         C = cholesky(B)
         @test C isa Cholesky{Float64, <:BlockDiagonal{Float64}}
         @test C.U ≈ cholesky(Matrix(B)).U
-        @test C.U ≈ BlockDiagonal([U, U])
-        @test C.L ≈ BlockDiagonal([U', U'])
+        @test C.U ≈ BlockDiagonal(T([U, U]))
+        @test C.L ≈ BlockDiagonal(T([U', U']))
         @test C.UL ≈ C.U
         @test C.uplo === 'U'
         @test C.info == 0
@@ -176,8 +181,8 @@ using Test
         M = BlockDiagonal(map(Matrix, blocks(C.L)))
         C = Cholesky(M, 'L', 0)
         @test C.U ≈ cholesky(Matrix(B)).U
-        @test C.U ≈ BlockDiagonal([U, U])
-        @test C.L ≈ BlockDiagonal([U', U'])
+        @test C.U ≈ BlockDiagonal(T([U, U]))
+        @test C.L ≈ BlockDiagonal(T([U', U']))
         @test C.UL ≈ C.L
         @test C.uplo === 'L'
         @test C.info == 0
@@ -187,7 +192,7 @@ using Test
         X = [  4  12 -16
               12  37 -43
              -16 -43  98]
-        B = BlockDiagonal([X, X])
+        B = BlockDiagonal(T([X, X]))
 
         @testset "full=$full" for full in (true, false)
 
@@ -230,3 +235,4 @@ using Test
         end
     end  # SVD
 end
+end

From 32f9d07995d8bfee83d386e3a525d671ad0d13a7 Mon Sep 17 00:00:00 2001
From: Nick Robinson <nicholas.robinson@invenialabs.co.uk>
Date: Mon, 15 Feb 2021 19:30:35 +0000
Subject: [PATCH 2/5] Simplify test setup for base_maths.jl

---
 test/base_maths.jl | 51 ++++++++++++++++++++++------------------------
 1 file changed, 24 insertions(+), 27 deletions(-)

diff --git a/test/base_maths.jl b/test/base_maths.jl
index caaedf6..2af948f 100644
--- a/test/base_maths.jl
+++ b/test/base_maths.jl
@@ -5,38 +5,31 @@ using Test
 
 @testset "base_maths.jl" begin
     rng = MersenneTwister(123456)
-    N1, N2, N3 = 3, 4, 5
-    N = N1 + N2 + N3
-    blocks1 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
-    blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
-    blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
-
-    @testset "$T" for (T, (b1, b2, b3)) in (
-        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
-        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
-    )
-    A = rand(rng, N, N + N1)
-    B = rand(rng, N + N1, N + N2)
-    A′, B′ = A', B'
-    a = rand(rng, N)
-    b = rand(rng, N + N1)
+    blocks1 = [rand(rng, 3, 3), rand(rng, 4, 4)]
+    blocks2 = [rand(rng, 3, 3), rand(rng, 5, 5)]
+
+    @testset for V in (Tuple, Vector)
+    b1 = BlockDiagonal(V(blocks1))
+    b2 = BlockDiagonal(V(blocks2))
+    N = size(b1, 1)
+    A = rand(rng, N, N + 1)
 
     @testset "Addition" begin
         @testset "BlockDiagonal + BlockDiagonal" begin
             @test b1 + b1 isa BlockDiagonal
             @test Matrix(b1 + b1) == Matrix(b1) + Matrix(b1)
-            @test_throws DimensionMismatch b1 + b3
+            @test_throws DimensionMismatch b1 + b2
         end
 
         @testset "BlockDiagonal + Matrix" begin
             @test b1 + Matrix(b1) isa Matrix
             @test b1 + Matrix(b1) == b1 + b1
-            @test_throws DimensionMismatch b1 + Matrix(b3)
+            @test_throws DimensionMismatch b1 + Matrix(b2)
 
             # Matrix + BlockDiagonal
             @test Matrix(b1) + b1 isa Matrix
             @test Matrix(b1) + b1 == b1 + b1
-            @test_throws DimensionMismatch Matrix(b1) + b3
+            @test_throws DimensionMismatch Matrix(b1) + b2
 
             # If the AbstractMatrix is diagonal, we should return a BlockDiagonal.
             # Test the StridedMatrix method.
@@ -50,7 +43,7 @@ using Test
 
         @testset "BlockDiagonal + Diagonal" begin
             D = Diagonal(randn(rng, N))
-            D′ = Diagonal(randn(rng, N + N1))
+            D′ = Diagonal(randn(rng, N + 1))
 
             @test b1 + D isa BlockDiagonal
             @test b1 + D == Matrix(b1) + D
@@ -73,11 +66,10 @@ using Test
     end  # Addition
 
     @testset "Multiplication" begin
-
         @testset "BlockDiagonal * BlockDiagonal" begin
             @test b1 * b1 isa BlockDiagonal
             @test Matrix(b1 * b1) ≈ Matrix(b1) * Matrix(b1)
-            @test_throws DimensionMismatch b3 * b1
+            @test_throws DimensionMismatch b2 * b1
         end
 
         @testset "BlockDiagonal * Number" begin
@@ -88,11 +80,14 @@ using Test
         end
 
         @testset "BlockDiagonal * Vector" begin
+            a = rand(rng, N)
             @test b1 * a isa Vector
             @test b1 * a ≈ Matrix(b1) * a
+            b = rand(rng, N + 1)
             @test_throws DimensionMismatch b1 * b
         end
         @testset "Vector^T * BlockDiagonal" begin
+            a = rand(rng, N)
             @test a' * b1 isa Adjoint{<:Number, <:Vector}
             @test transpose(a) * b1 isa Transpose{<:Number, <:Vector}
             @test a' * b1 ≈ a' * Matrix(b1)
@@ -102,11 +97,13 @@ using Test
         @testset "BlockDiagonal * Matrix" begin
             @test b1 * A isa Matrix
             @test b1 * A ≈ Matrix(b1) * A
+
+            B = rand(rng, N + 1, N)
             @test_throws DimensionMismatch b1 * B
 
             # Matrix * BlockDiagonal
-            @test A′ * b1 isa Matrix
-            @test A′ * b1 ≈ A′ * Matrix(b1)
+            @test A' * b1 isa Matrix
+            @test A' * b1 ≈ A' * Matrix(b1)
             @test_throws DimensionMismatch A * b1
 
             # degenerate cases
@@ -119,7 +116,7 @@ using Test
 
         @testset "BlockDiagonal * Diagonal" begin
             D = Diagonal(randn(rng, N))
-            D′ = Diagonal(randn(rng, N + N1))
+            D′ = Diagonal(randn(rng, N + 1))
 
             @test b1 * D isa BlockDiagonal
             @test b1 * D ≈ Matrix(b1) * D
@@ -132,8 +129,8 @@ using Test
         end
 
         @testset "Non-Square BlockDiagonal * Non-Square BlockDiagonal" begin
-            b4 = BlockDiagonal(T([ones(2, 4), 2 * ones(3, 2)]))
-            b5 = BlockDiagonal(T([3 * ones(2, 2), 2 * ones(4, 1)]))
+            b4 = BlockDiagonal(V([ones(2, 4), 2 * ones(3, 2)]))
+            b5 = BlockDiagonal(V([3 * ones(2, 2), 2 * ones(4, 1)]))
 
             @test b4 * b5 isa Array
             @test b4 * b5 == [6 * ones(2, 2) 4 * ones(2, 1); zeros(3, 2) 8 * ones(3, 1)]
@@ -142,5 +139,5 @@ using Test
             @test sum(size.(b5.blocks, 2)) == size(b4 * b5, 2)
         end
     end  # Multiplication
-end
+end  # V
 end

From 8b3235c606be8703cbb598567e156e9a40be71c7 Mon Sep 17 00:00:00 2001
From: Nick Robinson <nicholas.robinson@invenialabs.co.uk>
Date: Mon, 15 Feb 2021 19:40:29 +0000
Subject: [PATCH 3/5] Simplify test setup for linalg.jl

---
 test/linalg.jl | 30 +++++++++++-------------------
 1 file changed, 11 insertions(+), 19 deletions(-)

diff --git a/test/linalg.jl b/test/linalg.jl
index 03c1583..f5a31e4 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -7,20 +7,12 @@ using Test
 @testset "linalg.jl" begin
     rng = MersenneTwister(123456)
     N1, N2, N3 = 3, 4, 5
-    N = N1 + N2 + N3
     blocks1 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
     blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
-    blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
-
-    @testset "$T" for (T, (b1, b2, b3)) in (
-        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
-        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
-    )
-    A = rand(rng, N, N + N1)
-    B = rand(rng, N + N1, N + N2)
-    A′, B′ = A', B'
-    a = rand(rng, N)
-    b = rand(rng, N + N1)
+
+    @testset for V in (Tuple, Vector)
+    b1 = BlockDiagonal(V(blocks1))
+    b2 = BlockDiagonal(V(blocks2))
 
     @testset "mul!" begin
         c = similar(b1)
@@ -148,7 +140,7 @@ using Test
 
         @testset "eigvals on LinearAlgebra types" begin
             # `eigvals` has different methods for different types, e.g. Hermitian
-            b_herm = BlockDiagonal(T(Hermitian(rand(rng, 3, 3) + I) for _ in 1:3))
+            b_herm = BlockDiagonal(V([Hermitian(rand(rng, 3, 3) + I) for _ in 1:3]))
             @test eigvals(b_herm) ≈ eigvals(Matrix(b_herm))
             @test eigvals(b_herm, 1.0, 2.0) ≈ eigvals(Hermitian(Matrix(b_herm)), 1.0, 2.0)
         end
@@ -168,12 +160,12 @@ using Test
               0.0 1.0  5.0
               0.0 0.0  3.0]
 
-        B = BlockDiagonal(T([X, X]))
+        B = BlockDiagonal(V([X, X]))
         C = cholesky(B)
         @test C isa Cholesky{Float64, <:BlockDiagonal{Float64}}
         @test C.U ≈ cholesky(Matrix(B)).U
-        @test C.U ≈ BlockDiagonal(T([U, U]))
-        @test C.L ≈ BlockDiagonal(T([U', U']))
+        @test C.U ≈ BlockDiagonal(V([U, U]))
+        @test C.L ≈ BlockDiagonal(V([U', U']))
         @test C.UL ≈ C.U
         @test C.uplo === 'U'
         @test C.info == 0
@@ -181,8 +173,8 @@ using Test
         M = BlockDiagonal(map(Matrix, blocks(C.L)))
         C = Cholesky(M, 'L', 0)
         @test C.U ≈ cholesky(Matrix(B)).U
-        @test C.U ≈ BlockDiagonal(T([U, U]))
-        @test C.L ≈ BlockDiagonal(T([U', U']))
+        @test C.U ≈ BlockDiagonal(V([U, U]))
+        @test C.L ≈ BlockDiagonal(V([U', U']))
         @test C.UL ≈ C.L
         @test C.uplo === 'L'
         @test C.info == 0
@@ -192,7 +184,7 @@ using Test
         X = [  4  12 -16
               12  37 -43
              -16 -43  98]
-        B = BlockDiagonal(T([X, X]))
+        B = BlockDiagonal(V([X, X]))
 
         @testset "full=$full" for full in (true, false)
 

From c9bf91590f74b07309a7a34be458febc2b4f22a6 Mon Sep 17 00:00:00 2001
From: Nick Robinson <nicholas.robinson@invenialabs.co.uk>
Date: Mon, 15 Feb 2021 19:51:44 +0000
Subject: [PATCH 4/5] Improve tests in blockdiagonal.jl

---
 test/blockdiagonal.jl | 38 +++++++++++++++++---------------------
 1 file changed, 17 insertions(+), 21 deletions(-)

diff --git a/test/blockdiagonal.jl b/test/blockdiagonal.jl
index 16963ac..7373f20 100644
--- a/test/blockdiagonal.jl
+++ b/test/blockdiagonal.jl
@@ -11,18 +11,14 @@ using Test
     blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]
     blocks3 = [rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)]
 
-    @testset "$T" for (T, (b1, b2, b3)) in (
-        Tuple => (BlockDiagonal(Tuple(blocks1)), BlockDiagonal(Tuple(blocks2)), BlockDiagonal(Tuple(blocks3))),
-        Vector => (BlockDiagonal(blocks1), BlockDiagonal(blocks2), BlockDiagonal(blocks3)),
-    )
-    A = rand(rng, N, N + N1)
-    B = rand(rng, N + N1, N + N2)
-    A′, B′ = A', B'
-    a = rand(rng, N)
-    b = rand(rng, N + N1)
+    @testset for V in (Tuple, Vector)
+    b1 = BlockDiagonal(V(blocks1))
+    b2 = BlockDiagonal(V(blocks2))
+    N = size(b1, 1)
 
     @testset "AbstractArray" begin
-        X = rand(2, 2); Y = rand(3, 3)
+        X = rand(2, 2)
+        Y = rand(3, 3)
 
         @test size(b1) == (N, N)
         @test size(b1, 1) == N && size(b1, 2) == N
@@ -47,7 +43,7 @@ using Test
         end
 
         @testset "parent" begin
-            @test parent(b1) isa Union{Tuple,AbstractVector}
+            @test parent(b1) isa V
             @test eltype(parent(b1)) <: AbstractMatrix
             @test parent(BlockDiagonal([X, Y])) == [X, Y]
             @test parent(BlockDiagonal((X, Y))) == (X, Y)
@@ -60,7 +56,7 @@ using Test
         end
 
         @testset "setindex!" begin
-            X = BlockDiagonal([rand(Float32, 5, 5), rand(Float32, 3, 3)])
+            X = BlockDiagonal(V([rand(Float32, 5, 5), rand(Float32, 3, 3)]))
             X[10] = Int(10)
             @test X[10] === Float32(10.0)
             X[3, 3] = Int(9)
@@ -77,9 +73,9 @@ using Test
     end
 
     @testset "blocks size" begin
-        B = BlockDiagonal([rand(3, 3), rand(4, 4)])
+        B = BlockDiagonal(V([rand(3, 3), rand(4, 4)]))
         @test nblocks(B) == 2
-        @test blocksizes(B) == [(3, 3), (4, 4)]
+        @test blocksizes(B) == V([(3, 3), (4, 4)])
         @test blocksize(B, 2) == blocksizes(B)[2] == blocksize(B, 2, 2)
     end
 
@@ -101,15 +97,15 @@ using Test
     end  # Equality
 
     @testset "Non-Square Matrix" begin
-	A1 = ones(2, 4)
-	A2 = 2 * ones(3, 2)
-	B1 = BlockDiagonal([A1, A2])
+        A1 = ones(2, 4)
+        A2 = 2 * ones(3, 2)
+        B1 = BlockDiagonal(V([A1, A2]))
         B2 = [A1 zeros(2, 2); zeros(3, 4) A2]
 
-	@test B1 == B2
-	# Dimension check
-	@test sum(size.(B1.blocks, 1)) == size(B2, 1)
-	@test sum(size.(B1.blocks, 2)) == size(B2, 2)
+        @test B1 == B2
+        # Dimension check
+        @test sum(size.(B1.blocks, 1)) == size(B2, 1)
+        @test sum(size.(B1.blocks, 2)) == size(B2, 2)
     end  # Non-Square Matrix
 
     @testset "copy" begin

From d5d29da80fa335b623ecb03d9fd3e0284e6aba15 Mon Sep 17 00:00:00 2001
From: Nick Robinson <nicholas.robinson@invenialabs.co.uk>
Date: Mon, 15 Feb 2021 23:11:08 +0000
Subject: [PATCH 5/5] Update docs

---
 README.md            |  2 +-
 docs/Manifest.toml   | 63 +++++++++++++++++++++++++++++++++-----------
 docs/Project.toml    |  2 +-
 docs/make.jl         | 10 ++++---
 docs/src/index.md    |  5 ++--
 src/blockdiagonal.jl |  4 +--
 6 files changed, 60 insertions(+), 26 deletions(-)

diff --git a/README.md b/README.md
index f66d43c..3a0e946 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/dev)
-[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow%3ACI)
+[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow:CI)
 [![Codecov](https://codecov.io/gh/invenia/BlockDiagonals.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/invenia/BlockDiagonals.jl)
 [![code style blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)
 
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 4a20187..0a85124 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -4,36 +4,58 @@
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[BlockDiagonals]]
-deps = ["FillArrays", "LinearAlgebra"]
+deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra"]
 path = ".."
 uuid = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
-version = "0.1.0"
+version = "0.2.0"
+
+[[ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "de4f08843c332d355852721adb1592bce7924da3"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "0.9.29"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "919c7f3151e79ff196add81d7f4e45d91bbf420b"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.25.0"
 
 [[Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[DocStringExtensions]]
 deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "0513f1a8991e9d83255e0140aace0d0fc4486600"
+git-tree-sha1 = "50ddf44c53698f5e784bbebb3f4b21c5807401b1"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.0"
+version = "0.8.3"
 
 [[Documenter]]
-deps = ["Base64", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "4a84478277020abfff208cde31ba1aa68a5bc572"
+deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "21fb992ef1b28ff8f315354d3808ebf4a8fa6e45"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.23.0"
+version = "0.26.2"
 
 [[FillArrays]]
-deps = ["LinearAlgebra", "Random", "SparseArrays", "Test"]
-git-tree-sha1 = "9ab8f76758cbabba8d7f103c51dce7f73fcf8e92"
+deps = ["LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "ff537e5a3cba92fb48f30fec46723510450f2c0e"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "0.6.3"
+version = "0.10.2"
+
+[[IOCapture]]
+deps = ["Logging"]
+git-tree-sha1 = "377252859f740c217b936cebcd918a44f9b53b59"
+uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
+version = "0.1.1"
 
 [[InteractiveUtils]]
 deps = ["Markdown"]
@@ -41,11 +63,12 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
+git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.0"
+version = "0.21.1"
 
 [[LibGit2]]
+deps = ["Printf"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
 [[Libdl]]
@@ -66,13 +89,13 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[Parsers]]
-deps = ["Dates", "Test"]
-git-tree-sha1 = "db2b35dedab3c0e46dc15996d170af07a5ab91c9"
+deps = ["Dates"]
+git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "0.3.6"
+version = "1.0.15"
 
 [[Pkg]]
-deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[Printf]]
@@ -93,6 +116,10 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
@@ -100,6 +127,10 @@ uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [[Test]]
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/docs/Project.toml b/docs/Project.toml
index 31cea1d..5629748 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,4 +3,4 @@ BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 
 [compat]
-Documenter = "~0.23"
+Documenter = "0.26"
diff --git a/docs/make.jl b/docs/make.jl
index 870e79a..6844029 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,16 +1,20 @@
-using Documenter, BlockDiagonals
+using BlockDiagonals
+using Documenter
 
 makedocs(;
     modules=[BlockDiagonals],
-    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
+    format=Documenter.HTML(prettyurls=false),
     pages=[
         "Home" => "index.md",
     ],
-    repo="https://github.com/invenia/BlockDiagonals.jl/blob/{commit}{path}#L{line}",
+    repo="https://github.com/invenia/BlockDiagonals.jl/blob/{commit}{path}#{line}",
     sitename="BlockDiagonals.jl",
     authors="Invenia Technical Computing",
+    strict=true,
+    checkdocs=:exports,
 )
 
 deploydocs(;
     repo="github.com/invenia/BlockDiagonals.jl",
+    push_preview=true,
 )
diff --git a/docs/src/index.md b/docs/src/index.md
index 6fffd1a..aa1fcfa 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,11 +1,10 @@
 # BlockDiagonals.jl
 
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://invenia.github.io/BlockDiagonals.jl/dev)
-![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)
+[![CI](https://github.com/invenia/BlockDiagonals.jl/workflows/CI/badge.svg)](https://github.com/Invenia/BlockDiagonals.jl/actions?query=workflow:CI)
 [![Codecov](https://codecov.io/gh/invenia/BlockDiagonals.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/invenia/BlockDiagonals.jl)
 
-Functionality for working efficiently with [block diagonal matrices](https://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices).
+[BlockDiagonals.jl](https://github.com/invenia/BlockDiagonals.jl) provides functionality for working efficiently with [block diagonal matrices](https://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices).
 
 ```@autodocs
 Modules = [BlockDiagonals]
diff --git a/src/blockdiagonal.jl b/src/blockdiagonal.jl
index f974700..9842dfb 100644
--- a/src/blockdiagonal.jl
+++ b/src/blockdiagonal.jl
@@ -25,7 +25,7 @@ BlockDiagonal(B::BlockDiagonal) = B
 is_square(A::AbstractMatrix) = size(A, 1) == size(A, 2)
 
 """
-    blocks(B::BlockDiagonal{T, V}) -> Vector{V}
+    blocks(B::BlockDiagonal{T, V}) -> V
 
 Return the on-diagonal blocks of B.
 """
@@ -33,7 +33,7 @@ blocks(B::BlockDiagonal) = B.blocks
 
 # BlockArrays-like functions
 """
-    blocksizes(B::BlockDiagonal) -> Vector{Tuple}
+    blocksizes(B::BlockDiagonal{T, V}) -> V
 
 Return the size of each on-diagonal block in order.