Skip to content

Commit

Permalink
recursive vecnorm
Browse files Browse the repository at this point in the history
  • Loading branch information
Jutho committed Dec 15, 2017
1 parent df2616e commit e634972
Showing 1 changed file with 29 additions and 24 deletions.
53 changes: 29 additions & 24 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -280,10 +280,10 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons
function generic_vecnormMinusInf(x)
s = start(x)
(v, s) = next(x, s)
minabs = norm(v)
minabs = vecnormMinusInf(v)
while !done(x, s)
(v, s) = next(x, s)
vnorm = norm(v)
vnorm = vecnormMinusInf(v)
minabs = ifelse(isnan(minabs) | (minabs < vnorm), minabs, vnorm)
end
return float(minabs)
Expand All @@ -292,10 +292,10 @@ end
function generic_vecnormInf(x)
s = start(x)
(v, s) = next(x, s)
maxabs = norm(v)
maxabs = vecnormInf(v)
while !done(x, s)
(v, s) = next(x, s)
vnorm = norm(v)
vnorm = vecnormInf(v)
maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm)
end
return float(maxabs)
Expand All @@ -304,20 +304,20 @@ end
function generic_vecnorm1(x)
s = start(x)
(v, s) = next(x, s)
av = float(norm(v))
av = float(vecnorm1(v))
T = typeof(av)
sum::promote_type(Float64, T) = av
while !done(x, s)
(v, s) = next(x, s)
sum += norm(v)
sum += vecnorm1(v)
end
return convert(T, sum)
end

# 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))
vecnorm_sqr(x) = vecnorm(x)^2
vecnorm_sqr(x::Number) = abs2(x)
vecnorm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x))

function generic_vecnorm2(x)
maxabs = vecnormInf(x)
Expand All @@ -326,17 +326,17 @@ function generic_vecnorm2(x)
(v, s) = next(x, s)
T = typeof(maxabs)
if isfinite(_length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
sum::promote_type(Float64, T) = norm_sqr(v)
sum::promote_type(Float64, T) = vecnorm_sqr(v)
while !done(x, s)
(v, s) = next(x, s)
sum += norm_sqr(v)
sum += vecnorm_sqr(v)
end
return convert(T, sqrt(sum))
else
sum = abs2(norm(v)/maxabs)
sum = abs2(vecnorm(v)/maxabs)
while !done(x, s)
(v, s) = next(x, s)
sum += (norm(v)/maxabs)^2
sum += (vecnorm(v)/maxabs)^2
end
return convert(T, maxabs*sqrt(sum))
end
Expand All @@ -352,21 +352,21 @@ function generic_vecnormp(x, p)
(maxabs == 0 || isinf(maxabs)) && return maxabs
T = typeof(maxabs)
else
T = typeof(float(norm(v)))
T = typeof(float(vecnorm(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
sum::promote_type(Float64, T) = vecnormp(v, p)^spp
while !done(x, s)
(v, s) = next(x, s)
sum += norm(v)^spp
sum += vecnormp(v, p)^spp
end
return convert(T, sum^inv(spp))
else # rescaling
sum = (norm(v)/maxabs)^spp
sum = (vecnormp(v, p)/maxabs)^spp
while !done(x, s)
(v, s) = next(x, s)
sum += (norm(v)/maxabs)^spp
sum += (vecnorm(v, p)/maxabs)^spp
end
return convert(T, maxabs*sum^inv(spp))
end
Expand Down Expand Up @@ -449,6 +449,11 @@ julia> vecnorm(-2, Inf)
```
"""
@inline vecnorm(x::Number, p::Real=2) = p == 0 ? (x==0 ? zero(abs(x)) : oneunit(abs(x))) : abs(x)
@inline vecnormMinusInf(x::Number) = abs(x)
@inline vecnormInf(x::Number) = abs(x)
@inline vecnorm1(x::Number) = abs(x)
@inline vecnorm2(x::Number) = abs(x)
@inline vecnormp(x::Number, p) = p == 0 ? (x==0 ? zero(abs(x)) : oneunit(abs(x))) : abs(x)

function norm1(A::AbstractMatrix{T}) where T
m, n = size(A)
Expand Down Expand Up @@ -621,11 +626,11 @@ function vecdot(x::AbstractArray, y::AbstractArray)
throw(DimensionMismatch("first array has length $(lx) which does not match the length of the second, $(_length(y))."))
end
if lx == 0
return dot(zero(eltype(x)), zero(eltype(y)))
return vecdot(zero(eltype(x)), zero(eltype(y)))
end
s = zero(dot(first(x), first(y)))
s = zero(vecdot(first(x), first(y)))
for (Ix, Iy) in zip(eachindex(x), eachindex(y))
@inbounds s += dot(x[Ix], y[Iy])
@inbounds s += vecdot(x[Ix], y[Iy])
end
s
end
Expand Down Expand Up @@ -656,22 +661,22 @@ function vecdot(x, y) # arbitrary iterables
if !isempty(y)
throw(DimensionMismatch("x and y are of different lengths!"))
end
return dot(zero(eltype(x)), zero(eltype(y)))
return vecdot(zero(eltype(x)), zero(eltype(y)))
end
iy = start(y)
if done(y, iy)
throw(DimensionMismatch("x and y are of different lengths!"))
end
(vx, ix) = next(x, ix)
(vy, iy) = next(y, iy)
s = dot(vx, vy)
s = vecdot(vx, vy)
while !done(x, ix)
if done(y, iy)
throw(DimensionMismatch("x and y are of different lengths!"))
end
(vx, ix) = next(x, ix)
(vy, iy) = next(y, iy)
s += dot(vx, vy)
s += vecdot(vx, vy)
end
if !done(y, iy)
throw(DimensionMismatch("x and y are of different lengths!"))
Expand Down

0 comments on commit e634972

Please sign in to comment.