From 44ba6f0f33d0ca023e9fa02304e1f6da4a844dab Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 6 Dec 2021 01:11:12 -0500 Subject: [PATCH 01/20] Avoid underflow in norm() --- src/linalg.jl | 33 +++++++++++++++++++++------------ test/linalg.jl | 1 + 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index d0f072eb..d0ac3dbe 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -229,14 +229,20 @@ end return :(_init_zero(a)) end - expr = :(LinearAlgebra.norm_sqr(a[1])) + expr = :(LinearAlgebra.norm_sqr(a[1]/aₘ)) for j = 2:prod(S) - expr = :($expr + LinearAlgebra.norm_sqr(a[$j])) + expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/aₘ)) end return quote $(Expr(:meta, :inline)) - @inbounds return sqrt($expr) + zero_a = _init_zero(a) + aₘ = mapreduce(norm, max, a)::typeof(zero_a) + if iszero(aₘ) + return zero_a + else + @inbounds return aₘ * sqrt($expr) + end end end @@ -251,28 +257,32 @@ end return :(_init_zero(a)) end - expr = :(norm(a[1])^p) + expr = :(norm(a[1]/aₘ)^p) for j = 2:prod(S) - expr = :($expr + norm(a[$j])^p) + expr = :($expr + norm(a[$j]/aₘ)^p) end - expr_p1 = :(norm(a[1])) + expr_p1 = :(norm(a[1]/aₘ)) for j = 2:prod(S) - expr_p1 = :($expr_p1 + norm(a[$j])) + expr_p1 = :($expr_p1 + norm(a[$j]/aₘ)) end return quote $(Expr(:meta, :inline)) - if p == Inf - return mapreduce(norm, max, a) + zero_a = _init_zero(a) + aₘ = mapreduce(norm, max, a)::typeof(zero_a) + if iszero(aₘ) + return zero_a + elseif p == Inf + return aₘ elseif p == 1 - @inbounds return $expr_p1 + @inbounds return aₘ * $expr_p1 elseif p == 2 return norm(a) elseif p == 0 return mapreduce(_norm_p0, +, a) else - @inbounds return ($expr)^(inv(p)) + @inbounds return aₘ * ($expr)^(inv(p)) end end end @@ -466,4 +476,3 @@ end # Some shimming for special linear algebra matrix types @inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Symmetric{eltype(A),typeof(A)}(A, uplo)) @inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Hermitian{eltype(A),typeof(A)}(A, uplo)) - diff --git a/test/linalg.jl b/test/linalg.jl index 7b3923ee..a4b94407 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -208,6 +208,7 @@ StaticArrays.similar_type(::Union{RotMat2,Type{RotMat2}}) = SMatrix{2,2,Float64, end @testset "normalization" begin + @test norm(SVector(0.0,1e-180)) == 1e-180 @test norm(SVector(1.0,2.0,2.0)) ≈ 3.0 @test norm(SVector(1.0,2.0,2.0),2) ≈ 3.0 @test norm(SVector(1.0,2.0,2.0),Inf) ≈ 2.0 From 1b45376e7c41f718832620ad0b969eb4df83eba3 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 6 Dec 2021 12:12:04 -0500 Subject: [PATCH 02/20] Perform return-type assertion at last step --- src/linalg.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index d0ac3dbe..5276682e 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -237,11 +237,11 @@ end return quote $(Expr(:meta, :inline)) zero_a = _init_zero(a) - aₘ = mapreduce(norm, max, a)::typeof(zero_a) + aₘ = mapreduce(norm, max, a) if iszero(aₘ) return zero_a else - @inbounds return aₘ * sqrt($expr) + @inbounds return (aₘ * sqrt($expr))::typeof(zero_a) end end end @@ -270,19 +270,19 @@ end return quote $(Expr(:meta, :inline)) zero_a = _init_zero(a) - aₘ = mapreduce(norm, max, a)::typeof(zero_a) + aₘ = mapreduce(norm, max, a) if iszero(aₘ) return zero_a elseif p == Inf return aₘ elseif p == 1 - @inbounds return aₘ * $expr_p1 + @inbounds return (aₘ * $expr_p1)::typeof(zero_a) elseif p == 2 return norm(a) elseif p == 0 return mapreduce(_norm_p0, +, a) else - @inbounds return aₘ * ($expr)^(inv(p)) + @inbounds return (aₘ * ($expr)^(inv(p)))::typeof(zero_a) end end end From 7a8156139ded6d8d5b70913091344d6223068c97 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 19:55:39 -0500 Subject: [PATCH 03/20] Make _init_zero() take AbstractArray --- src/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index 5276682e..43818502 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -217,7 +217,7 @@ end # Norms _inner_eltype(v::AbstractArray) = isempty(v) ? eltype(v) : _inner_eltype(first(v)) _inner_eltype(x::Number) = typeof(x) -@inline _init_zero(v::StaticArray) = float(norm(zero(_inner_eltype(v)))) +@inline _init_zero(v::AbstractArray) = float(norm(zero(_inner_eltype(v)))) @inline function LinearAlgebra.norm_sqr(v::StaticArray) return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) From d553e15dd78be47e279b918f120f714395e055c8 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 19:59:25 -0500 Subject: [PATCH 04/20] Define maxabs_nested() to find maximum of absolute values of entries in nested array --- src/linalg.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index 43818502..9d8e34cb 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -219,6 +219,35 @@ _inner_eltype(v::AbstractArray) = isempty(v) ? eltype(v) : _inner_eltype(first(v _inner_eltype(x::Number) = typeof(x) @inline _init_zero(v::AbstractArray) = float(norm(zero(_inner_eltype(v)))) +@inline maxabs_nested(a::Number) = abs(a) + +function maxabs_nested(a::AbstractArray) + prod(size(a)) == 0 && (return _init_zero(a)) + + m = maxabs_nested(a[1]) + for j = 2:prod(size(a)) + m = max(m, maxabs_nested(a[j])) + end + + return m +end + +@generated function maxabs_nested(a::StaticArray) + if prod(Size(a)) == 0 + return :(StaticArrays._init_zero(a)) + end + + expr = :(maxabs_nested(a[1])) + for j = 2:prod(Size(a)) + expr = :(max($expr, maxabs_nested(a[$j]))) + end + + return quote + $(Expr(:meta, :inline)) + @inbounds return $expr + end +end + @inline function LinearAlgebra.norm_sqr(v::StaticArray) return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) end From 22739f7252fcef781f18e75d810848a853b39b0e Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 20:17:06 -0500 Subject: [PATCH 05/20] Remove _norm() and define norm() directly for type stability --- src/linalg.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 9d8e34cb..f3c51ae9 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -252,14 +252,13 @@ end return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) end -@inline norm(a::StaticArray) = _norm(Size(a), a) -@generated function _norm(::Size{S}, a::StaticArray) where {S} - if prod(S) == 0 +@generated function norm(a::StaticArray) + if prod(Size(a)) == 0 return :(_init_zero(a)) end expr = :(LinearAlgebra.norm_sqr(a[1]/aₘ)) - for j = 2:prod(S) + for j = 2:prod(Size(a)) expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/aₘ)) end @@ -280,19 +279,18 @@ function _norm_p0(x) return float(norm(iszero(x) ? zero(T) : one(T))) end -@inline norm(a::StaticArray, p::Real) = _norm(Size(a), a, p) -@generated function _norm(::Size{S}, a::StaticArray, p::Real) where {S} - if prod(S) == 0 +@generated function norm(a::StaticArray, p::Real) + if prod(Size(a)) == 0 return :(_init_zero(a)) end expr = :(norm(a[1]/aₘ)^p) - for j = 2:prod(S) + for j = 2:prod(Size(a)) expr = :($expr + norm(a[$j]/aₘ)^p) end expr_p1 = :(norm(a[1]/aₘ)) - for j = 2:prod(S) + for j = 2:prod(Size(a)) expr_p1 = :($expr_p1 + norm(a[$j]/aₘ)) end From 0f498a5b45987c05064a3f51f99b63ba6fe1ae8f Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 20:18:18 -0500 Subject: [PATCH 06/20] Use maxabs_nested() to define norm() stably without underflow --- src/linalg.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f3c51ae9..bc07a31a 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -265,11 +265,11 @@ end return quote $(Expr(:meta, :inline)) zero_a = _init_zero(a) - aₘ = mapreduce(norm, max, a) + aₘ = maxabs_nested(a) if iszero(aₘ) return zero_a else - @inbounds return (aₘ * sqrt($expr))::typeof(zero_a) + @inbounds return aₘ * sqrt($expr) end end end @@ -294,22 +294,27 @@ end expr_p1 = :($expr_p1 + norm(a[$j]/aₘ)) end + expr_pInf = :(norm(a[1]/aₘ)) + for j = 2:prod(Size(a)) + expr_pInf = :(max($expr_pInf, norm(a[$j]/aₘ))) + end + return quote $(Expr(:meta, :inline)) zero_a = _init_zero(a) - aₘ = mapreduce(norm, max, a) + aₘ = maxabs_nested(a) if iszero(aₘ) return zero_a elseif p == Inf - return aₘ + return aₘ * $expr_pInf elseif p == 1 - @inbounds return (aₘ * $expr_p1)::typeof(zero_a) + @inbounds return aₘ * $expr_p1 elseif p == 2 return norm(a) elseif p == 0 return mapreduce(_norm_p0, +, a) else - @inbounds return (aₘ * ($expr)^(inv(p)))::typeof(zero_a) + @inbounds return aₘ * ($expr)^(inv(p)) end end end From ff314a1fa842e5c6a1dde54d174df63b9195acde Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 20:18:54 -0500 Subject: [PATCH 07/20] Call _init_zero() without qualification by StaticArrays --- src/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index bc07a31a..17a1abf4 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -234,7 +234,7 @@ end @generated function maxabs_nested(a::StaticArray) if prod(Size(a)) == 0 - return :(StaticArrays._init_zero(a)) + return :(_init_zero(a)) end expr = :(maxabs_nested(a[1])) From e75bad8173b64d9ef6fda320cd44a08e986605a5 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 22:28:22 -0500 Subject: [PATCH 08/20] Define norm_scaled() and call it inside norm() --- src/linalg.jl | 81 +++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 62 insertions(+), 19 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 17a1abf4..c4587f05 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -252,11 +252,7 @@ end return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) end -@generated function norm(a::StaticArray) - if prod(Size(a)) == 0 - return :(_init_zero(a)) - end - +@generated function norm_scaled(a::StaticArray) expr = :(LinearAlgebra.norm_sqr(a[1]/aₘ)) for j = 2:prod(Size(a)) expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/aₘ)) @@ -274,16 +270,33 @@ end end end +@generated function norm(a::StaticArray) + if prod(Size(a)) == 0 + return :(_init_zero(a)) + end + + expr = :(LinearAlgebra.norm_sqr(a[1])) + for j = 2:prod(Size(a)) + expr = :($expr + LinearAlgebra.norm_sqr(a[$j])) + end + + return quote + $(Expr(:meta, :inline)) + l = @inbounds sqrt($expr) + if iszero(l) || isinf(l) + return norm_scaled(a) + else + return l + end + end +end + function _norm_p0(x) T = _inner_eltype(x) return float(norm(iszero(x) ? zero(T) : one(T))) end -@generated function norm(a::StaticArray, p::Real) - if prod(Size(a)) == 0 - return :(_init_zero(a)) - end - +@generated function norm_scaled(a::StaticArray, p::Real) expr = :(norm(a[1]/aₘ)^p) for j = 2:prod(Size(a)) expr = :($expr + norm(a[$j]/aₘ)^p) @@ -294,31 +307,61 @@ end expr_p1 = :($expr_p1 + norm(a[$j]/aₘ)) end - expr_pInf = :(norm(a[1]/aₘ)) - for j = 2:prod(Size(a)) - expr_pInf = :(max($expr_pInf, norm(a[$j]/aₘ))) - end - return quote $(Expr(:meta, :inline)) zero_a = _init_zero(a) aₘ = maxabs_nested(a) if iszero(aₘ) return zero_a - elseif p == Inf - return aₘ * $expr_pInf + elseif p == 0 + return mapreduce(_norm_p0, +, a) elseif p == 1 @inbounds return aₘ * $expr_p1 elseif p == 2 return norm(a) - elseif p == 0 - return mapreduce(_norm_p0, +, a) else @inbounds return aₘ * ($expr)^(inv(p)) end end end +@generated function norm(a::StaticArray, p::Real) + if prod(Size(a)) == 0 + return :(_init_zero(a)) + end + + expr = :(norm(a[1])^p) + for j = 2:prod(Size(a)) + expr = :($expr + norm(a[$j])^p) + end + + expr_p1 = :(norm(a[1])) + for j = 2:prod(Size(a)) + expr_p1 = :($expr_p1 + norm(a[$j])) + end + + return quote + $(Expr(:meta, :inline)) + p == Inf && return mapreduce(norm, max, a) # no need to scale + p == 2 && return norm(a) # norm(a) takes care of scaling + + local l + if p == 0 + l = mapreduce(_norm_p0, +, a) + elseif p == 1 + l = @inbounds $expr_p1 + else + l = @inbounds ($expr)^(inv(p)) + end + + if iszero(l) || isinf(l) + return norm_scaled(a, p) + else + return l + end + end +end + @inline normalize(a::StaticArray) = inv(norm(a))*a @inline normalize(a::StaticArray, p::Real) = inv(norm(a, p))*a From f33208aae4d3b67028e6b61cdd8a2ddb9e2ec88d Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 22:52:33 -0500 Subject: [PATCH 09/20] norm(a,p) does not call norm_scaled(a,p) if p == 0, 2, Inf --- src/linalg.jl | 37 ++++++++++--------------------------- 1 file changed, 10 insertions(+), 27 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index c4587f05..10c2dd52 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -309,19 +309,12 @@ end return quote $(Expr(:meta, :inline)) - zero_a = _init_zero(a) aₘ = maxabs_nested(a) - if iszero(aₘ) - return zero_a - elseif p == 0 - return mapreduce(_norm_p0, +, a) - elseif p == 1 - @inbounds return aₘ * $expr_p1 - elseif p == 2 - return norm(a) - else - @inbounds return aₘ * ($expr)^(inv(p)) - end + + iszero(aₘ) && return aₘ + p == 1 && return @inbounds aₘ * $expr_p1 + p == 2 && return norm(a) + return @inbounds aₘ * ($expr)^(inv(p)) end end @@ -342,23 +335,13 @@ end return quote $(Expr(:meta, :inline)) - p == Inf && return mapreduce(norm, max, a) # no need to scale + p == 0 && return mapreduce(_norm_p0, +, a) # no need for scaling p == 2 && return norm(a) # norm(a) takes care of scaling + p == Inf && return mapreduce(norm, max, a) # no need for scaling - local l - if p == 0 - l = mapreduce(_norm_p0, +, a) - elseif p == 1 - l = @inbounds $expr_p1 - else - l = @inbounds ($expr)^(inv(p)) - end - - if iszero(l) || isinf(l) - return norm_scaled(a, p) - else - return l - end + l = p==1 ? @inbounds($expr_p1) : @inbounds(($expr)^(inv(p))) + (iszero(l) || isinf(l)) && return norm_scaled(a, p) + return l end end From 629d020987e01b83e6f018c55fd3b355e58a29e3 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 22:53:03 -0500 Subject: [PATCH 10/20] Make code more succinct --- src/linalg.jl | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 10c2dd52..d405a35d 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -260,13 +260,10 @@ end return quote $(Expr(:meta, :inline)) - zero_a = _init_zero(a) aₘ = maxabs_nested(a) - if iszero(aₘ) - return zero_a - else - @inbounds return aₘ * sqrt($expr) - end + + iszero(aₘ) && return aₘ + return @inbounds aₘ * sqrt($expr) end end @@ -283,11 +280,9 @@ end return quote $(Expr(:meta, :inline)) l = @inbounds sqrt($expr) - if iszero(l) || isinf(l) - return norm_scaled(a) - else - return l - end + + (iszero(l) || isinf(l)) && return norm_scaled(a) + return l end end From 2e3c6b6c7b0bda616351efbdda91e65f4d4ec7d2 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 22:54:02 -0500 Subject: [PATCH 11/20] Use normInf() instead of maxabs_nested() --- src/linalg.jl | 33 ++------------------------------- 1 file changed, 2 insertions(+), 31 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index d405a35d..379b203c 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -219,35 +219,6 @@ _inner_eltype(v::AbstractArray) = isempty(v) ? eltype(v) : _inner_eltype(first(v _inner_eltype(x::Number) = typeof(x) @inline _init_zero(v::AbstractArray) = float(norm(zero(_inner_eltype(v)))) -@inline maxabs_nested(a::Number) = abs(a) - -function maxabs_nested(a::AbstractArray) - prod(size(a)) == 0 && (return _init_zero(a)) - - m = maxabs_nested(a[1]) - for j = 2:prod(size(a)) - m = max(m, maxabs_nested(a[j])) - end - - return m -end - -@generated function maxabs_nested(a::StaticArray) - if prod(Size(a)) == 0 - return :(_init_zero(a)) - end - - expr = :(maxabs_nested(a[1])) - for j = 2:prod(Size(a)) - expr = :(max($expr, maxabs_nested(a[$j]))) - end - - return quote - $(Expr(:meta, :inline)) - @inbounds return $expr - end -end - @inline function LinearAlgebra.norm_sqr(v::StaticArray) return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) end @@ -260,7 +231,7 @@ end return quote $(Expr(:meta, :inline)) - aₘ = maxabs_nested(a) + aₘ = LinearAlgebra.normInf(a) iszero(aₘ) && return aₘ return @inbounds aₘ * sqrt($expr) @@ -304,7 +275,7 @@ end return quote $(Expr(:meta, :inline)) - aₘ = maxabs_nested(a) + aₘ = LinearAlgebra.normInf(a) iszero(aₘ) && return aₘ p == 1 && return @inbounds aₘ * $expr_p1 From 146d0e9b79efacd038d32077fa49ecc168d55a79 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sun, 2 Jan 2022 22:55:03 -0500 Subject: [PATCH 12/20] =?UTF-8?q?Rename=20a=E2=82=98=20to=20scale?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/linalg.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 379b203c..f93f0bd4 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -224,17 +224,17 @@ _inner_eltype(x::Number) = typeof(x) end @generated function norm_scaled(a::StaticArray) - expr = :(LinearAlgebra.norm_sqr(a[1]/aₘ)) + expr = :(LinearAlgebra.norm_sqr(a[1]/scale)) for j = 2:prod(Size(a)) - expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/aₘ)) + expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/scale)) end return quote $(Expr(:meta, :inline)) - aₘ = LinearAlgebra.normInf(a) + scale = LinearAlgebra.normInf(a) - iszero(aₘ) && return aₘ - return @inbounds aₘ * sqrt($expr) + iszero(scale) && return scale + return @inbounds scale * sqrt($expr) end end @@ -263,24 +263,24 @@ function _norm_p0(x) end @generated function norm_scaled(a::StaticArray, p::Real) - expr = :(norm(a[1]/aₘ)^p) + expr = :(norm(a[1]/scale)^p) for j = 2:prod(Size(a)) - expr = :($expr + norm(a[$j]/aₘ)^p) + expr = :($expr + norm(a[$j]/scale)^p) end - expr_p1 = :(norm(a[1]/aₘ)) + expr_p1 = :(norm(a[1]/scale)) for j = 2:prod(Size(a)) - expr_p1 = :($expr_p1 + norm(a[$j]/aₘ)) + expr_p1 = :($expr_p1 + norm(a[$j]/scale)) end return quote $(Expr(:meta, :inline)) - aₘ = LinearAlgebra.normInf(a) + scale = LinearAlgebra.normInf(a) - iszero(aₘ) && return aₘ - p == 1 && return @inbounds aₘ * $expr_p1 + iszero(scale) && return scale + p == 1 && return @inbounds scale * $expr_p1 p == 2 && return norm(a) - return @inbounds aₘ * ($expr)^(inv(p)) + return @inbounds scale * ($expr)^(inv(p)) end end From 90e7c70dbbfd7bebf68d104466c39c9cbc8a59a8 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Mon, 3 Jan 2022 15:57:30 -0500 Subject: [PATCH 13/20] Calculate scale with mapreduce() --- src/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f93f0bd4..f4c92c62 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -231,7 +231,7 @@ end return quote $(Expr(:meta, :inline)) - scale = LinearAlgebra.normInf(a) + scale = mapreduce(norm, @fastmath(max), a) iszero(scale) && return scale return @inbounds scale * sqrt($expr) @@ -275,7 +275,7 @@ end return quote $(Expr(:meta, :inline)) - scale = LinearAlgebra.normInf(a) + scale = mapreduce(norm, @fastmath(max), a) iszero(scale) && return scale p == 1 && return @inbounds scale * $expr_p1 From e44647c6076f7220ff7616b9545ab70135c326af Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 22:43:22 -0500 Subject: [PATCH 14/20] Re-introduce maxabs_nested to avoid type instability --- src/linalg.jl | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f4c92c62..f103d844 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -223,6 +223,18 @@ _inner_eltype(x::Number) = typeof(x) return mapreduce(LinearAlgebra.norm_sqr, +, v; init=_init_zero(v)) end +@inline maxabs_nested(a::Number) = abs(a) +function maxabs_nested(a::AbstractArray) + prod(size(a)) == 0 && (return _init_zero(a)) + + m = maxabs_nested(a[1]) + for j = 2:prod(size(a)) + m = @fastmath max(m, maxabs_nested(a[j])) + end + + return m +end + @generated function norm_scaled(a::StaticArray) expr = :(LinearAlgebra.norm_sqr(a[1]/scale)) for j = 2:prod(Size(a)) @@ -231,9 +243,9 @@ end return quote $(Expr(:meta, :inline)) - scale = mapreduce(norm, @fastmath(max), a) + scale = maxabs_nested(a) - iszero(scale) && return scale + iszero(scale) && return _init_zero(a) return @inbounds scale * sqrt($expr) end end @@ -275,9 +287,9 @@ end return quote $(Expr(:meta, :inline)) - scale = mapreduce(norm, @fastmath(max), a) + scale = maxabs_nested(a) - iszero(scale) && return scale + iszero(scale) && return _init_zero(a) p == 1 && return @inbounds scale * $expr_p1 p == 2 && return norm(a) return @inbounds scale * ($expr)^(inv(p)) From 02172031b97cbff3cdb373810e5a8ffd8e68f636 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 23:23:07 -0500 Subject: [PATCH 15/20] norm() calls _norm() and _norm_scaled() --- src/linalg.jl | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index f103d844..c074aec9 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -235,9 +235,9 @@ function maxabs_nested(a::AbstractArray) return m end -@generated function norm_scaled(a::StaticArray) +@generated function _norm_scaled(::Size{S}, a::StaticArray) where {S} expr = :(LinearAlgebra.norm_sqr(a[1]/scale)) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr = :($expr + LinearAlgebra.norm_sqr(a[$j]/scale)) end @@ -250,13 +250,12 @@ end end end -@generated function norm(a::StaticArray) - if prod(Size(a)) == 0 - return :(_init_zero(a)) - end +@inline norm(a::StaticArray) = _norm(Size(a), a) +@generated function _norm(::Size{S}, a::StaticArray) where {S} + prod(S) == 0 && return :(_init_zero(a)) expr = :(LinearAlgebra.norm_sqr(a[1])) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr = :($expr + LinearAlgebra.norm_sqr(a[$j])) end @@ -264,7 +263,7 @@ end $(Expr(:meta, :inline)) l = @inbounds sqrt($expr) - (iszero(l) || isinf(l)) && return norm_scaled(a) + (iszero(l) || isinf(l)) && return _norm_scaled(Size(a), a) return l end end @@ -274,14 +273,15 @@ function _norm_p0(x) return float(norm(iszero(x) ? zero(T) : one(T))) end -@generated function norm_scaled(a::StaticArray, p::Real) +# Do not need to deal with p == 0, 2, Inf; see norm(a, p). +@generated function _norm_scaled(::Size{S}, a::StaticArray, p::Real) where {S} expr = :(norm(a[1]/scale)^p) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr = :($expr + norm(a[$j]/scale)^p) end expr_p1 = :(norm(a[1]/scale)) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr_p1 = :($expr_p1 + norm(a[$j]/scale)) end @@ -296,18 +296,17 @@ end end end -@generated function norm(a::StaticArray, p::Real) - if prod(Size(a)) == 0 - return :(_init_zero(a)) - end +@inline norm(a::StaticArray, p::Real) = _norm(Size(a), a, p) +@generated function _norm(::Size{S}, a::StaticArray, p::Real) where {S} + prod(S) == 0 && return :(_init_zero(a)) expr = :(norm(a[1])^p) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr = :($expr + norm(a[$j])^p) end expr_p1 = :(norm(a[1])) - for j = 2:prod(Size(a)) + for j = 2:prod(S) expr_p1 = :($expr_p1 + norm(a[$j])) end @@ -318,7 +317,7 @@ end p == Inf && return mapreduce(norm, max, a) # no need for scaling l = p==1 ? @inbounds($expr_p1) : @inbounds(($expr)^(inv(p))) - (iszero(l) || isinf(l)) && return norm_scaled(a, p) + (iszero(l) || isinf(l)) && return _norm_scaled(Size(a), a, p) # p != 0, 2, Inf return l end end From 2a66d8ef1e55e8b608b1eed34d6a0d7d197b4aa2 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 23:24:00 -0500 Subject: [PATCH 16/20] _norm_scaled(Size(a), a, p) does not need to deal with p==2 --- src/linalg.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index c074aec9..8e7ec0f6 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -291,7 +291,6 @@ end iszero(scale) && return _init_zero(a) p == 1 && return @inbounds scale * $expr_p1 - p == 2 && return norm(a) return @inbounds scale * ($expr)^(inv(p)) end end From 51cdb73a3f61aedd63ebaa5ee9d57afc0fcdfeae Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 23:24:49 -0500 Subject: [PATCH 17/20] Add test for overflow --- test/linalg.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/linalg.jl b/test/linalg.jl index a4b94407..689bab11 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -208,7 +208,10 @@ StaticArrays.similar_type(::Union{RotMat2,Type{RotMat2}}) = SMatrix{2,2,Float64, end @testset "normalization" begin - @test norm(SVector(0.0,1e-180)) == 1e-180 + @test norm(SVector(0.0,1e-180)) == 1e-180 # avoid underflow + @test all([norm(SVector(0.0,1e-180), p) == 1e-180 for p = [2,3,Inf]]) # avoid underflow + @test norm(SVector(0.0,1e155)) == 1e155 # avoid overflow + @test all([norm(SVector(0.0,1e155), p) == 1e155 for p = [2,3,Inf]]) # avoid overflow @test norm(SVector(1.0,2.0,2.0)) ≈ 3.0 @test norm(SVector(1.0,2.0,2.0),2) ≈ 3.0 @test norm(SVector(1.0,2.0,2.0),Inf) ≈ 2.0 From 921da5980dffbc5f562a78c0af600a54aca70a27 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 23:25:48 -0500 Subject: [PATCH 18/20] =?UTF-8?q?Add=20tests=20for=20underflow=20and=20ove?= =?UTF-8?q?rflow=20for=20p=20=E2=89=A5=202?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/linalg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg.jl b/test/linalg.jl index 689bab11..ae459f60 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -209,8 +209,8 @@ StaticArrays.similar_type(::Union{RotMat2,Type{RotMat2}}) = SMatrix{2,2,Float64, @testset "normalization" begin @test norm(SVector(0.0,1e-180)) == 1e-180 # avoid underflow - @test all([norm(SVector(0.0,1e-180), p) == 1e-180 for p = [2,3,Inf]]) # avoid underflow @test norm(SVector(0.0,1e155)) == 1e155 # avoid overflow + @test all([norm(SVector(0.0,1e-180), p) == 1e-180 for p = [2,3,Inf]]) # avoid underflow @test all([norm(SVector(0.0,1e155), p) == 1e155 for p = [2,3,Inf]]) # avoid overflow @test norm(SVector(1.0,2.0,2.0)) ≈ 3.0 @test norm(SVector(1.0,2.0,2.0),2) ≈ 3.0 From 9ea186b6e1b8c9de327c50f03ac96db6c6581866 Mon Sep 17 00:00:00 2001 From: Wonseok Shin Date: Sat, 8 Jan 2022 23:50:09 -0500 Subject: [PATCH 19/20] Avoid using iszero(l) and isinf(l) for speedup --- src/linalg.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 8e7ec0f6..6bd435cd 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -263,8 +263,8 @@ end $(Expr(:meta, :inline)) l = @inbounds sqrt($expr) - (iszero(l) || isinf(l)) && return _norm_scaled(Size(a), a) - return l + 0 Date: Sat, 8 Jan 2022 23:53:35 -0500 Subject: [PATCH 20/20] Avoid using iszero(scale) for speedup --- src/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 6bd435cd..c5c0bcf3 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -245,7 +245,7 @@ end $(Expr(:meta, :inline)) scale = maxabs_nested(a) - iszero(scale) && return _init_zero(a) + scale==0 && return _init_zero(a) return @inbounds scale * sqrt($expr) end end @@ -289,7 +289,7 @@ end $(Expr(:meta, :inline)) scale = maxabs_nested(a) - iszero(scale) && return _init_zero(a) + scale==0 && return _init_zero(a) p == 1 && return @inbounds scale * $expr_p1 return @inbounds scale * ($expr)^(inv(p)) end