Specialize norm-related functions for Diagonal

Author: jishnub

src/diagonal.jl

@@ -1246,3 +1246,14 @@ function fillband!(D::Diagonal, x, l, u)
     end
     return D
 end
+
+# norm
+function generic_normMinusInf(D::Diagonal)
+    norm_diag = generic_normMinusInf(D.diag)
+    min(norm_diag, zero(norm_diag))
+end
+generic_normInf(D::Diagonal) = generic_normInf(D.diag)
+generic_norm1(D::Diagonal) = generic_norm1(D.diag)
+_generic_norm2(D::Diagonal, maxabs) = _generic_norm2(D.diag, maxabs)
+_generic_normp(D::Diagonal, p, maxabs) = _generic_normp(D.diag, p, maxabs)
+norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag)

src/generic.jl

@@ -575,24 +575,21 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x))
 function generic_norm2(x)
     maxabs = normInf(x)
     (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs
+    return _generic_norm2(x, maxabs)
+end
+
+function _generic_norm2(x, maxabs::T) where {T}
     (v, s) = iterate(x)::Tuple
-    T = typeof(maxabs)
     if isfinite(length(x)*maxabs*maxabs) && !iszero(maxabs*maxabs) # Scaling not necessary
         sum::promote_type(Float64, T) = norm_sqr(v)
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             sum += norm_sqr(v)
         end
         ismissing(sum) && return missing
         return convert(T, sqrt(sum))
     else
         sum = abs2(norm(v)/maxabs)
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             sum += (norm(v)/maxabs)^2
         end
         ismissing(sum) && return missing
@@ -607,28 +604,27 @@ function generic_normp(x, p)
     if p > 1 || p < -1 # might need to rescale to avoid overflow
         maxabs = p > 1 ? normInf(x) : normMinusInf(x)
         (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs
-        T = typeof(maxabs)
+        return _generic_normp(x, p, maxabs)
     else
-        T = typeof(float(norm(v)))
+        # in this case, only the type of the last argument is used
+        # there is no scaling involved
+        return _generic_normp(x, p, float(norm(v)))
     end
+end
+
+function _generic_normp(x, p, maxabs::T) where {T}
     spp::promote_type(Float64, T) = p
     if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && !iszero(maxabs^spp)) # scaling not necessary
         sum::promote_type(Float64, T) = norm(v)^spp
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             ismissing(v) && return missing
             sum += norm(v)^spp
         end
         return convert(T, sum^inv(spp))
     else # rescaling
         sum = (norm(v)/maxabs)^spp
         ismissing(sum) && return missing
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             ismissing(v) && return missing
             sum += (norm(v)/maxabs)^spp
         end