From 3e696ed95b5f3c687d7dd344730bf3f6fb44d4ca Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 29 Nov 2021 12:22:19 -0500 Subject: [PATCH 01/11] faster and simpler generic_norm2 I'm not 100% sure if this works. Probably needs a PkgEval --- stdlib/LinearAlgebra/src/generic.jl | 21 +++------------------ 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 5db2b525ee584..5f6eeaa06a6f7 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -462,27 +462,12 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x)) function generic_norm2(x) maxabs = normInf(x) - (maxabs == 0 || isinf(maxabs)) && return maxabs - (v, s) = iterate(x)::Tuple + (ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs T = typeof(maxabs) if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary - sum::promote_type(Float64, T) = norm_sqr(v) - while true - y = iterate(x, s) - y === nothing && break - (v, s) = y - sum += norm_sqr(v) - end - return convert(T, sqrt(sum)) + return convert(T, sqrt(mapreduce(norm_sqr, +, x))) else - sum = abs2(norm(v)/maxabs) - while true - y = iterate(x, s) - y === nothing && break - (v, s) = y - sum += (norm(v)/maxabs)^2 - end - return convert(T, maxabs*sqrt(sum)) + return convert(T, maxabs*sqrt(mapreduce(v -> abs2(norm(v)/maxabs), +, x))) end end From 604adc36f1674046600dd62c212dd58f1267754c Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 2 Dec 2021 16:31:01 -0500 Subject: [PATCH 02/11] accumulate norm2 in at least Float64 precision and convert generic_normp --- stdlib/LinearAlgebra/src/generic.jl | 41 ++++++++++------------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 5f6eeaa06a6f7..4027893ef215d 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -449,7 +449,7 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons # Dot products and norms # special cases of norm; note that they don't need to handle isempty(x) -generic_normMinusInf(x) = float(mapreduce(norm, min, x)) +2MinusInf(x) = float(mapreduce(norm, min, x)) generic_normInf(x) = float(mapreduce(norm, max, x)) @@ -464,44 +464,31 @@ function generic_norm2(x) maxabs = normInf(x) (ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs T = typeof(maxabs) - if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary - return convert(T, sqrt(mapreduce(norm_sqr, +, x))) + sT = promote_type(Float64, T) + maxabs_st = sT(maxabs) + if isfinite(length(x)*maxabs_st*maxabs_st) && maxabs_st*maxabs_st != 0 # Scaling not necessary + return convert(T, sqrt(mapreduce(v -> sT(norm_sqr(v)), +, x))) else - return convert(T, maxabs*sqrt(mapreduce(v -> abs2(norm(v)/maxabs), +, x))) + invmaxabs = inv(maxabs_st) + return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))*invmaxabs), +, x))) end end # Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p) # (Not technically a "norm" for p < 1.) function generic_normp(x, p) - (v, s) = iterate(x)::Tuple + T = typeof(float(first(x))) + spp::promote_type(Float64, T) = p if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) (maxabs == 0 || isinf(maxabs)) && return maxabs - T = typeof(maxabs) - else - T = typeof(float(norm(v))) - end - spp::promote_type(Float64, T) = p - if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && maxabs^spp != 0) # scaling not necessary - sum::promote_type(Float64, T) = norm(v)^spp - while true - y = iterate(x, s) - y === nothing && break - (v, s) = y - sum += norm(v)^spp - end - return convert(T, sum^inv(spp)) - else # rescaling - sum = (norm(v)/maxabs)^spp - while true - y = iterate(x, s) - y === nothing && break - (v, s) = y - sum += (norm(v)/maxabs)^spp + maxelnorm = maxabs^spp + if (isfinite(length(x)*maxelnorm) && maxelnorm != 0) # rescaling necessary + invmaxabs = inv(maxabs_st) + return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) end - return convert(T, maxabs*sum^inv(spp)) end + return convert(T, mapreduce(v -> norm(v)^spp, +, x))^inv(spp)) end normMinusInf(x) = generic_normMinusInf(x) From 44509daf2302ac0f92c1e6ae7073a7a515739031 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 3 Dec 2021 08:03:24 -0500 Subject: [PATCH 03/11] Update stdlib/LinearAlgebra/src/generic.jl Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 4027893ef215d..0072b855ad570 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -477,7 +477,7 @@ end # Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p) # (Not technically a "norm" for p < 1.) function generic_normp(x, p) - T = typeof(float(first(x))) + T = typeof(float(norm(first(x)))) spp::promote_type(Float64, T) = p if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) From aedb5199fbe269199b2127c039a184513e5bad5a Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 3 Dec 2021 09:46:26 -0500 Subject: [PATCH 04/11] fix typo --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 0072b855ad570..874fa74503636 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -449,7 +449,7 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons # Dot products and norms # special cases of norm; note that they don't need to handle isempty(x) -2MinusInf(x) = float(mapreduce(norm, min, x)) +generic_normMinusInf(x) = float(mapreduce(norm, min, x)) generic_normInf(x) = float(mapreduce(norm, max, x)) From b8c8b85673b277f776e1cd12d280bc975cec9e45 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 3 Dec 2021 10:30:54 -0500 Subject: [PATCH 05/11] fix typo --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 874fa74503636..363b39af03722 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -488,7 +488,7 @@ function generic_normp(x, p) return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) end end - return convert(T, mapreduce(v -> norm(v)^spp, +, x))^inv(spp)) + return convert(T, mapreduce(v -> norm(v)^spp, +, x)^inv(spp)) end normMinusInf(x) = generic_normMinusInf(x) From 689032c7f4813c3f50d27c7980a4e66c16c0be07 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 3 Dec 2021 11:48:15 -0500 Subject: [PATCH 06/11] Update stdlib/LinearAlgebra/src/generic.jl Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 363b39af03722..a30d673606862 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -484,7 +484,7 @@ function generic_normp(x, p) (maxabs == 0 || isinf(maxabs)) && return maxabs maxelnorm = maxabs^spp if (isfinite(length(x)*maxelnorm) && maxelnorm != 0) # rescaling necessary - invmaxabs = inv(maxabs_st) + invmaxabs = inv(maxabs) return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) end end From 812bd577b37585b5eaa380a9ea86adedc0e25a22 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 3 Dec 2021 14:13:16 -0500 Subject: [PATCH 07/11] fix norm of subnormals --- stdlib/LinearAlgebra/src/generic.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index a30d673606862..e5126af4f6b3b 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -470,7 +470,10 @@ function generic_norm2(x) return convert(T, sqrt(mapreduce(v -> sT(norm_sqr(v)), +, x))) else invmaxabs = inv(maxabs_st) - return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))*invmaxabs), +, x))) + if isfinite(invmaxabs) + return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))*invmaxabs), +, x))) + end + return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))/maxabs_st), +, x))) end end @@ -482,10 +485,13 @@ function generic_normp(x, p) if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) (maxabs == 0 || isinf(maxabs)) && return maxabs - maxelnorm = maxabs^spp - if (isfinite(length(x)*maxelnorm) && maxelnorm != 0) # rescaling necessary + max_el_norm = maxabs^spp + if (isfinite(length(x)*max_el_norm) && max_el_norm != 0) # rescaling necessary invmaxabs = inv(maxabs) - return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) + if isfinite(invmaxabs) + return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) + end + return convert(T, maxabs*mapreduce(v -> (norm(v)/maxabs)^spp, +, x)^inv(spp)) end end return convert(T, mapreduce(v -> norm(v)^spp, +, x)^inv(spp)) From 43066430eeeec46d3b31fd55cbfcb81cc6a4b329 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 3 Dec 2021 15:30:30 -0500 Subject: [PATCH 08/11] test was backwards --- stdlib/LinearAlgebra/src/generic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index e5126af4f6b3b..cd0215d869523 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -486,7 +486,7 @@ function generic_normp(x, p) maxabs = p > 1 ? normInf(x) : normMinusInf(x) (maxabs == 0 || isinf(maxabs)) && return maxabs max_el_norm = maxabs^spp - if (isfinite(length(x)*max_el_norm) && max_el_norm != 0) # rescaling necessary + if !(isfinite(length(x)*max_el_norm) && max_el_norm != 0) # rescaling necessary invmaxabs = inv(maxabs) if isfinite(invmaxabs) return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) From 3de392abdf87b910f94c81633d4b47b1e78e14f4 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 3 Dec 2021 17:17:55 -0500 Subject: [PATCH 09/11] fix type stability --- stdlib/LinearAlgebra/src/generic.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index cd0215d869523..fd2b79c0ee109 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -467,13 +467,13 @@ function generic_norm2(x) sT = promote_type(Float64, T) maxabs_st = sT(maxabs) if isfinite(length(x)*maxabs_st*maxabs_st) && maxabs_st*maxabs_st != 0 # Scaling not necessary - return convert(T, sqrt(mapreduce(v -> sT(norm_sqr(v)), +, x))) + return convert(T, sqrt(mapreduce(v -> convert(sT, norm_sqr(v)), +, x))) else invmaxabs = inv(maxabs_st) if isfinite(invmaxabs) - return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))*invmaxabs), +, x))) + return convert(T, maxabs*sqrt(mapreduce(v -> abs2((norm(v))*invmaxabs), +, x))) end - return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))/maxabs_st), +, x))) + return convert(T, maxabs*sqrt(mapreduce(v -> abs2((norm(v))/maxabs_st), +, x))) end end From d9de30b8d1f7d0c359af83bf270cfb147caa899b Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 28 Dec 2021 13:38:13 -0500 Subject: [PATCH 10/11] be more lazy --- stdlib/LinearAlgebra/src/generic.jl | 51 +++++++++++++---------------- 1 file changed, 23 insertions(+), 28 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index fd2b79c0ee109..b23e860dfe8f9 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -457,44 +457,39 @@ generic_norm1(x) = mapreduce(float ∘ norm, +, x) # faster computation of norm(x)^2, avoiding overflow for integers norm_sqr(x) = norm(x)^2 -norm_sqr(x::Number) = abs2(x) -norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x)) +norm_sqr(x::T) where {T<:Number} = abs2(promote_type(T,Float64)(x)) function generic_norm2(x) - maxabs = normInf(x) - (ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs - T = typeof(maxabs) - sT = promote_type(Float64, T) - maxabs_st = sT(maxabs) - if isfinite(length(x)*maxabs_st*maxabs_st) && maxabs_st*maxabs_st != 0 # Scaling not necessary - return convert(T, sqrt(mapreduce(v -> convert(sT, norm_sqr(v)), +, x))) - else - invmaxabs = inv(maxabs_st) - if isfinite(invmaxabs) - return convert(T, maxabs*sqrt(mapreduce(v -> abs2((norm(v))*invmaxabs), +, x))) - end - return convert(T, maxabs*sqrt(mapreduce(v -> abs2((norm(v))/maxabs_st), +, x))) + isempty(x) && return float(norm(zero(eltype(itr)))) + T = typeof(float(norm(first(x)))) + sT = promote_type(T, Float64) + ans = mapreduce(norm_sqr, +, x) + ans in (0, Inf) || return convert(T, sqrt(ans)) + maxabs = sT(normInf(x)) + ans = sT(0) + for v in x + ans += (norm(v)/maxabs)^2 end + return convert(T, maxabs*sqrt(ans)) end # Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p) # (Not technically a "norm" for p < 1.) function generic_normp(x, p) + isempty(x) && return float(norm(zero(eltype(itr)))) T = typeof(float(norm(first(x)))) - spp::promote_type(Float64, T) = p - if p > 1 || p < -1 # might need to rescale to avoid overflow - maxabs = p > 1 ? normInf(x) : normMinusInf(x) - (maxabs == 0 || isinf(maxabs)) && return maxabs - max_el_norm = maxabs^spp - if !(isfinite(length(x)*max_el_norm) && max_el_norm != 0) # rescaling necessary - invmaxabs = inv(maxabs) - if isfinite(invmaxabs) - return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp)) - end - return convert(T, maxabs*mapreduce(v -> (norm(v)/maxabs)^spp, +, x)^inv(spp)) - end + sT = promote_type(T, Float64) + ans = sT(0) + for v in x + ans += sT(norm(v))^p + end + ans in (0, Inf) || return convert(T, ans^inv(spp)) + maxabs = p > 1 ? normInf(x) : normMinusInf(x) + ans = sT(0) + for v in x + ans += (sT(norm(v))/maxabs)^p end - return convert(T, mapreduce(v -> norm(v)^spp, +, x)^inv(spp)) + return convert(T, ans^inv(p)) end normMinusInf(x) = generic_normMinusInf(x) From 543f5632f23740c9ab9a46ada7d9ba9220de4524 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 28 Dec 2021 14:59:18 -0500 Subject: [PATCH 11/11] Update stdlib/LinearAlgebra/src/generic.jl Co-authored-by: Michael Abbott <32575566+mcabbott@users.noreply.github.com> --- stdlib/LinearAlgebra/src/generic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index b23e860dfe8f9..9925f565550d0 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -460,8 +460,8 @@ norm_sqr(x) = norm(x)^2 norm_sqr(x::T) where {T<:Number} = abs2(promote_type(T,Float64)(x)) function generic_norm2(x) - isempty(x) && return float(norm(zero(eltype(itr)))) - T = typeof(float(norm(first(x)))) + isempty(x) && return norm(zero(eltype(x))) + T = typeof(norm(first(x))) sT = promote_type(T, Float64) ans = mapreduce(norm_sqr, +, x) ans in (0, Inf) || return convert(T, sqrt(ans))