Skip to content

RFC: Allows blocks to be <:AbstractVector or <:Tuple #56

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 5 commits into from
+151 −96
Closed

RFC: Allows blocks to be <:AbstractVector or <:Tuple #56

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
ac3ddc2
Allows blocks to be <:AbstractVector or <:Tuple
nickrobinson251 Feb 2, 2021
32f9d07
Simplify test setup for base_maths.jl
nickrobinson251 Feb 15, 2021
8b3235c
Simplify test setup for linalg.jl
nickrobinson251 Feb 15, 2021
c9bf915
Improve tests in blockdiagonal.jl
nickrobinson251 Feb 15, 2021
d5d29da
Update docs
nickrobinson251 Feb 15, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
23 changes: 14 additions & 9 deletions src/blockdiagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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`.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we want to constrain T to at least a Number, or are we happy with anything?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think if T isn't Number-like, some methods will fail. But 🤷 Would there be any gain from constraining things?

"""
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}}}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you really want all block matrices to have the same eltype? Or could you live with them having different eltypes, and T being promote_type(map(eltype, blocks)...)? See LinearMaps.jl. If you need them to be homogeneous, then you may wish to include some promotion mechanism.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could you live with them having different eltypes, and T being promote_type(map(eltype, blocks)...)?

that's a great suggestion

to be honest, i hadn't considered it. This currently just keeps T with the same meaning as on master.
Are there any downsides to T being promote_type(map(eltype, blocks)...)?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

trying this locally, and running tests, it does cause some inference issues for rrule

BlockDiagonal * Vector: Error During Test at /Users/npr/repos/BlockDiagonals.jl/test/chainrules.jl:12
  Got exception outside of a @test
  return type Array{Float64,1} does not match inferred return type Any
  Stacktrace:
   [1] error(::String) at ./error.jl:33
   [2] _test_inferred(::Function, ::InplaceableThunk{Thunk{BlockDiagonals.var"#13#19"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}},BlockDiagonals.var"#14#20"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:228
   [3] _test_inferred(::Function, ::InplaceableThunk{Thunk{BlockDiagonals.var"#13#19"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}},BlockDiagonals.var"#14#20"{Array{Float64,1},BlockDiagonal{Float64,Array{Array{Float64,2},1}}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:227
   [4] test_rrule(::typeof(*), ::BlockDiagonal{Float64,Array{Array{Float64,2},1}}, ::Vararg{Any,N} where N; output_tangent::ChainRulesTestUtils.Auto, fdm::FiniteDifferences.AdaptedFiniteDifferenceMethod{5,1,FiniteDifferences.UnadaptedFiniteDifferenceMethod{7,5}}, check_inferred::Bool, fkwargs::NamedTuple{(),Tuple{}}, rtol::Float64, atol::Float64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:187
   [5] test_rrule(::Function, ::BlockDiagonal{Float64,Array{Array{Float64,2},1}}, ::Vararg{Any,N} where N) at /Users/npr/.julia/packages/ChainRulesTestUtils/2QRer/src/testers.jl:153
   [6] top-level scope at /Users/npr/repos/BlockDiagonals.jl/test/chainrules.jl:15

perhaps someone more familiar with ChainRules has some insight on why this might be (@willtebbutt , @mzgubic)?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know why this is (I assume the answer contains the word compiler). But I'd rather live with the same eltype than with this I think

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, this change would be non-breaking anyway, so let's leave it as a potentially follow-up #62

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that's perfectly fine. But then you may wish to consider to include a constructor that by itself does not yet require equal eltypes, but first determines T via promotion and then map(B -> convert(AbstractMatrix{T}, B), blocks). This should be a no-op if all matrices have the same eltype, but otherwise relieve you from annoying method errors. Just a suggestion.

}
return BlockDiagonal{T, V}(blocks)
end

Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/chainrules.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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}
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think AbstractVector is fine here, but i'll rely on @willtebbutt to tell me if not

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we exclude the tuple we won't have an rrule for constructors of BlockDiagonals with tuple blocks. Are you asking if it's ok not to have rrule for tuple blocks, or am I missing the point of the question?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we won't have an rrule for constructors of BlockDiagonals with tuple blocks

correct, that'd be a nice feature to add as a follow-up, but out-of-scope for this PR

But my question is a different one. Before this PR, blocks were always a Vector, so we had rrule(::Type{<:BlockDiagonal}, blocks::Vector) whereas this PR allows them to be any AbstractVector (or Tuple but ignore that, since we're not adding that rrule here). So my question is, now that blocks can be any AbstractVector, are we okay with broadening the rrule from rrule(::Type{<:BlockDiagonal}, blocks::Vector) to rrule(::Type{<:BlockDiagonal}, blocks::AbstractVector) ?

(it's a question about how tightly to constrain rrules #56 (comment))

BlockDiagonal_pullback(Δ::Composite) = (NO_FIELDS, Δ.blocks)
return BlockDiagonal(blocks), BlockDiagonal_pullback
end
Expand Down Expand Up @@ -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}}}
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this okay, @willtebbutt ?


y = bm * v

Expand Down
13 changes: 12 additions & 1 deletion src/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)}[]
Expand All @@ -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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How come we aren't push!ing like above? Or vice versa if first.()/last.() is better?

Copy link
Contributor Author

@nickrobinson251 nickrobinson251 Feb 16, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

because above we want U (likewise Vt) to be a Vector, so create a Vector and add the elements. Whereas here we want it to be a Tuple hence ntuple, and first.(ntup) will still give us a tuple.

rather than push! we should probably be preallocating a Vector of the right length and overwriting the values, but the two functions would still be different (since Tuples are immutable)

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)
Expand Down
45 changes: 24 additions & 21 deletions test/base_maths.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,33 +5,31 @@ using Test

@testset "base_maths.jl" begin
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)])
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.
Expand All @@ -45,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
Expand All @@ -68,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
Expand All @@ -83,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)
Expand All @@ -97,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
Expand All @@ -114,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
Expand All @@ -127,8 +129,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(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)]
Expand All @@ -137,4 +139,5 @@ using Test
@test sum(size.(b5.blocks, 2)) == size(b4 * b5, 2)
end
end # Multiplication
end # V
end
44 changes: 24 additions & 20 deletions test/blockdiagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,17 +7,18 @@ 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)])
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, 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 for V in (Tuple, Vector)
b1 = BlockDiagonal(V(blocks1))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we shold introduce an indent under this new @testset for V in (Tuple, Vector) loop, but that'd give us huge diffs, and make it harder to see what's changed in the tests, so to make review easier i'll make that aesthetic change in a follow-up PR

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
Expand All @@ -42,8 +43,10 @@ using Test
end

@testset "parent" begin
@test parent(b1) isa Vector{<:AbstractMatrix}
@test parent(b1) isa V
@test eltype(parent(b1)) <: AbstractMatrix
@test parent(BlockDiagonal([X, Y])) == [X, Y]
@test parent(BlockDiagonal((X, Y))) == (X, Y)
end

@testset "similar" begin
Expand All @@ -53,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)
Expand All @@ -70,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

Expand All @@ -94,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
Expand All @@ -117,3 +120,4 @@ using Test
@test_throws DimensionMismatch copy!(b2, b1)
end
end
end
30 changes: 14 additions & 16 deletions test/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,12 @@ using Test
@testset "linalg.jl" begin
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)])
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, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)]
blocks2 = [rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)]

@testset for V in (Tuple, Vector)
b1 = BlockDiagonal(V(blocks1))
b2 = BlockDiagonal(V(blocks2))

@testset "mul!" begin
c = similar(b1)
Expand Down Expand Up @@ -143,7 +140,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(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
Expand All @@ -163,21 +160,21 @@ using Test
0.0 1.0 5.0
0.0 0.0 3.0]

B = BlockDiagonal([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([U, U])
@test C.L ≈ BlockDiagonal([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

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(V([U, U]))
@test C.L ≈ BlockDiagonal(V([U', U']))
@test C.UL ≈ C.L
@test C.uplo === 'L'
@test C.info == 0
Expand All @@ -187,7 +184,7 @@ using Test
X = [ 4 12 -16
12 37 -43
-16 -43 98]
B = BlockDiagonal([X, X])
B = BlockDiagonal(V([X, X]))

@testset "full=$full" for full in (true, false)

Expand Down Expand Up @@ -230,3 +227,4 @@ using Test
end
end # SVD
end
end