From 7ca98fb94f6be02105ab49613a673092b1cf481b Mon Sep 17 00:00:00 2001 From: Keita Nakamura Date: Sat, 31 Aug 2024 18:45:42 +0900 Subject: [PATCH] Use eigen functions in StaticArrays.jl (#211) * Use eigen functions in StaticArrays.jl * Set version to 0.15.0 --- Project.toml | 2 +- src/Tensorial.jl | 1 - src/eigen.jl | 97 ------------------------------------------------ src/ops.jl | 10 +++++ test/eigen.jl | 86 ------------------------------------------ test/ops.jl | 7 ++++ test/runtests.jl | 1 - 7 files changed, 18 insertions(+), 186 deletions(-) delete mode 100644 src/eigen.jl delete mode 100644 test/eigen.jl diff --git a/Project.toml b/Project.toml index 72ea104e..db7f4242 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Tensorial" uuid = "98f94333-fa9f-48a9-ad80-1c66397b2b38" authors = ["Keita Nakamura "] -version = "0.14.1" +version = "0.15.0" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" diff --git a/src/Tensorial.jl b/src/Tensorial.jl index b643f1ad..f4cd857d 100644 --- a/src/Tensorial.jl +++ b/src/Tensorial.jl @@ -91,7 +91,6 @@ include("permute.jl") include("ops.jl") include("continuum_mechanics.jl") include("inv.jl") -include("eigen.jl") include("voigt.jl") include("ad.jl") # include("simd.jl") diff --git a/src/eigen.jl b/src/eigen.jl deleted file mode 100644 index 94dd20f8..00000000 --- a/src/eigen.jl +++ /dev/null @@ -1,97 +0,0 @@ -@inline function eigvals(x::AbstractSymmetricSecondOrderTensor; permute::Bool=true, scale::Bool=true) - Tensor(eigvals(Symmetric(SArray(x)); permute=permute, scale = scale)) -end - -@inline function eigen(x::AbstractSymmetricSecondOrderTensor; permute::Bool=true, scale::Bool=true) - _eig(x; permute=permute, scale=scale) -end - -function eigen(x::AbstractSymmetricSecondOrderTensor{2, T}; permute::Bool=true, scale::Bool=true) where {T} - _iszero(x) = abs(x) < eps(typeof(x)) - a, c, b = Tuple(x) - if _iszero(c) - # check diagonal case - # `eigen(Symmetric(SMatrix{2,2}([1.0 1.9932760451045367e-17; 1.9932760451045367e-17 1.0])))` fails in StaticArrays.jl - return Eigen(Vec(a,b), one(Mat{2,2,T})) - else - return _eig(x; permute=permute, scale=scale) - end -end - -function eigen(x::AbstractSymmetricSecondOrderTensor{3, T}; permute::Bool=true, scale::Bool=true) where {T} - _isapproxzero(x) = abs(x) < cbrt(eps(typeof(x))) # sqrt fails, so use `cbrt` - _iszero(x) = abs(x) < eps(typeof(x)) - - a, d, f, b, e, c = Tuple(x) - - # block diagonal - if !_isapproxzero(a) && _iszero(d) && _iszero(f) - vals, vecs = eigen(@Tensor(x[[2,3], [2,3]]); permute=permute, scale=scale) - λ₁ = a - λ₂ = vals[1] - λ₃ = vals[2] - v₁ = Vec{3,T}(1,0,0) - v₂ = vcat(0, vecs[:,1]) - v₃ = vcat(0, vecs[:,2]) - return _eig33_construct((λ₁, λ₂, λ₃), (v₁, v₂, v₃), true) - elseif !_isapproxzero(b) && _iszero(d) && _iszero(e) - vals, vecs = eigen(@Tensor(x[[1,3], [1,3]]); permute=permute, scale=scale) - λ₁ = vals[1] - λ₂ = b - λ₃ = vals[2] - v₁ = Vec(vecs[1,1], 0, vecs[2,1]) - v₂ = Vec{3,T}(0,1,0) - v₃ = Vec(vecs[1,2], 0, vecs[2,2]) - return _eig33_construct((λ₁, λ₂, λ₃), (v₁, v₂, v₃), true) - elseif !_isapproxzero(c) && _iszero(f) && _iszero(e) - vals, vecs = eigen(@Tensor(x[[1,2], [1,2]]); permute=permute, scale=scale) - λ₁ = vals[1] - λ₂ = vals[2] - λ₃ = c - v₁ = vcat(vecs[:,1], 0) - v₂ = vcat(vecs[:,2], 0) - v₃ = Vec{3,T}(0,0,1) - return _eig33_construct((λ₁, λ₂, λ₃), (v₁, v₂, v₃), true) - else - # special implementation for 3x3 case (https://hal.science/hal-01501221/document) - _isapproxzero(f) && return _eig(x; permute=permute, scale=scale) - λ₁, λ₂, λ₃ = eigvals(x; permute=permute, scale=scale) - if _isapproxzero(λ₁) || _isapproxzero(λ₂) || _isapproxzero(λ₃) || - _isapproxzero(λ₁-λ₂) || _isapproxzero(λ₁-λ₃) || _isapproxzero(λ₂-λ₃) - return _eig(x; permute=permute, scale=scale) - end - v₁, v₂, v₃ = map((λ₁, λ₂, λ₃)) do λ - m = (d*(c-λ) - e*f) / (f*(b-λ) - d*e) - Vec((λ-c-e*m)/f, m, 1) - end - if !isfinite(v₁[2]) || !isfinite(v₂[2]) || !isfinite(v₃[2]) - return _eig(x; permute=permute, scale=scale) - end - return _eig33_construct((λ₁, λ₂, λ₃), (v₁, v₂, v₃), false) - end -end - -@inline function _eig33_construct((λ₁,λ₂,λ₃), (v₁,v₂,v₃), permute) - function _sortperm3(v) - local perm = SVector(1,2,3) - # unrolled bubble-sort - (v[perm[1]] > v[perm[2]]) && (perm = SVector(perm[2], perm[1], perm[3])) - (v[perm[2]] > v[perm[3]]) && (perm = SVector(perm[1], perm[3], perm[2])) - (v[perm[1]] > v[perm[2]]) && (perm = SVector(perm[2], perm[1], perm[3])) - perm - end - values = Vec(λ₁, λ₂, λ₃) - vectors = hcat(normalize(v₁), normalize(v₂), normalize(v₃)) - if permute - perm = _sortperm3(values) - return Eigen(values[perm], vectors[:,perm]) - else - return Eigen(values, vectors) - end -end - -# fallback to StaticArrays.jl -@inline function _eig(x::AbstractSymmetricSecondOrderTensor; permute::Bool, scale::Bool) - eig = eigen(Symmetric(SArray(x)); permute=permute, scale=scale) - Eigen(Tensor(eig.values), Tensor(eig.vectors)) -end diff --git a/src/ops.jl b/src/ops.jl index ce6b0d0e..eb2f81d1 100644 --- a/src/ops.jl +++ b/src/ops.jl @@ -775,6 +775,16 @@ for pv in (:true, :false) end lu(A::AbstractMat; check = true) = lu(A, Val(true); check = check) +# eigen +@inline function eigvals(x::AbstractSymmetricSecondOrderTensor; permute::Bool=true, scale::Bool=true) + Tensor(eigvals(Symmetric(SArray(x)); permute=permute, scale = scale)) +end + +@inline function eigen(x::AbstractSymmetricSecondOrderTensor; permute::Bool=true, scale::Bool=true) + eig = eigen(Symmetric(SArray(x)); permute=permute, scale=scale) + Eigen(Tensor(eig.values), Tensor(eig.vectors)) +end + # svd struct SVD{T, TU, TS, TVt} <: Factorization{T} U::TU diff --git a/test/eigen.jl b/test/eigen.jl deleted file mode 100644 index d9b6f824..00000000 --- a/test/eigen.jl +++ /dev/null @@ -1,86 +0,0 @@ -@testset "Eigen" begin - for T in (Float32, Float64) - for dim in (2, 3) - x = rand(SymmetricSecondOrderTensor{dim, T}) - @test (@inferred eigvals(x)) ≈ eigvals(Array(x)) - @test (@inferred eigen(x)).values ≈ eigen(Array(x)).values - @test (@inferred eigen(x)).vectors ≈ (@inferred eigvecs(x)) - end - end - @testset "2x2" begin - # `eigen(Symmetric(SMatrix{2,2}([1.0 1.9932760451045367e-17; 1.9932760451045367e-17 1.0])))` fails in StaticArrays.jl - @test eigvals(SymmetricSecondOrderTensor{2}([1.0 1.9932760451045367e-17; 1.9932760451045367e-17 1.0])) == Vec(1,1) - @test eigvecs(SymmetricSecondOrderTensor{2}([1.0 1.9932760451045367e-17; 1.9932760451045367e-17 1.0])) == one(Mat{2,2}) - end - # copyied from StaticArrays.jl - @testset "3x3 degenerate cases" begin - # Rank 1 - v = randn(Vec{3,Float64}) - m = symmetric(v ⊗ v, :U) - vv = sum(abs2, v) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vecs' ⋅ vecs ≈ one(Mat{3,3,Float64}) - @test vals ≈ Vec(0.0, 0.0, vv) - @test eigvals(m) ≈ vals - - # Rank 2 - v2 = randn(Vec{3,Float64}) - v2 -= dot(v,v2)*v/(vv) - v2v2 = sum(abs2, v2) - m += symmetric(v2 ⊗ v2, :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vecs' ⋅ vecs ≈ one(Mat{3,3,Float64}) - if vv < v2v2 - @test vals ≈ Vec(0.0, vv, v2v2) - else - @test vals ≈ Vec(0.0, v2v2, vv) - end - @test eigvals(m) ≈ vals - - # Degeneracy (2 large) - m = symmetric(-99*(v ⊗ v)/vv + 100*one(Mat{3,3,Float64}), :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vecs' ⋅ vecs ≈ one(Mat{3,3,Float64}) - @test vals ≈ Vec(1.0, 100.0, 100.0) - @test eigvals(m) ≈ vals - - # Degeneracy (2 small) - m = symmetric((v ⊗ v)/vv + 1e-2*one(Mat{3,3,Float64}), :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vecs' ⋅ vecs ≈ one(Mat{3,3,Float64}) - @test vals ≈ Vec(1e-2, 1e-2, 1.01) - @test eigvals(m) ≈ vals - - # Block diagonal - m = symmetric(@Mat([1.0 0.0 0.0 - 0.0 1.0 1.0 - 0.0 1.0 1.0]), :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vals ≈ [0.0, 1.0, 2.0] - @test vecs⋅diagm(Val(0) => vals)⋅vecs' ≈ m - @test eigvals(m) ≈ vals - - m = symmetric(@Mat([1.0 0.0 1.0 - 0.0 1.0 0.0 - 1.0 0.0 1.0]), :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vals ≈ [0.0, 1.0, 2.0] - @test vecs⋅diagm(Val(0) => vals)⋅vecs' ≈ m - @test eigvals(m) ≈ vals - - m = symmetric(@Mat([1.0 1.0 0.0 - 1.0 1.0 0.0 - 0.0 0.0 1.0]), :U) - vals, vecs = eigen(m)::Eigen{<:Any,<:Any,<:Mat,<:Vec} - - @test vals ≈ [0.0, 1.0, 2.0] - @test vecs⋅diagm(Val(0) => vals)⋅vecs' ≈ m - @test eigvals(m) ≈ vals - end -end diff --git a/test/ops.jl b/test/ops.jl index d11d39ee..03d57b86 100644 --- a/test/ops.jl +++ b/test/ops.jl @@ -320,6 +320,13 @@ end @test U == Fs.U @test p == Fs.p end + @testset "eigen" begin + Random.seed!(1234) + x = rand(SymmetricSecondOrderTensor{3}) + @test (@inferred eigvals(x)) ≈ eigvals(Array(x)) + @test (@inferred eigen(x)).values ≈ eigen(Array(x)).values + @test (@inferred eigen(x)).vectors ≈ (@inferred eigvecs(x)) + end @testset "svd" begin Random.seed!(1234) At = rand(Mat{3,3}) diff --git a/test/runtests.jl b/test/runtests.jl index 9e92dbcf..b2722e50 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,7 +19,6 @@ include("einsum.jl") include("ops.jl") include("continuum_mechanics.jl") include("inv.jl") -include("eigen.jl") include("voigt.jl") include("ad.jl") include("broadcast.jl")