diff --git a/src/structuredbroadcast.jl b/src/structuredbroadcast.jl index f2cfc74a..a8eeaa94 100644 --- a/src/structuredbroadcast.jl +++ b/src/structuredbroadcast.jl @@ -84,9 +84,9 @@ function structured_broadcast_alloc(bc, ::Type{Bidiagonal}, ::Type{ElType}, n) w return Bidiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n1), uplo) end structured_broadcast_alloc(bc, ::Type{SymTridiagonal}, ::Type{ElType}, n) where {ElType} = - SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, n-1)) + SymTridiagonal(Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1))) structured_broadcast_alloc(bc, ::Type{Tridiagonal}, ::Type{ElType}, n) where {ElType} = - Tridiagonal(Array{ElType}(undef, n-1),Array{ElType}(undef, n),Array{ElType}(undef, n-1)) + Tridiagonal(Array{ElType}(undef, max(0,n-1)),Array{ElType}(undef, n),Array{ElType}(undef, max(0,n-1))) structured_broadcast_alloc(bc, ::Type{LowerTriangular}, ::Type{ElType}, n) where {ElType} = LowerTriangular(Array{ElType}(undef, n, n)) structured_broadcast_alloc(bc, ::Type{UpperTriangular}, ::Type{ElType}, n) where {ElType} = diff --git a/test/structuredbroadcast.jl b/test/structuredbroadcast.jl index b8b16a8d..398d0ff5 100644 --- a/test/structuredbroadcast.jl +++ b/test/structuredbroadcast.jl @@ -11,160 +11,164 @@ isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), using .Main.SizedArrays @testset "broadcast[!] over combinations of scalars, structured matrices, and dense vectors/matrices" begin - N = 10 - s = rand() - fV = rand(N) - fA = rand(N, N) - Z = copy(fA) - D = Diagonal(rand(N)) - B = Bidiagonal(rand(N), rand(N - 1), :U) - T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1)) - S = SymTridiagonal(rand(N), rand(N - 1)) - U = UpperTriangular(rand(N,N)) - L = LowerTriangular(rand(N,N)) - M = Matrix(rand(N,N)) - structuredarrays = (D, B, T, U, L, M, S) - fstructuredarrays = map(Array, structuredarrays) - for (X, fX) in zip(structuredarrays, fstructuredarrays) - @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX)) - @test broadcast!(sin, Z, X) == broadcast(sin, fX) - @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) - @test broadcast!(cos, Z, X) == broadcast(cos, fX) - @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX)) - @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) - @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) - @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) - @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) - @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) - - @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 - @test X .* 2.0 isa typeof(X) - @test X .* (2.0,) isa typeof(X) - @test isequal(X .* Inf, fX .* Inf) - - two = 2 - @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two - @test X .^ 2 isa typeof(X) - @test X .^ (2,) isa typeof(X) - @test X .^ two isa typeof(X) - @test X .^ 0 == fX .^ 0 - @test X .^ -1 == fX .^ -1 - - for (Y, fY) in zip(structuredarrays, fstructuredarrays) - @test broadcast(+, X, Y) == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y) == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + @testset for N in (0,1,2,10) # some edge cases, and a structured case + s = rand() + fV = rand(N) + fA = rand(N, N) + Z = copy(fA) + D = Diagonal(rand(N)) + B = Bidiagonal(rand(N), rand(max(0,N-1)), :U) + T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1))) + S = SymTridiagonal(rand(N), rand(max(0,N-1))) + U = UpperTriangular(rand(N,N)) + L = LowerTriangular(rand(N,N)) + M = Matrix(rand(N,N)) + structuredarrays = (D, B, T, U, L, M, S) + fstructuredarrays = map(Array, structuredarrays) + @testset "$(nameof(typeof(X)))" for (X, fX) in zip(structuredarrays, fstructuredarrays) + @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX)) + @test broadcast!(sin, Z, X) == broadcast(sin, fX) + @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) + @test broadcast!(cos, Z, X) == broadcast(cos, fX) + @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX)) + @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) + @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) + @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) + @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) + @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) + + @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 + @test X .* 2.0 isa typeof(X) + @test X .* (2.0,) isa typeof(X) + @test isequal(X .* Inf, fX .* Inf) + + two = 2 + @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two + @test X .^ 2 isa typeof(X) + @test X .^ (2,) isa typeof(X) + @test X .^ two isa typeof(X) + @test X .^ 0 == fX .^ 0 + @test X .^ -1 == fX .^ -1 + + for (Y, fY) in zip(structuredarrays, fstructuredarrays) + @test broadcast(+, X, Y) == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y) == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - diagonals = (D, B, T) - fdiagonals = map(Array, diagonals) - for (X, fX) in zip(diagonals, fdiagonals) - for (Y, fY) in zip(diagonals, fdiagonals) - @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + diagonals = (D, B, T) + fdiagonals = map(Array, diagonals) + for (X, fX) in zip(diagonals, fdiagonals) + for (Y, fY) in zip(diagonals, fdiagonals) + @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - UU = UnitUpperTriangular(rand(N,N)) - UL = UnitLowerTriangular(rand(N,N)) - unittriangulars = (UU, UL) - Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU)))) - funittriangulars = map(Array, unittriangulars) - for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris) - @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX)) - @test broadcast!(sin, Z, X) == broadcast(sin, fX) - @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) - @test broadcast!(cos, Z, X) == broadcast(cos, fX) - @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX)) - @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) - @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) - @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) - @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) - @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) - - @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 - @test X .* 2.0 isa Ttri - @test X .* (2.0,) isa Ttri - @test isequal(X .* Inf, fX .* Inf) - - two = 2 - @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two - @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving - @test X .^ (2,) isa Ttri - @test X .^ two isa Ttri - @test X .^ 0 == fX .^ 0 - @test X .^ -1 == fX .^ -1 - - for (Y, fY) in zip(unittriangulars, funittriangulars) - @test broadcast(+, X, Y) == broadcast(+, fX, fY) - @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) - @test broadcast(*, X, Y) == broadcast(*, fX, fY) - @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + UU = UnitUpperTriangular(rand(N,N)) + UL = UnitLowerTriangular(rand(N,N)) + unittriangulars = (UU, UL) + Ttris = typeof.((UpperTriangular(parent(UU)), LowerTriangular(parent(UU)))) + funittriangulars = map(Array, unittriangulars) + for (X, fX, Ttri) in zip(unittriangulars, funittriangulars, Ttris) + @test (Q = broadcast(sin, X); typeof(Q) == Ttri && Q == broadcast(sin, fX)) + @test broadcast!(sin, Z, X) == broadcast(sin, fX) + @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX)) + @test broadcast!(cos, Z, X) == broadcast(cos, fX) + @test (Q = broadcast(*, s, X); typeof(Q) == Ttri && Q == broadcast(*, s, fX)) + @test broadcast!(*, Z, s, X) == broadcast(*, s, fX) + @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX)) + @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX) + @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX)) + @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX) + + @test X .* 2.0 == X .* (2.0,) == fX .* 2.0 + @test X .* 2.0 isa Ttri + @test X .* (2.0,) isa Ttri + @test isequal(X .* Inf, fX .* Inf) + + two = 2 + @test X .^ 2 == X .^ (2,) == fX .^ 2 == X .^ two + @test X .^ 2 isa typeof(X) # special cased, as isstructurepreserving + @test X .^ (2,) isa Ttri + @test X .^ two isa Ttri + @test X .^ 0 == fX .^ 0 + @test X .^ -1 == fX .^ -1 + + for (Y, fY) in zip(unittriangulars, funittriangulars) + @test broadcast(+, X, Y) == broadcast(+, fX, fY) + @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY) + @test broadcast(*, X, Y) == broadcast(*, fX, fY) + @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY) + end end - end - @testset "type-stability in Bidiagonal" begin - B2 = @inferred (B -> .- B)(B) - @test B2 isa Bidiagonal - @test B2 == -1 * B - B2 = @inferred (B -> B .* 2)(B) - @test B2 isa Bidiagonal - @test B2 == B + B - B2 = @inferred (B -> 2 .* B)(B) - @test B2 isa Bidiagonal - @test B2 == B + B - B2 = @inferred (B -> B ./ 1)(B) - @test B2 isa Bidiagonal - @test B2 == B - B2 = @inferred (B -> 1 .\ B)(B) - @test B2 isa Bidiagonal - @test B2 == B + @testset "type-stability in Bidiagonal" begin + B2 = @inferred (B -> .- B)(B) + @test B2 isa Bidiagonal + @test B2 == -1 * B + B2 = @inferred (B -> B .* 2)(B) + @test B2 isa Bidiagonal + @test B2 == B + B + B2 = @inferred (B -> 2 .* B)(B) + @test B2 isa Bidiagonal + @test B2 == B + B + B2 = @inferred (B -> B ./ 1)(B) + @test B2 isa Bidiagonal + @test B2 == B + B2 = @inferred (B -> 1 .\ B)(B) + @test B2 isa Bidiagonal + @test B2 == B + end end end @testset "broadcast! where the destination is a structured matrix" begin - N = 5 - A = rand(N, N) - sA = A + copy(A') - D = Diagonal(rand(N)) - Bu = Bidiagonal(rand(N), rand(N - 1), :U) - Bl = Bidiagonal(rand(N), rand(N - 1), :L) - T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1)) - ◣ = LowerTriangular(rand(N,N)) - ◥ = UpperTriangular(rand(N,N)) - M = Matrix(rand(N,N)) - - @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D)) - @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) - @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) - @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T)) - @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣)) - @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥)) - @test broadcast!(sin, copy(M), M) == Matrix(sin.(M)) - @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A)) - @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) - @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) - @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) - @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) - @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) - @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A)) - - @test_throws ArgumentError broadcast!(cos, copy(D), D) == Diagonal(sin.(D)) - @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) - @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) - @test_throws ArgumentError broadcast!(cos, copy(T), T) == Tridiagonal(sin.(T)) - @test_throws ArgumentError broadcast!(cos, copy(◣), ◣) == LowerTriangular(sin.(◣)) - @test_throws ArgumentError broadcast!(cos, copy(◥), ◥) == UpperTriangular(sin.(◥)) - @test_throws ArgumentError broadcast!(+, copy(D), D, A) == Diagonal(broadcast(*, D, A)) - @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) - @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) - @test_throws ArgumentError broadcast!(+, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) - @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) - @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) - @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2) - @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2) + @testset for N in (0,1,2,5) + A = rand(N, N) + sA = A + copy(A') + D = Diagonal(rand(N)) + Bu = Bidiagonal(rand(N), rand(max(0,N-1)), :U) + Bl = Bidiagonal(rand(N), rand(max(0,N-1)), :L) + T = Tridiagonal(rand(max(0,N-1)), rand(N), rand(max(0,N-1))) + ◣ = LowerTriangular(rand(N,N)) + ◥ = UpperTriangular(rand(N,N)) + M = Matrix(rand(N,N)) + + @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D)) + @test broadcast!(sin, copy(Bu), Bu) == Bidiagonal(sin.(Bu), :U) + @test broadcast!(sin, copy(Bl), Bl) == Bidiagonal(sin.(Bl), :L) + @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T)) + @test broadcast!(sin, copy(◣), ◣) == LowerTriangular(sin.(◣)) + @test broadcast!(sin, copy(◥), ◥) == UpperTriangular(sin.(◥)) + @test broadcast!(sin, copy(M), M) == Matrix(sin.(M)) + @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A)) + @test broadcast!(*, copy(Bu), Bu, A) == Bidiagonal(broadcast(*, Bu, A), :U) + @test broadcast!(*, copy(Bl), Bl, A) == Bidiagonal(broadcast(*, Bl, A), :L) + @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A)) + @test broadcast!(*, copy(◣), ◣, A) == LowerTriangular(broadcast(*, ◣, A)) + @test broadcast!(*, copy(◥), ◥, A) == UpperTriangular(broadcast(*, ◥, A)) + @test broadcast!(*, copy(M), M, A) == Matrix(broadcast(*, M, A)) + + if N > 2 + @test_throws ArgumentError broadcast!(cos, copy(D), D) + @test_throws ArgumentError broadcast!(cos, copy(Bu), Bu) + @test_throws ArgumentError broadcast!(cos, copy(Bl), Bl) + @test_throws ArgumentError broadcast!(cos, copy(T), T) + @test_throws ArgumentError broadcast!(cos, copy(◣), ◣) + @test_throws ArgumentError broadcast!(cos, copy(◥), ◥) + @test_throws ArgumentError broadcast!(+, copy(D), D, A) + @test_throws ArgumentError broadcast!(+, copy(Bu), Bu, A) + @test_throws ArgumentError broadcast!(+, copy(Bl), Bl, A) + @test_throws ArgumentError broadcast!(+, copy(T), T, A) + @test_throws ArgumentError broadcast!(+, copy(◣), ◣, A) + @test_throws ArgumentError broadcast!(+, copy(◥), ◥, A) + @test_throws ArgumentError broadcast!(*, copy(◥), ◣, 2) + @test_throws ArgumentError broadcast!(*, copy(Bu), Bl, 2) + end + end end @testset "map[!] over combinations of structured matrices" begin