Devectorize norm calculation and avoid code dublication.

base/linalg/blas.jl

@@ -161,14 +161,14 @@ for (fname, elty, ret_type) in ((:dnrm2_,:Float64,:Float64),
                                 (:scnrm2_,:Complex64,:Float32))
     @eval begin
         # SUBROUTINE DNRM2(N,X,INCX)
-        function nrm2(n::Integer, X::Union(Ptr{$elty},Array{$elty}), incx::Integer)
+        function nrm2(n::Integer, X::Union(Ptr{$elty},StridedVector{$elty}), incx::Integer)
             ccall(($(string(fname)),libblas), $ret_type,
                 (Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &n, X, &incx)
         end
     end
 end
-nrm2(A::Array) = nrm2(length(A), A, 1)
+nrm2(x::StridedVector) = nrm2(length(x), x, stride(x,1))
 
 ## asum
 for (fname, elty, ret_type) in ((:dasum_,:Float64,:Float64),
@@ -177,14 +177,14 @@ for (fname, elty, ret_type) in ((:dasum_,:Float64,:Float64),
                                 (:scasum_,:Complex64,:Float32))
      @eval begin
         # SUBROUTINE ASUM(N, X, INCX)
-        function asum(n::Integer, X::Union(Ptr{$elty},Array{$elty}), incx::Integer)
+        function asum(n::Integer, X::Union(Ptr{$elty},StridedVector{$elty}), incx::Integer)
             ccall(($(string(fname)),libblas), $ret_type,
                 (Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                  &n, X, &incx)
         end
     end
 end
-asum(A::Array) = asum(length(A), A, 1)
+asum(x::StridedVector) = asum(length(x), x, stride(x,1))
 
 ## axpy
 for (fname, elty) in ((:daxpy_,:Float64),

base/linalg/dense.jl

@@ -14,34 +14,16 @@ isposdef{T<:Number}(A::Matrix{T}, UL::Char) = isposdef!(float64(A), UL)
 isposdef{T<:Number}(A::Matrix{T}) = isposdef!(float64(A))
 isposdef(x::Number) = imag(x)==0 && real(x) > 0
 
-norm{T<:BlasFloat}(x::Vector{T}) = BLAS.nrm2(length(x), x, 1)
-
-function norm{T<:BlasFloat, TI<:Integer}(x::Vector{T}, rx::Union(Range1{TI},Range{TI}))
+function norm{T<:BlasFloat, TI<:Integer}(x::StridedVector{T}, rx::Union(Range1{TI},Range{TI}))
     (minimum(rx) < 1 || maximum(rx) > length(x)) && throw(BoundsError())
     BLAS.nrm2(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx))
 end
 
-function norm{T<:BlasFloat}(x::Vector{T}, p::Number)
-    n = length(x)
-    if n == 0
-        return zero(T)
-    elseif p == 2
-        return BLAS.nrm2(n, x, 1)
-    elseif p == 1
-        return BLAS.asum(n, x, 1)
-    elseif p == Inf
-        return maximum(abs(x))
-    elseif p == -Inf
-        return minimum(abs(x))
-    elseif p == 0
-        return convert(T, nnz(x))
-    else
-        absx = abs(x)
-        dx = maximum(absx)
-        dx==zero(T) && return zero(T)
-        scale!(absx, 1/dx)
-        return dx * (sum(absx.^p).^(1/p))
-    end
+function norm{T<:BlasFloat}(x::StridedVector{T}, p::Number=2)
+    length(x) == 0 && return zero(T)
+    p == 1 && T <: Real && return BLAS.asum(x)
+    p == 2 && return BLAS.nrm2(x)
+    invoke(norm, (AbstractVector, Number), x, p)
 end
 
 function triu!{T}(M::Matrix{T}, k::Integer)

base/linalg/generic.jl

@@ -43,59 +43,93 @@ diag(A::AbstractVector) = error("use diagm instead of diag to construct a diagon
 
 #diagm{T}(v::AbstractVecOrMat{T})
 
-function norm{T}(x::AbstractVector{T}, p::Number)
-    if length(x) == 0
-        a = zero(T)
-    elseif p == Inf
-        a = maximum(abs(x))
-    elseif p == -Inf
-        a = minimum(abs(x))
-    else
-        absx = abs(x)
-        dx = maximum(absx)
-        if dx != zero(T)
-            scale!(absx, 1/dx)
-            a = dx * (sum(absx.^p).^(1/p))
-        else
-            a = sum(absx.^p).^(1/p)
-        end
+function normMinusInf(x::AbstractVector)
+    n = length(x)
+    n > 0 || return real(zero(eltype(x)))
+    a = abs(x[1])
+    @inbounds for i = 2:n
+        a = min(a, abs(x[i]))
+    end
+    return a
+end
+function normInf(x::AbstractVector)
+    a = real(zero(eltype(x)))
+    @inbounds for i = 1:length(x)
+        a = max(a, abs(x[i]))
+    end
+    return a
+end
+function norm1(x::AbstractVector)
+    a = real(zero(eltype(x)))
+    @inbounds for i = 1:length(x)
+        a += abs(x[i])
+    end
+    return a
+end
+function norm2(x::AbstractVector)
+    nrmInfInv = inv(norm(x,Inf))
+    isinf(nrmInfInv) && return zero(nrmInfInv)
+    a = abs2(x[1]*nrmInfInv)
+    @inbounds for i = 2:length(x)
+        a += abs2(x[i]*nrmInfInv)
     end
-    float(a)
+    return sqrt(a)/nrmInfInv
+end
+function normp(x::AbstractVector,p::Number)
+    absx = convert(Vector{typeof(sqrt(abs(x[1])))}, abs(x))
+    dx = maximum(absx)
+    dx == 0 && return zero(typeof(absx))
+    scale!(absx, 1/dx)
+    pp = convert(eltype(absx), p)
+    return dx*sum(absx.^pp)^inv(pp)
+end
+function norm(x::AbstractVector, p::Number=2)
+    p == 0 && return nnz(x)
+    p == Inf && return normInf(x)
+    p == -Inf && return normMinusInf(x)
+    p == 1 && return norm1(x)
+    p == 2 && return norm2(x)
+    normp(x,p)
 end
-norm{T<:Integer}(x::AbstractVector{T}, p::Number) = norm(float(x), p)
-norm(x::AbstractVector) = norm(x, 2)
 
-function norm{T}(A::AbstractMatrix{T}, p::Number=2)
-    m, n = size(A)
-    if m == 0 || n == 0
-        return zero(real(zero(T)))
-    elseif m == 1 || n == 1
-        return norm(reshape(A, length(A)), p)
-    elseif p == 1
-        nrm = zero(real(zero(T)))
+function norm1{T}(A::AbstractMatrix{T})
+    m,n = size(A)
+    nrm = zero(real(zero(T)))
+    @inbounds begin
         for j = 1:n
             nrmj = zero(real(zero(T)))
             for i = 1:m
                 nrmj += abs(A[i,j])
             end
             nrm = max(nrm,nrmj)
         end
-        return nrm
-    elseif p == 2
-        return svdvals(A)[1]
-    elseif p == Inf
-        nrm = zero(real(zero(T)))
+    end
+    return nrm
+end
+function norm2(A::AbstractMatrix)
+    m,n = size(A)
+    if m == 0 || n == 0 return real(zero(eltype(A))) end
+    svdvals(A)[1]
+end
+function normInf{T}(A::AbstractMatrix{T})
+    m,n = size(A)
+    nrm = zero(real(zero(T)))
+    @inbounds begin
         for i = 1:m
             nrmi = zero(real(zero(T)))
             for j = 1:n
                 nrmi += abs(A[i,j])
             end
             nrm = max(nrm,nrmi)
         end
-        return nrm
-    else
-        throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
     end
+    return nrm
+end
+function norm{T}(A::AbstractMatrix{T}, p::Number=2)
+    p == 1 && return norm1(A)
+    p == 2 && return norm2(A)
+    p == Inf && return normInf(A)
+    throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf"))
 end
 
 norm(x::Number, p=nothing) = abs(x)

base/subarray.jl

@@ -403,6 +403,7 @@ stride(s::SubArray, i::Integer) = i <= length(s.strides) ? s.strides[i] : s.stri
 
 convert{T}(::Type{Ptr{T}}, x::SubArray{T}) =
     pointer(x.parent) + (x.first_index-1)*sizeof(T)
+convert{T,S,N}(::Type{Array{T,N}},A::SubArray{S,N}) = copy!(Array(T,size(A)), A)
 
 pointer(s::SubArray, i::Int) = pointer(s, ind2sub(size(s), i))
 

doc/stdlib/linalg.rst

@@ -270,9 +270,9 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: norm(A, [p])
 
-   Compute the ``p``-norm of a vector or a matrix ``A``, defaulting to the ``p=2``-norm.
+   Compute the ``p``-norm of a vector or the operator norm of a matrix ``A``, defaulting to the ``p=2``-norm.
 
-   For vectors, ``p`` can be any numeric value is valid (even though not all values produce a mathematically valid vector norm). In particular, `norm(A, Inf)`` returns the largest value in ``abs(A)``, whereas ``norm(A, -Inf)`` returns the smallest.
+   For vectors, ``p`` can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular, `norm(A, Inf)`` returns the largest value in ``abs(A)``, whereas ``norm(A, -Inf)`` returns the smallest.
 
    For matrices, valid values of ``p`` are ``1``, ``2``, or ``Inf``. Use :func:`normfro` to compute the Frobenius norm.
 

test/linalg.jl

@@ -908,6 +908,120 @@ for elty in (Float32, Float64, Complex64, Complex128)
     @test norm(F[:vectors]*Diagonal(F[:values])/F[:vectors] - A) > 0.01
 end
 
+# Tests norms
+nnorm = 1000
+mmat = 100
+nmat = 80
+for elty in (Float16, Float32, Float64, BigFloat, Complex{Float16}, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
+    debug && println(elty)
+    ## Vector
+    x = ones(elty,10)
+    xs = sub(x,1:2:10)
+    @test_approx_eq norm(x, -Inf) 1
+    @test_approx_eq norm(x, -1) 1/10
+    @test_approx_eq norm(x, 0) 10
+    @test_approx_eq norm(x, 1) 10
+    @test_approx_eq norm(x, 2) sqrt(10)
+    @test_approx_eq norm(x, 3) cbrt(10)
+    @test_approx_eq norm(x, Inf) 1
+    @test_approx_eq norm(xs, -Inf) 1
+    @test_approx_eq norm(xs, -1) 1/5
+    @test_approx_eq norm(xs, 0) 5
+    @test_approx_eq norm(xs, 1) 5
+    @test_approx_eq norm(xs, 2) sqrt(5)
+    @test_approx_eq norm(xs, 3) cbrt(5)
+    @test_approx_eq norm(xs, Inf) 1
+
+    for i = 1:10    
+        x = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : 
+            elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : 
+            convert(Vector{elty}, randn(nnorm))
+        xs = sub(x,1:2:nnorm)
+        y = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : 
+            elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : 
+            convert(Vector{elty}, randn(nnorm))        
+        ys = sub(y,1:2:nnorm)
+        α = elty <: Integer ? randn() : 
+            elty <: Complex ? convert(elty, complex(randn(),randn())) : 
+            convert(elty, randn())
+        # Absolute homogeneity
+        @test_approx_eq norm(α*x,-Inf) abs(α)*norm(x,-Inf)
+        @test_approx_eq norm(α*x,-1) abs(α)*norm(x,-1)
+        @test_approx_eq norm(α*x,1) abs(α)*norm(x,1)
+        @test_approx_eq norm(α*x) abs(α)*norm(x) # two is default
+        @test_approx_eq norm(α*x,3) abs(α)*norm(x,3)
+        @test_approx_eq norm(α*x,Inf) abs(α)*norm(x,Inf)
+
+        @test_approx_eq norm(α*xs,-Inf) abs(α)*norm(xs,-Inf)
+        @test_approx_eq norm(α*xs,-1) abs(α)*norm(xs,-1)
+        @test_approx_eq norm(α*xs,1) abs(α)*norm(xs,1)
+        @test_approx_eq norm(α*xs) abs(α)*norm(xs) # two is default
+        @test_approx_eq norm(α*xs,3) abs(α)*norm(xs,3)
+        @test_approx_eq norm(α*xs,Inf) abs(α)*norm(xs,Inf)
+
+        # Triangle inequality
+        @test norm(x + y,1) <= norm(x,1) + norm(y,1)
+        @test norm(x + y) <= norm(x) + norm(y) # two is default
+        @test norm(x + y,3) <= norm(x,3) + norm(y,3)
+        @test norm(x + y,Inf) <= norm(x,Inf) + norm(y,Inf)
+
+        @test norm(xs + ys,1) <= norm(xs,1) + norm(ys,1)
+        @test norm(xs + ys) <= norm(xs) + norm(ys) # two is default
+        @test norm(xs + ys,3) <= norm(xs,3) + norm(ys,3)
+        @test norm(xs + ys,Inf) <= norm(xs,Inf) + norm(ys,Inf)
+
+        # Against vectorized versions
+        @test_approx_eq norm(x,-Inf) minimum(abs(x))
+        @test_approx_eq norm(x,-1) inv(sum(1./abs(x)))
+        @test_approx_eq norm(x,0) sum(x .!= 0)
+        @test_approx_eq norm(x,1) sum(abs(x))
+        @test_approx_eq norm(x) sqrt(sum(abs2(x)))
+        @test_approx_eq norm(x,3) cbrt(sum(abs(x).^3.))
+        @test_approx_eq norm(x,Inf) maximum(abs(x))
+    end
+    ## Matrix (Operator)
+        A = ones(elty,10,10)
+        As = sub(A,1:5,1:5)
+        @test_approx_eq norm(A, 1) 10
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test_approx_eq norm(A, 2) 10
+        @test_approx_eq norm(A, Inf) 10
+        @test_approx_eq norm(As, 1) 5
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test_approx_eq norm(As, 2) 5
+        @test_approx_eq norm(As, Inf) 5
+
+    for i = 1:10
+        A = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : 
+            elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : 
+            convert(Matrix{elty}, randn(mmat, nmat))
+        As = sub(A,1:nmat,1:nmat)
+        B = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : 
+            elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : 
+            convert(Matrix{elty}, randn(mmat, nmat))        
+        Bs = sub(B,1:nmat,1:nmat)
+        α = elty <: Integer ? randn() :
+            elty <: Complex ? convert(elty, complex(randn(),randn())) : 
+            convert(elty, randn())
+
+        # Absolute homogeneity
+        @test_approx_eq norm(α*A,1) abs(α)*norm(A,1)
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test_approx_eq norm(α*A) abs(α)*norm(A) # two is default
+        @test_approx_eq norm(α*A,Inf) abs(α)*norm(A,Inf)
+
+        @test_approx_eq norm(α*As,1) abs(α)*norm(As,1)
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test_approx_eq norm(α*As) abs(α)*norm(As) # two is default
+        @test_approx_eq norm(α*As,Inf) abs(α)*norm(As,Inf)
+
+        # Triangle inequality
+        @test norm(A + B,1) <= norm(A,1) + norm(B,1)
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test norm(A + B) <= norm(A) + norm(B) # two is default
+        @test norm(A + B,Inf) <= norm(A,Inf) + norm(B,Inf)
+
+        @test norm(As + Bs,1) <= norm(As,1) + norm(Bs,1)
+        elty <: Union(BigFloat,Complex{BigFloat},BigInt) || @test norm(As + Bs) <= norm(As) + norm(Bs) # two is default
+        @test norm(As + Bs,Inf) <= norm(As,Inf) + norm(Bs,Inf)
+    end
+end
+
 ## Issue related tests
 # issue 1447
 let