Merge branch 'master' into patch-1

Project.toml

@@ -1,12 +1,14 @@
 name = "ForwardDiff"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.16"
+version = "0.10.18"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -26,9 +28,8 @@ julia = "1"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 DiffTests = "de460e47-3fe3-5279-bb4a-814414816d5d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Calculus", "DiffTests", "LinearAlgebra", "SparseArrays", "Test", "InteractiveUtils"]
+test = ["Calculus", "DiffTests", "SparseArrays", "Test", "InteractiveUtils"]

src/ForwardDiff.jl

@@ -4,7 +4,9 @@ using DiffRules, DiffResults
 using DiffResults: DiffResult, MutableDiffResult, ImmutableDiffResult
 using StaticArrays
 using Random
+using LinearAlgebra
 
+import Printf
 import NaNMath
 import SpecialFunctions
 import CommonSubexpressions

src/dual.jl

@@ -34,7 +34,7 @@ Base.showerror(io::IO, e::DualMismatchError{A,B}) where {A,B} =
 
 @noinline function throw_cannot_dual(V::Type)
     throw(ArgumentError("Cannot create a dual over scalar type $V." *
-        " If the type behaves as a scalar, define FowardDiff.can_dual."))
+        " If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{$V}) = true."))
 end
 
 """
@@ -346,7 +346,7 @@ end
 # Promotion/Conversion #
 ########################
 
-Base.@pure function Base.promote_rule(::Type{Dual{T1,V1,N1}},
+function Base.promote_rule(::Type{Dual{T1,V1,N1}},
                                       ::Type{Dual{T2,V2,N2}}) where {T1,V1,N1,T2,V2,N2}
     # V1 and V2 might themselves be Dual types
     if T2 ≺ T1
@@ -634,6 +634,52 @@ if VERSION >= v"1.6.0-DEV.292"
     end
 end
 
+# Symmetric eigvals #
+#-------------------#
+
+function LinearAlgebra.eigvals(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ,Q = eigen(Symmetric(value.(parent(A))))
+    parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N)
+    Dual{Tg}.(λ, tuple.(parts...))
+end
+
+function LinearAlgebra.eigvals(A::Hermitian{<:Complex{<:Dual{Tg,T,N}}}) where {Tg,T<:Real,N}
+    λ,Q = eigen(Hermitian(value.(real.(parent(A))) .+ im .* value.(imag.(parent(A)))))
+    parts = ntuple(j -> diag(real.(Q' * (getindex.(partials.(real.(A)) .+ im .* partials.(imag.(A)), j)) * Q)), N)
+    Dual{Tg}.(λ, tuple.(parts...))
+end
+
+function LinearAlgebra.eigvals(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev)))
+    parts = ntuple(j -> diag(Q' * getindex.(partials.(A), j) * Q), N)
+    Dual{Tg}.(λ, tuple.(parts...))
+end
+
+# A ./ (λ - λ') but with diag special cased
+function _lyap_div!(A, λ)
+    for (j,μ) in enumerate(λ), (k,λ) in enumerate(λ)
+        if k ≠ j
+            A[k,j] /= μ - λ
+        end
+    end
+    A
+end
+
+function LinearAlgebra.eigen(A::Symmetric{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ = eigvals(A)
+    _,Q = eigen(SymTridiagonal(value.(parent(A).dv),value.(parent(A).ev)))
+    parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
+    Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
+end
+
+function LinearAlgebra.eigen(A::SymTridiagonal{<:Dual{Tg,T,N}}) where {Tg,T<:Real,N}
+    λ = eigvals(A)
+    _,Q = eigen(SymTridiagonal(value.(parent(A))))
+    parts = ntuple(j -> Q*_lyap_div!(Q' * getindex.(partials.(A), j) * Q - Diagonal(getindex.(partials.(λ), j)), value.(λ)), N)
+    Eigen(λ,Dual{Tg}.(Q, tuple.(parts...)))
+end
+
+
 ###################
 # Pretty Printing #
 ###################
@@ -653,3 +699,7 @@ end
 function Base.typemax(::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
     ForwardDiff.Dual{T,V,N}(typemax(V))
 end
+
+if VERSION >= v"1.6.0-rc1"
+    Printf.tofloat(d::Dual) = Printf.tofloat(value(d))
+end

test/DualTest.jl

@@ -1,6 +1,7 @@
 module DualTest
 
 using Test
+using Printf
 using Random
 using ForwardDiff
 using ForwardDiff: Partials, Dual, value, partials
@@ -516,5 +517,15 @@ end
     @test length(UnitRange(Dual(1.5), Dual(3.5))) == 3
     @test length(UnitRange(Dual(1.5,1), Dual(3.5,3))) == 3
 end
+
+if VERSION >= v"1.6.0-rc1"
+    @testset "@printf" begin
+        for T in (Float16, Float32, Float64, BigFloat)
+            d1 = Dual(one(T))
+            @test_nowarn @printf("Testing @printf: %.2e\n", d1)
+            @test @sprintf("Testing @sprintf: %.2e\n", d1) == "Testing @sprintf: 1.00e+00\n"
+        end
+    end
+end
 
 end # module

test/JacobianTest.jl

@@ -7,6 +7,7 @@ using ForwardDiff
 using ForwardDiff: Dual, Tag, JacobianConfig
 using StaticArrays
 using DiffTests
+using LinearAlgebra
 
 include(joinpath(dirname(@__FILE__), "utils.jl"))
 
@@ -231,4 +232,9 @@ end
     @test_throws DimensionMismatch ForwardDiff.jacobian(sum, fill(2pi, 10^6)) # chunk_mode_jacobian
 end
 
+@testset "eigen" begin
+    @test ForwardDiff.jacobian(x -> eigvals(SymTridiagonal(x, x[1:end-1])), [1.,2.]) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2]
+    @test ForwardDiff.jacobian(x -> eigvals(Symmetric(x*x')), [1.,2.]) ≈ [0 0; 2 4]
+end
+
 end # module