From 75a340894a30fed8a5ce7ce67cb3343bb794a080 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 20 May 2021 06:28:43 +0200 Subject: [PATCH] Use numerically robust norm This PR is mostly a FYI regarding https://github.com/JuliaArrays/StaticArrays.jl/issues/913 Depending on whether or not that issue is closed, you may want to switch to explicitly calling `generic_norm2` on static vectors to circumvent the accuracy issue. The issue arises when combining static arrays with ForwardDiff, in our case it occured in `exp(::SkewSymmetric)`. I did some benchmarking, and `generic_norm2` is faster than `norm` for all standard arrays up to at least length 9. For static arrays, it appears to do okay as well, while avoiding the accuracy issue. ```julia julia> a 3-element SVector{3, Int64} with indices SOneTo(3): 1 2 3 julia> @btime norm($(Ref(a))[]) # standard norm of static array 5.135 ns (0 allocations: 0 bytes) 3.74166 julia> @btime norm($(Ref(Vector(a)))[]) # standard norm of array 8.545 ns (0 allocations: 0 bytes) 3.74166 julia> @btime LinearAlgebra.generic_norm2($(Ref((a)))[]) # generic_norm2 of static array 4.490 ns (0 allocations: 0 bytes) 3.74166 julia> @btime LinearAlgebra.generic_norm2($(Ref(Vector(a)))[]) # generic_norm2 of array 7.671 ns (0 allocations: 0 bytes) 3.74166 ``` --- src/spatial/threevectors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spatial/threevectors.jl b/src/spatial/threevectors.jl index e6894f7d..b2b29537 100644 --- a/src/spatial/threevectors.jl +++ b/src/spatial/threevectors.jl @@ -75,7 +75,7 @@ LinearAlgebra.cross(v1::FreeVector3D, v2::FreeVector3D) = begin @framecheck(v1.f LinearAlgebra.dot(v1::FreeVector3D, v2::FreeVector3D) = begin @framecheck(v1.frame, v2.frame); v1.v ⋅ v2.v end Base.:*(t::Transform3D, vector::FreeVector3D) = begin @framecheck(t.from, vector.frame); FreeVector3D(t.to, rotation(t) * vector.v) end Base.:\(t::Transform3D, point::FreeVector3D) = begin @framecheck point.frame t.to; FreeVector3D(t.from, rotation(t) \ point.v) end -LinearAlgebra.norm(v::FreeVector3D) = norm(v.v) +LinearAlgebra.norm(v::FreeVector3D) = LinearAlgebra.generic_norm2(v.v) LinearAlgebra.normalize(v::FreeVector3D, p = 2) = FreeVector3D(v.frame, normalize(v.v, p)) # Mixed Point3D and FreeVector3D