From 421c25f701d95a0e326d6e05e0180b8347409d0e Mon Sep 17 00:00:00 2001 From: araujoms Date: Mon, 21 Sep 2026 20:52:52 +0200 Subject: [PATCH 1/5] fix crash of eigen with #undef elements --- src/dual.jl | 12 +++++------ test/JacobianTest.jl | 48 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 6 deletions(-) diff --git a/src/dual.jl b/src/dual.jl index 6d13dec3..8a9eeb28 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -784,15 +784,15 @@ 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::Symmetric{Dual{T,V,N}}) where {T,V,N} = Symmetric(value.(A), A.uplo === 'U' ? :U : :L) +_structured_value(A::Hermitian{Dual{T,V,N}}) where {T,V,N} = Hermitian(value.(A), A.uplo === 'U' ? :U : :L) +_structured_value(A::Hermitian{Complex{Dual{T,V,N}}}) where {T,V,N} = Hermitian(broadcast(z -> splat(complex)(map(value, reim(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, A.ev)) -_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_partials(A::Symmetric{Dual{T,V,N}}, j::Int) where {T,V,N} = Symmetric(partials.(A, j), A.uplo === 'U' ? :U : :L) +_structured_partials(A::Hermitian{Dual{T,V,N}}, j::Int) where {T,V,N} = Hermitian(partials.(A, j), 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(complex.(partials.(real.(A), j), partials.(imag.(A), j)), 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)) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index 5050e23b..b3bd2b93 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -335,6 +335,54 @@ end end end end + + # the M[2, 1] element is undef, we're testing whether ForwardDiff tries to read from it + @testset "#undef" begin + function raw_real(x) + M = similar(x, 2, 2) + M[1, 1] = x[1] + M[1, 2] = x[2] + M[2, 2] = x[3] + return M + end + #mock function so we don't need to import GenericLinearAlgebra + function LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) + a, b, c = M[1], M[4], M[1, 2] + rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) + return [(a+b)/2-rt, (a+b)/2+rt] + end + function LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) + values = eigvals(M) + x = (values[1] - M[1])/M[1, 2] + vectors = [-1 -x; -x 1]/sqrt(1+x^2) + return Eigen(values, vectors) + end + LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(Symmetric(M)) + LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigen(Symmetric(M)) + function LinearAlgebra.eigvals(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) + a, b, c = M[1], M[4], M[1, 2] + rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) + return [(a+b)/2-rt, (a+b)/2+rt] + end + function LinearAlgebra.eigen(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) + values, vectors = eigen(Hermitian(real(Matrix(M)))) #workaround for a bug in julia 1.10 + return Eigen(values, complex(vectors)) + end + + x0 = BigFloat[1, 2, 3] + + for wrap in ( + (x -> Symmetric(raw_real(x))), + (x -> Hermitian(raw_real(x))), + (x -> Hermitian(complex(raw_real(x)))), + ) + for ev in (x -> eigvals(wrap(x)), + x -> eigen(wrap(x)).values, + x -> map(abs2, eigen(wrap(x)).vectors[:, 1])) + @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) + end + end + end end @testset "type stability" begin From cb10c4ea57c344306ef8544db67dcf980ae8a294 Mon Sep 17 00:00:00 2001 From: araujoms Date: Tue, 22 Sep 2026 11:51:46 +0200 Subject: [PATCH 2/5] address review comments --- src/dual.jl | 15 ++++++----- test/JacobianTest.jl | 63 ++++++++++++++++++++++++++------------------ 2 files changed, 46 insertions(+), 32 deletions(-) diff --git a/src/dual.jl b/src/dual.jl index 8a9eeb28..50f9f3ba 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -784,15 +784,18 @@ end #------------------------------------------------# # Extract structured matrices of primal values and partials -_structured_value(A::Symmetric{Dual{T,V,N}}) where {T,V,N} = Symmetric(value.(A), A.uplo === 'U' ? :U : :L) -_structured_value(A::Hermitian{Dual{T,V,N}}) where {T,V,N} = Hermitian(value.(A), A.uplo === 'U' ? :U : :L) -_structured_value(A::Hermitian{Complex{Dual{T,V,N}}}) where {T,V,N} = Hermitian(broadcast(z -> splat(complex)(map(value, reim(z))), A), A.uplo === 'U' ? :U : :L) +_maptri(f, A::Union{Symmetric,Hermitian}) = _maptri(f, A, parent(A)) +_maptri(f, A, P::AbstractArray{V}) where {V} = isbitstype(V) ? map(f, P) : broadcast(f, A) + +_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, A.ev)) -_structured_partials(A::Symmetric{Dual{T,V,N}}, j::Int) where {T,V,N} = Symmetric(partials.(A, j), A.uplo === 'U' ? :U : :L) -_structured_partials(A::Hermitian{Dual{T,V,N}}, j::Int) where {T,V,N} = Hermitian(partials.(A, j), A.uplo === 'U' ? :U : :L) +_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.(A), j), partials.(imag.(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)) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index b3bd2b93..fbbf6db7 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -338,48 +338,59 @@ end # the M[2, 1] element is undef, we're testing whether ForwardDiff tries to read from it @testset "#undef" begin - function raw_real(x) + function raw_real(x, uplo) + offdiagonal = uplo == :U ? 3 : 2 M = similar(x, 2, 2) M[1, 1] = x[1] - M[1, 2] = x[2] + M[offdiagonal] = x[2] M[2, 2] = x[3] return M end - #mock function so we don't need to import GenericLinearAlgebra - function LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) - a, b, c = M[1], M[4], M[1, 2] - rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) - return [(a+b)/2-rt, (a+b)/2+rt] + function raw_complex(x, uplo) + offdiagonal = uplo == :U ? 3 : 2 + M = complex(similar(x, 2, 2)) + M[1, 1] = x[1] + M[offdiagonal] = x[2] + im*x[2] + M[2, 2] = x[3] + return M end - function LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) - values = eigvals(M) - x = (values[1] - M[1])/M[1, 2] - vectors = [-1 -x; -x 1]/sqrt(1+x^2) - return Eigen(values, vectors) + + #mock function so we don't need to import GenericLinearAlgebra + LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigvals(Hermitian(M)) + LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigen(Hermitian(M)) + LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(complex(M)) + function LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) + values, vectors = eigen(complex(M)) + return Eigen(values, real(vectors)) end - LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(Symmetric(M)) - LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigen(Symmetric(M)) function LinearAlgebra.eigvals(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - a, b, c = M[1], M[4], M[1, 2] + a, b, c = real(M[1]), real(M[4]), M[1, 2] rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) return [(a+b)/2-rt, (a+b)/2+rt] end function LinearAlgebra.eigen(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - values, vectors = eigen(Hermitian(real(Matrix(M)))) #workaround for a bug in julia 1.10 - return Eigen(values, complex(vectors)) + values = eigvals(M) + x = (values[1] - M[1])/M[1, 2] + vectors = [conj(sign(x)) -conj(x); + abs(x) 1 ]/sqrt(1+abs2(x)) + return Eigen(values, vectors) end x0 = BigFloat[1, 2, 3] - for wrap in ( - (x -> Symmetric(raw_real(x))), - (x -> Hermitian(raw_real(x))), - (x -> Hermitian(complex(raw_real(x)))), - ) - for ev in (x -> eigvals(wrap(x)), - x -> eigen(wrap(x)).values, - x -> map(abs2, eigen(wrap(x)).vectors[:, 1])) - @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) + for uplo in (:U, :L) + for wrap in ( + (x -> Symmetric(raw_real(x, uplo), uplo)), + (x -> Hermitian(raw_real(x, uplo), uplo)), + (x -> Hermitian(raw_complex(x, uplo), uplo)), + ) + for ev in ( + x -> eigvals(wrap(x)), + x -> eigen(wrap(x)).values, + x -> map(abs2, eigen(wrap(x)).vectors[:, 1]) + ) + @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) + end end end end From 36e41e71a4cb5a9dca8b639998ad59b87c0e4e55 Mon Sep 17 00:00:00 2001 From: araujoms Date: Wed, 23 Sep 2026 12:54:41 +0200 Subject: [PATCH 3/5] address review comments --- src/dual.jl | 4 +-- test/JacobianTest.jl | 81 ++++++++++++++++---------------------------- 2 files changed, 31 insertions(+), 54 deletions(-) diff --git a/src/dual.jl b/src/dual.jl index 50f9f3ba..5fa3fe5a 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -784,8 +784,8 @@ end #------------------------------------------------# # Extract structured matrices of primal values and partials -_maptri(f, A::Union{Symmetric,Hermitian}) = _maptri(f, A, parent(A)) -_maptri(f, A, P::AbstractArray{V}) where {V} = isbitstype(V) ? map(f, P) : broadcast(f, A) +# 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_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) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index fbbf6db7..6fa14ec5 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -313,32 +313,31 @@ end end end - # The matrices above are all of the form `x*x'` and hence symmetric, so `:U` and `:L` wrap the - # same matrix. Here the raw storage is deliberately not symmetric/Hermitian: the values and the - # partials both have to be read from the triangle that `uplo` selects. - @testset "uplo = :$uplo" for uplo in (:U, :L) - # `2*x[1]+x[2]^2` rather than something proportional to the off-diagonal entry, so that the - # eigenvector direction actually depends on `x` and the eigenvector test is not vacuous - raw_real(x) = [x[1] x[2]; 3*x[2] 2*x[1]+x[2]^2] - raw_complex(x) = complex.(raw_real(x), [0 x[1]; -2*x[2] 0]) - x0 = [1.0, 2.0] - - @testset "$name" for (name, wrap) in ( - ("Symmetric{<:Real}", x -> Symmetric(raw_real(x), uplo)), - ("Hermitian{<:Real}", x -> Hermitian(raw_real(x), uplo)), - ("Hermitian{<:Complex}", x -> Hermitian(raw_complex(x), uplo)), - ) - for ev in (x -> eigvals(wrap(x)), - x -> eigen(wrap(x)).values, - x -> map(abs2, eigen(wrap(x)).vectors[:, 1])) - @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) - end - end + #mock functions for the following thest so we don't need to import GenericLinearAlgebra + + LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigvals(Hermitian(M)) + LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigen(Hermitian(M)) + LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(complex(M)) + function LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) + values, vectors = eigen(complex(M)) + return Eigen(values, real(vectors)) + end + function LinearAlgebra.eigvals(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) + a, b, c = real(M[1]), real(M[4]), M[1, 2] + rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) + return [(a+b)/2-rt, (a+b)/2+rt] + end + function LinearAlgebra.eigen(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) + values = eigvals(M) + x = (values[1] - M[1])/M[1, 2] + vectors = [conj(sign(x)) -conj(x); + abs(x) 1 ]/sqrt(1+abs2(x)) + return Eigen(values, vectors) end - # the M[2, 1] element is undef, we're testing whether ForwardDiff tries to read from it - @testset "#undef" begin - function raw_real(x, uplo) + # one of the off-diagonal elements is undef, we're testing whether ForwardDiff tries to read from it + @testset "#undef uplo = :$uplo" for uplo in (:U, :L) + function raw_real(x) offdiagonal = uplo == :U ? 3 : 2 M = similar(x, 2, 2) M[1, 1] = x[1] @@ -346,43 +345,21 @@ end M[2, 2] = x[3] return M end - function raw_complex(x, uplo) + function raw_complex(x) offdiagonal = uplo == :U ? 3 : 2 M = complex(similar(x, 2, 2)) M[1, 1] = x[1] - M[offdiagonal] = x[2] + im*x[2] + M[offdiagonal] = x[2] + im*x[1] M[2, 2] = x[3] return M end - #mock function so we don't need to import GenericLinearAlgebra - LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigvals(Hermitian(M)) - LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigen(Hermitian(M)) - LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(complex(M)) - function LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) - values, vectors = eigen(complex(M)) - return Eigen(values, real(vectors)) - end - function LinearAlgebra.eigvals(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - a, b, c = real(M[1]), real(M[4]), M[1, 2] - rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) - return [(a+b)/2-rt, (a+b)/2+rt] - end - function LinearAlgebra.eigen(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - values = eigvals(M) - x = (values[1] - M[1])/M[1, 2] - vectors = [conj(sign(x)) -conj(x); - abs(x) 1 ]/sqrt(1+abs2(x)) - return Eigen(values, vectors) - end - - x0 = BigFloat[1, 2, 3] - for uplo in (:U, :L) + for x0 in (BigFloat[1, 2, 3], Float64[1, 2, 3]) for wrap in ( - (x -> Symmetric(raw_real(x, uplo), uplo)), - (x -> Hermitian(raw_real(x, uplo), uplo)), - (x -> Hermitian(raw_complex(x, uplo), uplo)), + (x -> Symmetric(raw_real(x), uplo)), + (x -> Hermitian(raw_real(x), uplo)), + (x -> Hermitian(raw_complex(x), uplo)), ) for ev in ( x -> eigvals(wrap(x)), From f479f9850879e37716cbcc02989c9ffef310c9cb Mon Sep 17 00:00:00 2001 From: araujoms Date: Fri, 25 Sep 2026 14:46:43 +0200 Subject: [PATCH 4/5] address review comments --- src/dual.jl | 4 +- test/JacobianTest.jl | 87 ++++++++++++++++++-------------------------- 2 files changed, 37 insertions(+), 54 deletions(-) diff --git a/src/dual.jl b/src/dual.jl index 5fa3fe5a..097f1d90 100644 --- a/src/dual.jl +++ b/src/dual.jl @@ -790,14 +790,14 @@ _maptri(f, A::Union{Symmetric,Hermitian}) = isbitstype(eltype(A)) ? map(f, paren _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, A.ev)) +_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(_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 6fa14ec5..e2c21c25 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -313,61 +313,44 @@ end end end - #mock functions for the following thest so we don't need to import GenericLinearAlgebra - - LinearAlgebra.eigvals(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigvals(Hermitian(M)) - LinearAlgebra.eigen(M::Symmetric{BigFloat, Matrix{BigFloat}}) = eigen(Hermitian(M)) - LinearAlgebra.eigvals(M::Hermitian{BigFloat, Matrix{BigFloat}}) = eigvals(complex(M)) - function LinearAlgebra.eigen(M::Hermitian{BigFloat, Matrix{BigFloat}}) - values, vectors = eigen(complex(M)) - return Eigen(values, real(vectors)) - end - function LinearAlgebra.eigvals(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - a, b, c = real(M[1]), real(M[4]), M[1, 2] - rt = sqrt(((a+b)/2)^2 + abs2(c) - a * b) - return [(a+b)/2-rt, (a+b)/2+rt] - end - function LinearAlgebra.eigen(M::Hermitian{Complex{BigFloat}, Matrix{Complex{BigFloat}}}) - values = eigvals(M) - x = (values[1] - M[1])/M[1, 2] - vectors = [conj(sign(x)) -conj(x); - abs(x) 1 ]/sqrt(1+abs2(x)) - return Eigen(values, vectors) - end - - # one of the off-diagonal elements is undef, we're testing whether ForwardDiff tries to read from it - @testset "#undef uplo = :$uplo" for uplo in (:U, :L) - function raw_real(x) - offdiagonal = uplo == :U ? 3 : 2 - M = similar(x, 2, 2) - M[1, 1] = x[1] - M[offdiagonal] = x[2] - M[2, 2] = x[3] - return M - end - function raw_complex(x) - offdiagonal = uplo == :U ? 3 : 2 - M = complex(similar(x, 2, 2)) - M[1, 1] = x[1] - M[offdiagonal] = x[2] + im*x[1] - M[2, 2] = x[3] - return M + # The matrices above are all of the form `x*x'` and hence symmetric, so `:U` and `:L` wrap the + # same matrix. Here the raw storage is deliberately not symmetric/Hermitian: the values and the + # partials both have to be read from the triangle that `uplo` selects. + @testset "uplo = :$uplo" for uplo in (:U, :L) + # `2*x[1]+x[2]^2` rather than something proportional to the off-diagonal entry, so that the + # eigenvector direction actually depends on `x` and the eigenvector test is not vacuous + raw_real(x) = [x[1] x[2]; 3*x[2] 2*x[1]+x[2]^2] + raw_complex(x) = complex.(raw_real(x), [0 x[1]; -2*x[2] 0]) + x0 = [1.0, 2.0] + + @testset "$name" for (name, wrap) in ( + ("Symmetric{<:Real}", x -> Symmetric(raw_real(x), uplo)), + ("Hermitian{<:Real}", x -> Hermitian(raw_real(x), uplo)), + ("Hermitian{<:Complex}", x -> Hermitian(raw_complex(x), uplo)), + ) + for ev in (x -> eigvals(wrap(x)), + x -> eigen(wrap(x)).values, + x -> map(abs2, eigen(wrap(x)).vectors[:, 1])) + @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) + end end - - for x0 in (BigFloat[1, 2, 3], Float64[1, 2, 3]) - for wrap in ( - (x -> Symmetric(raw_real(x), uplo)), - (x -> Hermitian(raw_real(x), uplo)), - (x -> Hermitian(raw_complex(x), uplo)), + @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]), ) - for ev in ( - x -> eigvals(wrap(x)), - x -> eigen(wrap(x)).values, - x -> map(abs2, eigen(wrap(x)).vectors[:, 1]) - ) - @test ForwardDiff.jacobian(ev, x0) ≈ Calculus.finite_difference_jacobian(ev, x0) - end + @test ForwardDiff._structured_value(A) == value + @test ForwardDiff._structured_partials(A, 1) == ∂1 + @test ForwardDiff._structured_partials(A, 2) == ∂2 end end end From 3520531e49142e719e80dd9a878295133e640c0a Mon Sep 17 00:00:00 2001 From: araujoms Date: Fri, 25 Sep 2026 14:50:05 +0200 Subject: [PATCH 5/5] add SymTridiagonal test --- test/JacobianTest.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/JacobianTest.jl b/test/JacobianTest.jl index e2c21c25..ab3be9cf 100644 --- a/test/JacobianTest.jl +++ b/test/JacobianTest.jl @@ -353,6 +353,16 @@ end @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