Skip to content

Commit eaf6a55

Browse files
author
Andy Ferris
committed
Added more norm() functions.
Should be relatively performant for a range of (vector) p-norms. Possible warrants benchmarking?
1 parent 8db1be2 commit eaf6a55

File tree

3 files changed

+49
-2
lines changed

3 files changed

+49
-2
lines changed

src/StaticArrays.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ import Base: @pure, @propagate_inbounds, getindex, setindex!, size, similar,
66
length, convert, promote_op, map, map!, reduce, mapreduce,
77
broadcast, broadcast!, conj, transpose, ctranspose, hcat, vcat,
88
ones, zeros, eye, cross, vecdot, reshape, fill, fill!, det, inv,
9-
eig, trace, vecnorm, dot
9+
eig, trace, vecnorm, norm, dot
1010

1111
export StaticScalar, StaticArray, StaticVector, StaticMatrix
1212
export Scalar, SArray, SVector, SMatrix

src/linalg.jl

Lines changed: 42 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -285,9 +285,14 @@ end
285285
end
286286
end
287287

288+
@inline norm(v::StaticVector) = vecnorm(v)
289+
@inline norm(v::StaticVector, p::Real) = vecnorm(v, p)
290+
291+
@inline Base.LinAlg.norm_sqr(v::StaticVector) = mapreduce(abs2, +, zero(real(eltype(v))), v)
292+
288293
@generated function vecnorm(a::StaticArray)
289294
if length(a) == 0
290-
return zero(eltype(a))
295+
return zero(real(eltype(a)))
291296
end
292297

293298
expr = :(real(conj(a[1]) * a[1]))
@@ -301,7 +306,43 @@ end
301306
end
302307
end
303308

309+
@generated function vecnorm(a::StaticArray, p::Real)
310+
if length(a) == 0
311+
return zero(real(eltype(a)))
312+
end
313+
314+
expr = :(abs(a[1])^p)
315+
for j = 2:length(a)
316+
expr = :($expr + abs(a[$j])^p)
317+
end
318+
319+
expr_p1 = :(abs(a[1]))
320+
for j = 2:length(a)
321+
expr_p1 = :($expr_p1 + abs(a[$j]))
322+
end
323+
324+
return quote
325+
$(Expr(:meta, :inline))
326+
if p == Inf
327+
return mapreduce(abs, max, $(zero(real(eltype(a)))), a)
328+
elseif p == 1
329+
@inbounds return $expr_p1
330+
elseif p == 2
331+
return vecnorm(a)
332+
elseif p == 0
333+
return mapreduce(x -> (x == 0 ? zero(real(eltype(a))) : one(real(eltype(a)))), +, $(zero(real(eltype(a)))), a)
334+
else
335+
@inbounds return ($expr)^(inv(p))
336+
end
337+
end
338+
end
339+
304340
@inline Base.normalize(a::StaticVector) = inv(vecnorm(a))*a
341+
@inline Base.normalize(a::StaticVector, p::Real) = inv(vecnorm(a, p))*a
342+
343+
@inline Base.normalize!(a::StaticVector) = (a .*= inv(vecnorm(a)))
344+
@inline Base.normalize!(a::StaticVector, p::Real) = (a.*= inv(vecnorm(a, p)))
345+
305346

306347
@generated function trace(a::StaticMatrix)
307348
s = size(a)

test/linalg.jl

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,12 @@
5454
end
5555

5656
@testset "normalization" begin
57+
@test norm(SVector(1.0,2.0,2.0)) 3.0
58+
@test norm(SVector(1.0,2.0,2.0),2) 3.0
59+
@test norm(SVector(1.0,2.0,2.0),Inf) 2.0
60+
@test norm(SVector(1.0,2.0,2.0),1) 5.0
61+
@test norm(SVector(1.0,2.0,0.0),0) 2.0
62+
5763
@test vecnorm(SVector(1.0,2.0)) vecnorm([1.0,2.0])
5864
@test vecnorm(@SMatrix [1 2; 3 4.0+im]) vecnorm([1 2; 3 4.0+im])
5965

0 commit comments

Comments
 (0)