diff --git a/src/dual.jl b/src/dual.jl index 6d13dec3..097f1d90 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -784,17 +784,20 @@ end #------------------------------------------------# # Extract structured matrices of primal values and partials -_structured_value(A::Symmetric{Dual{T,V,N}}) where {T,V,N} = Symmetric(map(value, parent(A)), A.uplo === 'U' ? :U : :L) -_structured_value(A::Hermitian{Dual{T,V,N}}) where {T,V,N} = Hermitian(map(value, parent(A)), A.uplo === 'U' ? :U : :L) -_structured_value(A::Hermitian{Complex{Dual{T,V,N}}}) where {T,V,N} = Hermitian(map(z -> splat(complex)(map(value, reim(z))), parent(A)), A.uplo === 'U' ? :U : :L) -_structured_value(A::SymTridiagonal{Dual{T,V,N}}) where {T,V,N} = SymTridiagonal(map(value, A.dv), map(value, A.ev)) +# non-isbits storage may be #undef outside of the `uplo` triangle +_maptri(f, A::Union{Symmetric,Hermitian}) = isbitstype(eltype(A)) ? map(f, parent(A)) : broadcast(f, A) -_structured_partials(A::Symmetric{Dual{T,V,N}}, j::Int) where {T,V,N} = Symmetric(partials.(parent(A), j), A.uplo === 'U' ? :U : :L) -_structured_partials(A::Hermitian{Dual{T,V,N}}, j::Int) where {T,V,N} = Hermitian(partials.(parent(A), j), A.uplo === 'U' ? :U : :L) +_structured_value(A::Symmetric{Dual{T,V,N}}) where {T,V,N} = Symmetric(_maptri(value, A), A.uplo === 'U' ? :U : :L) +_structured_value(A::Hermitian{Dual{T,V,N}}) where {T,V,N} = Hermitian(_maptri(value, A), A.uplo === 'U' ? :U : :L) +_structured_value(A::Hermitian{Complex{Dual{T,V,N}}}) where {T,V,N} = Hermitian(_maptri(z -> complex(value(real(z)), value(imag(z))), A), A.uplo === 'U' ? :U : :L) +_structured_value(A::SymTridiagonal{Dual{T,V,N}}) where {T,V,N} = SymTridiagonal(map(value, A.dv), map(value, view(A.ev, 1:length(A.dv)-1))) + +_structured_partials(A::Symmetric{Dual{T,V,N}}, j::Int) where {T,V,N} = Symmetric(_maptri(a -> partials(a, j), A), A.uplo === 'U' ? :U : :L) +_structured_partials(A::Hermitian{Dual{T,V,N}}, j::Int) where {T,V,N} = Hermitian(_maptri(a -> partials(a, j), A), A.uplo === 'U' ? :U : :L) function _structured_partials(A::Hermitian{Complex{Dual{T,V,N}}}, j::Int) where {T,V,N} - return Hermitian(complex.(partials.(real.(parent(A)), j), partials.(imag.(parent(A)), j)), A.uplo === 'U' ? :U : :L) + return Hermitian(_maptri(z -> complex(partials(real(z), j), partials(imag(z), j)), A), A.uplo === 'U' ? :U : :L) end -_structured_partials(A::SymTridiagonal{Dual{T,V,N}}, j::Int) where {T,V,N} = SymTridiagonal(partials.(A.dv, j), partials.(A.ev, j)) +_structured_partials(A::SymTridiagonal{Dual{T,V,N}}, j::Int) where {T,V,N} = SymTridiagonal(partials.(A.dv, j), partials.(view(A.ev, 1:length(A.dv)-1), j)) # Convert arrays of primal values and partials to arrays of Duals function _to_duals(::Val{T}, values::AbstractArray{<:Real}, partials::Tuple{Vararg{AbstractArray{<:Real}}}) where {T} diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 5050e23b..ab3be9cf 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -334,6 +334,35 @@ end @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) end end + + @testset "#undef outside of the `uplo` triangle" begin + x = ForwardDiff.Dual{Nothing}.(BigFloat[1, 2, 3], BigFloat[1, 0, 0], BigFloat[0, 1, 0]) + k = uplo === :U ? 3 : 2 # linear index of the stored off-diagonal element + M = similar(x, 2, 2) + M[1, 1], M[k], M[2, 2] = x[1], x[2], x[3] + Mc = similar(x, Complex{eltype(x)}, 2, 2) + Mc[1, 1], Mc[k], Mc[2, 2] = x[1], x[2] + im * x[1], x[3] + s = uplo === :U ? 1 : -1 # sign of imag(A[1, 2]) + @testset "$name" for (name, A, value, ∂1, ∂2) in ( + ("Symmetric{<:Real}", Symmetric(M, uplo), [1 2; 2 3], [1 0; 0 0], [0 1; 1 0]), + ("Hermitian{<:Real}", Hermitian(M, uplo), [1 2; 2 3], [1 0; 0 0], [0 1; 1 0]), + ("Hermitian{<:Complex}", Hermitian(Mc, uplo), [1 2+s*im; 2-s*im 3], [1 s*im; -s*im 0], [0 1; 1 0]), + ) + @test ForwardDiff._structured_value(A) == value + @test ForwardDiff._structured_partials(A, 1) == ∂1 + @test ForwardDiff._structured_partials(A, 2) == ∂2 + end + end + + @testset "#undef in the ignored element of SymTridiagonal" begin + dv = ForwardDiff.Dual{Nothing}.(BigFloat[1, 2, 3], BigFloat[1, 0, 0], BigFloat[0, 1, 0]) + ev = similar(dv, 3) + ev[1], ev[2] = dv[1], dv[2] + A = SymTridiagonal(dv, ev) + @test ForwardDiff._structured_value(A) == [1 1 0; 1 2 2; 0 2 3] + @test ForwardDiff._structured_partials(A, 1) == [1 1 0; 1 0 0; 0 0 0] + @test ForwardDiff._structured_partials(A, 2) == [0 0 0; 0 1 1; 0 1 0] + end end end