From cb2c2721907840aaad8f024130c9203042a37e66 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 3 Jun 2018 14:26:48 +0200 Subject: [PATCH 01/17] change vecnorm to norm; introduce opnorm for matrix norms --- stdlib/LinearAlgebra/docs/src/index.md | 2 +- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 2 +- stdlib/LinearAlgebra/src/dense.jl | 14 +- stdlib/LinearAlgebra/src/deprecated.jl | 7 +- stdlib/LinearAlgebra/src/generic.jl | 240 +++++++++++--------- stdlib/LinearAlgebra/src/triangular.jl | 22 +- stdlib/LinearAlgebra/src/uniformscaling.jl | 10 +- stdlib/LinearAlgebra/test/adjtrans.jl | 25 +- stdlib/LinearAlgebra/test/bidiag.jl | 2 +- stdlib/LinearAlgebra/test/dense.jl | 24 +- stdlib/LinearAlgebra/test/generic.jl | 8 +- stdlib/LinearAlgebra/test/matmul.jl | 2 +- stdlib/LinearAlgebra/test/pinv.jl | 8 +- stdlib/LinearAlgebra/test/qr.jl | 2 +- stdlib/LinearAlgebra/test/triangular.jl | 2 +- stdlib/LinearAlgebra/test/uniformscaling.jl | 2 +- stdlib/SparseArrays/src/SparseArrays.jl | 4 +- stdlib/SparseArrays/src/linalg.jl | 20 +- stdlib/SparseArrays/src/sparsevector.jl | 2 +- stdlib/SparseArrays/test/sparse.jl | 64 +++--- stdlib/SparseArrays/test/sparsevector.jl | 10 +- stdlib/SuiteSparse/test/cholmod.jl | 4 +- 22 files changed, 253 insertions(+), 223 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 3b3d6c95d8a0b..f8acc46722bef 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -358,7 +358,7 @@ LinearAlgebra.diag LinearAlgebra.diagm LinearAlgebra.rank LinearAlgebra.norm -LinearAlgebra.vecnorm +LinearAlgebra.opnorm LinearAlgebra.normalize! LinearAlgebra.normalize LinearAlgebra.cond diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 69fe8cb3139d8..6d4e082914350 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -126,6 +126,7 @@ export qr!, lq, lq!, + opnorm, rank, rdiv!, schur, @@ -144,7 +145,6 @@ export tril!, triu!, vecdot, - vecnorm, # Operators \, diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index f670bf5e59b61..9c121abe26514 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -137,11 +137,11 @@ function norm(x::StridedVector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}) w GC.@preserve x BLAS.nrm2(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx)) end -vecnorm1(x::Union{Array{T},StridedVector{T}}) where {T<:BlasReal} = - length(x) < ASUM_CUTOFF ? generic_vecnorm1(x) : BLAS.asum(x) +norm1(x::Union{Array{T},StridedVector{T}}) where {T<:BlasReal} = + length(x) < ASUM_CUTOFF ? generic_norm1(x) : BLAS.asum(x) -vecnorm2(x::Union{Array{T},StridedVector{T}}) where {T<:BlasFloat} = - length(x) < NRM2_CUTOFF ? generic_vecnorm2(x) : BLAS.nrm2(x) +norm2(x::Union{Array{T},StridedVector{T}}) where {T<:BlasFloat} = + length(x) < NRM2_CUTOFF ? generic_norm2(x) : BLAS.nrm2(x) """ triu!(M, k::Integer) @@ -509,7 +509,7 @@ function exp!(A::StridedMatrix{T}) where T<:BlasFloat return copytri!(parent(exp(Hermitian(A))), 'U', true) end ilo, ihi, scale = LAPACK.gebal!('B', A) # modifies A - nA = norm(A, 1) + nA = opnorm(A, 1) Inn = Matrix{T}(I, n, n) ## For sufficiently small nA, use lower order Padé-Approximations if (nA <= 2.1) @@ -1370,8 +1370,8 @@ function cond(A::AbstractMatrix, p::Real=2) end throw(ArgumentError("p-norm must be 1, 2 or Inf, got $p")) end -_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A), p, norm(A, p)) -_cond1Inf(A::AbstractMatrix, p::Real) = norm(A, p)*norm(inv(A), p) +_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lu(A), p, opnorm(A, p)) +_cond1Inf(A::AbstractMatrix, p::Real) = opnorm(A, p)*opnorm(inv(A), p) ## Lyapunov and Sylvester equation diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index f79afa785b553..4ed5e5f8a2734 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -6,6 +6,9 @@ using Base: @deprecate, depwarn @deprecate cond(F::LinearAlgebra.LU, p::Integer) cond(convert(AbstractArray, F), p) +# deprecate vecnorm in favor of norm +@deprecate vecnorm norm + # PR #22188 export cholfact, cholfact! @deprecate cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) cholesky!(Hermitian(A, uplo), Val(false)) @@ -1201,8 +1204,8 @@ _mat_ldiv_rowvec_error() = throw(DimensionMismatch("Cannot left-divide matrix by *(A::RowVector, B::Adjoint{<:Any,<:AbstractRotation}) = A * adjoint(B.parent) # methods involving RowVector from base/linalg/generic.jl, to deprecate -norm(tv::RowVector, q::Real) = q == Inf ? norm(rvtranspose(tv), 1) : norm(rvtranspose(tv), q/(q-1)) -norm(tv::RowVector) = norm(rvtranspose(tv)) +norm(tv::RowVector, q::Real) = q == Inf ? opnorm(rvtranspose(tv), 1) : opnorm(rvtranspose(tv), q/(q-1)) +norm(tv::RowVector) = opnorm(rvtranspose(tv)) # methods involving RowVector from base/linalg/factorization.jl, to deprecate \(A::Adjoint{<:Any,<:Factorization}, B::RowVector) = adjoint(A.parent) \ B diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 78bf499553ef0..7f76a32ed9016 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -249,43 +249,43 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons ########################################################################################### # Inner products and norms -# special cases of vecnorm; note that they don't need to handle isempty(x) -function generic_vecnormMinusInf(x) +# special cases of norm; note that they don't need to handle isempty(x) +function generic_normMinusInf(x) (v, s) = iterate(x)::Tuple - minabs = norm(v) + minabs = norm(v, -Inf) while true y = iterate(x, s) y === nothing && break (v, s) = y - vnorm = norm(v) + vnorm = norm(v, -Inf) minabs = ifelse(isnan(minabs) | (minabs < vnorm), minabs, vnorm) end return float(minabs) end -function generic_vecnormInf(x) +function generic_normInf(x) (v, s) = iterate(x)::Tuple - maxabs = norm(v) + maxabs = norm(v, Inf) while true y = iterate(x, s) y === nothing && break (v, s) = y - vnorm = norm(v) + vnorm = norm(v, Inf) maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm) end return float(maxabs) end -function generic_vecnorm1(x) +function generic_norm1(x) (v, s) = iterate(x)::Tuple - av = float(norm(v)) + av = float(norm(v, 1)) T = typeof(av) sum::promote_type(Float64, T) = av while true y = iterate(x, s) y === nothing && break (v, s) = y - sum += norm(v) + sum += norm(v, 1) end return convert(T, sum) end @@ -295,8 +295,8 @@ 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)) -function generic_vecnorm2(x) - maxabs = vecnormInf(x) +function generic_norm2(x) + maxabs = normInf(x) (maxabs == 0 || isinf(maxabs)) && return maxabs (v, s) = iterate(x)::Tuple T = typeof(maxabs) @@ -323,45 +323,46 @@ end # Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p) # (Not technically a "norm" for p < 1.) -function generic_vecnormp(x, p) +function generic_normp(x, p) (v, s) = iterate(x)::Tuple if p > 1 || p < -1 # might need to rescale to avoid overflow - maxabs = p > 1 ? vecnormInf(x) : vecnormMinusInf(x) + maxabs = p > 1 ? normInf(x) : normMinusInf(x) (maxabs == 0 || isinf(maxabs)) && return maxabs T = typeof(maxabs) else - T = typeof(float(norm(v))) + T = typeof(float(norm(v,p))) 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) = norm(v,p)^spp while true y = iterate(x, s) y === nothing && break (v, s) = y - sum += norm(v)^spp + sum += norm(v,p)^spp end return convert(T, sum^inv(spp)) else # rescaling - sum = (norm(v)/maxabs)^spp + sum = (norm(v,p)/maxabs)^spp while true y = iterate(x, s) y == nothing && break (v, s) = y - sum += (norm(v)/maxabs)^spp + sum += (norm(v,p)/maxabs)^spp end return convert(T, maxabs*sum^inv(spp)) end end -vecnormMinusInf(x) = generic_vecnormMinusInf(x) -vecnormInf(x) = generic_vecnormInf(x) -vecnorm1(x) = generic_vecnorm1(x) -vecnorm2(x) = generic_vecnorm2(x) -vecnormp(x, p) = generic_vecnormp(x, p) +normMinusInf(x) = generic_normMinusInf(x) +normInf(x) = generic_normInf(x) +norm1(x) = generic_norm1(x) +norm2(x) = generic_norm2(x) +normp(x, p) = generic_normp(x, p) + """ - vecnorm(A, p::Real=2) + norm(A, p::Real=2) For any iterable container `A` (including arrays of any dimension) of numbers (or any element type for which `norm` is defined), compute the `p`-norm (defaulting to `p=2`) as if @@ -374,65 +375,96 @@ The `p`-norm is defined as: with ``a_i`` the entries of ``A`` and ``n`` its length. `p` can assume any numeric value (even though not all values produce a -mathematically valid vector norm). In particular, `vecnorm(A, Inf)` returns the largest value -in `abs(A)`, whereas `vecnorm(A, -Inf)` returns the smallest. If `A` is a matrix and `p=2`, +mathematically valid vector norm). In particular, `norm(A, Inf)` returns the largest value +in `abs.(A)`, whereas `norm(A, -Inf)` returns the smallest. If `A` is a matrix and `p=2`, then this is equivalent to the Frobenius norm. +Use [`opnorm`](@ref) to compute the operator norm of a matrix. + # Examples ```jldoctest -julia> vecnorm([1 2 3; 4 5 6; 7 8 9]) +julia> v = [3, -2, 6] +3-element Array{Int64,1}: + 3 + -2 + 6 + +julia> norm(v) +7.0 + +julia> norm(v, 1) +11.0 + +julia> norm(v, Inf) +6.0 + +julia> norm([1 2 3; 4 5 6; 7 8 9]) 16.881943016134134 -julia> vecnorm([1 2 3 4 5 6 7 8 9]) +julia> norm([1 2 3 4 5 6 7 8 9]) 16.881943016134134 + +julia> norm([1 2 3 4 5 6 7 8 9]') +16.881943016134134 + +julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) == norm([v,v], 1) +true + +julia> norm(hcat(v,v), 2) == norm(vcat(v,v), 2) == norm([v,v], 2) +true + +julia> norm(hcat(v,v), Inf) == norm(vcat(v,v), Inf) == norm([v,v], Inf) +true ``` """ -function vecnorm(itr, p::Real=2) +function norm(itr, p::Real=2) isempty(itr) && return float(norm(zero(eltype(itr)))) if p == 2 - return vecnorm2(itr) + return norm2(itr) elseif p == 1 - return vecnorm1(itr) + return norm1(itr) elseif p == Inf - return vecnormInf(itr) + return normInf(itr) elseif p == 0 return typeof(float(norm(first(itr))))(count(!iszero, itr)) elseif p == -Inf - return vecnormMinusInf(itr) + return normMinusInf(itr) else - vecnormp(itr,p) + normp(itr,p) end end """ - vecnorm(x::Number, p::Real=2) + norm(x::Number, p::Real=2) -For numbers, return ``\\left( |x|^p \\right) ^{1/p}``. +For numbers, return ``\\left( |x|^p \\right)^{1/p}``. # Examples ```jldoctest -julia> vecnorm(2, 1) +julia> norm(2, 1) 2 -julia> vecnorm(-2, 1) +julia> norm(-2, 1) 2 -julia> vecnorm(2, 2) +julia> norm(2, 2) 2 -julia> vecnorm(-2, 2) +julia> norm(-2, 2) 2 -julia> vecnorm(2, Inf) +julia> norm(2, Inf) 2 -julia> vecnorm(-2, Inf) +julia> norm(-2, Inf) 2 ``` """ -@inline vecnorm(x::Number, p::Real=2) = p == 0 ? (x==0 ? zero(abs(x)) : oneunit(abs(x))) : abs(x) +@inline norm(x::Number, p::Real=2) = p == 0 ? (x==0 ? zero(abs(x)) : oneunit(abs(x))) : abs(x) -function norm1(A::AbstractMatrix{T}) where T + +# special cases of opnorm +function opnorm1(A::AbstractMatrix{T}) where T m, n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) @@ -448,13 +480,15 @@ function norm1(A::AbstractMatrix{T}) where T end return convert(Tnorm, nrm) end -function norm2(A::AbstractMatrix{T}) where T + +function opnorm2(A::AbstractMatrix{T}) where T m,n = size(A) - if m == 1 || n == 1 return vecnorm2(A) end + if m == 1 || n == 1 return norm2(A) end Tnorm = typeof(float(real(zero(T)))) (m == 0 || n == 0) ? zero(Tnorm) : convert(Tnorm, svdvals(A)[1]) end -function normInf(A::AbstractMatrix{T}) where T + +function opnormInf(A::AbstractMatrix{T}) where T m,n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) @@ -471,58 +505,25 @@ function normInf(A::AbstractMatrix{T}) where T return convert(Tnorm, nrm) end -""" - norm(A::AbstractArray, p::Real=2) - -Compute the `p`-norm of a vector or the operator norm of a matrix `A`, -defaulting to the 2-norm. - - norm(A::AbstractVector, p::Real=2) - -For vectors, this is equivalent to [`vecnorm`](@ref) and equal to: -```math -\\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p} -``` -with ``a_i`` the entries of ``A`` and ``n`` its length. - -`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. - -# Examples -```jldoctest -julia> v = [3, -2, 6] -3-element Array{Int64,1}: - 3 - -2 - 6 - -julia> norm(v) -7.0 -julia> norm(v, Inf) -6.0 -``` """ -norm(x::AbstractVector, p::Real=2) = vecnorm(x, p) + opnorm(A::AbstractMatrix, p::Real=2) -""" - norm(A::AbstractMatrix, p::Real=2) +Computes the operator norm (or matrix norm) induced by the vector `p`-norm, +where valid values of `p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, +`p=2` is currently not implemented.) Use [`norm`](@ref) to compute the Frobenius +norm. -For matrices, the matrix norm induced by the vector `p`-norm is used, where valid values of -`p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, `p=2` is currently not -implemented.) Use [`vecnorm`](@ref) to compute the Frobenius norm. - -When `p=1`, the matrix norm is the maximum absolute column sum of `A`: +When `p=1`, the operator norm is the maximum absolute column sum of `A`: ```math \\|A\\|_1 = \\max_{1 ≤ j ≤ n} \\sum_{i=1}^m | a_{ij} | ``` with ``a_{ij}`` the entries of ``A``, and ``m`` and ``n`` its dimensions. -When `p=2`, the matrix norm is the spectral norm, equal to the largest +When `p=2`, the operator norm is the spectral norm, equal to the largest singular value of `A`. -When `p=Inf`, the matrix norm is the maximum absolute row sum of `A`: +When `p=Inf`, the operator norm is the maximum absolute row sum of `A`: ```math \\|A\\|_\\infty = \\max_{1 ≤ i ≤ m} \\sum _{j=1}^n | a_{ij} | ``` @@ -534,40 +535,44 @@ julia> A = [1 -2 -3; 2 3 -1] 1 -2 -3 2 3 -1 -julia> norm(A, Inf) +julia> opnorm(A, Inf) 6.0 + +julia> opnorm(A, 1) +5.0 ``` """ -function norm(A::AbstractMatrix, p::Real=2) +function opnorm(A::AbstractMatrix, p::Real=2) if p == 2 - return norm2(A) + return opnorm2(A) elseif p == 1 - return norm1(A) + return opnorm1(A) elseif p == Inf - return normInf(A) + return opnormInf(A) else throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf")) end end """ - norm(x::Number, p::Real=2) + opnorm(x::Number, p::Real=2) For numbers, return ``\\left( |x|^p \\right)^{1/p}``. -This is equivalent to [`vecnorm`](@ref). +This is equivalent to [`norm`](@ref). """ -@inline norm(x::Number, p::Real=2) = vecnorm(x, p) +@inline opnorm(x::Number, p::Real=2) = norm(x, p) """ - norm(A::Adjoint{<:Any,<:AbstracVector}, q::Real=2) - norm(A::Transpose{<:Any,<:AbstracVector}, q::Real=2) + opnorm(A::Adjoint{<:Any,<:AbstracVector}, q::Real=2) + opnorm(A::Transpose{<:Any,<:AbstracVector}, q::Real=2) -For Adjoint/Transpose-wrapped vectors, return the ``q``-norm of `A`, which is -equivalent to the p-norm with value `p = q/(q-1)`. They coincide at `p = q = 2`. +For Adjoint/Transpose-wrapped vectors, return the operator ``q``-norm of `A`, which is +equivalent to the `p`-norm with value `p = q/(q-1)`. They coincide at `p = q = 2`. +Use [`norm`](@ref) to compute the `p` norm of `A` as a vector. The difference in norm between a vector space and its dual arises to preserve the relationship between duality and the inner product, and the result is -consistent with the p-norm of `1 × n` matrix. +consistent with the operator `p`-norm of `1 × n` matrix. # Examples ```jldoctest @@ -575,29 +580,40 @@ julia> v = [1; im]; julia> vc = v'; -julia> norm(vc, 1) +julia> opnorm(vc, 1) 1.0 +julia> norm(vc, 1) +2.0 + julia> norm(v, 1) 2.0 +julia> opnorm(vc, 2) +1.4142135623730951 + julia> norm(vc, 2) 1.4142135623730951 julia> norm(v, 2) 1.4142135623730951 -julia> norm(vc, Inf) +julia> opnorm(vc, Inf) 2.0 +julia> norm(vc, Inf) +1.0 + julia> norm(v, Inf) 1.0 ``` """ -norm(v::TransposeAbsVec, q::Real) = q == Inf ? norm(v.parent, 1) : norm(v.parent, q/(q-1)) -norm(v::AdjointAbsVec, q::Real) = q == Inf ? norm(conj(v.parent), 1) : norm(conj(v.parent), q/(q-1)) -norm(v::AdjointAbsVec) = norm(conj(v.parent)) -norm(v::TransposeAbsVec) = norm(v.parent) +opnorm(v::TransposeAbsVec, q::Real) = q == Inf ? norm(v.parent, 1) : norm(v.parent, q/(q-1)) +opnorm(v::AdjointAbsVec, q::Real) = q == Inf ? norm(conj(v.parent), 1) : norm(conj(v.parent), q/(q-1)) +opnorm(v::AdjointAbsVec) = norm(conj(v.parent)) +opnorm(v::TransposeAbsVec) = norm(v.parent) + +norm(v::Union{TransposeAbsVec,AdjointAbsVec}, p::Real) = norm(v.parent, p) function vecdot(x::AbstractArray, y::AbstractArray) lx = _length(x) @@ -877,7 +893,7 @@ cond(x::Number) = x == 0 ? Inf : 1.0 cond(x::Number, p) = cond(x) #Skeel condition numbers -condskeel(A::AbstractMatrix, p::Real=Inf) = norm(abs.(inv(A))*abs.(A), p) +condskeel(A::AbstractMatrix, p::Real=Inf) = opnorm(abs.(inv(A))*abs.(A), p) """ condskeel(M, [x, p::Real=Inf]) @@ -1326,7 +1342,7 @@ promote_leaf_eltypes(x::Union{AbstractArray,Tuple}) = mapreduce(promote_leaf_elt function isapprox(x::AbstractArray, y::AbstractArray; atol::Real=0, rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y),atol), - nans::Bool=false, norm::Function=vecnorm) + nans::Bool=false, norm::Function=norm) d = norm(x - y) if isfinite(d) return d <= max(atol, rtol*max(norm(x), norm(y))) @@ -1341,7 +1357,7 @@ end Normalize the vector `v` in-place so that its `p`-norm equals unity, i.e. `norm(v, p) == 1`. -See also [`normalize`](@ref) and [`vecnorm`](@ref). +See also [`normalize`](@ref) and [`norm`](@ref). """ function normalize!(v::AbstractVector, p::Real=2) nrm = norm(v, p) @@ -1370,7 +1386,7 @@ end Normalize the vector `v` so that its `p`-norm equals unity, i.e. `norm(v, p) == vecnorm(v, p) == 1`. -See also [`normalize!`](@ref) and [`vecnorm`](@ref). +See also [`normalize!`](@ref) and [`norm`](@ref). # Examples ```jldoctest diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index eb21fefd8e33d..029501cc49bd2 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -1952,7 +1952,7 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real) ArgumentError("p must be a real number in (-1,1), got $p") end - normA0 = norm(A0, 1) + normA0 = opnorm(A0, 1) rmul!(A0, 1/normA0) theta = [1.53e-5, 2.25e-3, 1.92e-2, 6.08e-2, 1.25e-1, 2.03e-1, 2.84e-1] @@ -2050,8 +2050,8 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat end AmI = A - I - d2 = sqrt(norm(AmI^2, 1)) - d3 = cbrt(norm(AmI^3, 1)) + d2 = sqrt(opnorm(AmI^2, 1)) + d3 = cbrt(opnorm(AmI^3, 1)) alpha2 = max(d2, d3) foundm = false if alpha2 <= theta[2] @@ -2062,9 +2062,9 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat while !foundm more = false if s > s0 - d3 = cbrt(norm(AmI^3, 1)) + d3 = cbrt(opnorm(AmI^3, 1)) end - d4 = norm(AmI^4, 1)^(1/4) + d4 = opnorm(AmI^4, 1)^(1/4) alpha3 = max(d3, d4) if alpha3 <= theta[tmax] local j @@ -2083,7 +2083,7 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat end if !more - d5 = norm(AmI^5, 1)^(1/5) + d5 = opnorm(AmI^5, 1)^(1/5) alpha4 = max(d4, d5) eta = min(alpha3, alpha4) if eta <= theta[tmax] @@ -2250,8 +2250,8 @@ function invsquaring(A0::UpperTriangular, theta) end AmI = A - I - d2 = sqrt(norm(AmI^2, 1)) - d3 = cbrt(norm(AmI^3, 1)) + d2 = sqrt(opnorm(AmI^2, 1)) + d3 = cbrt(opnorm(AmI^3, 1)) alpha2 = max(d2, d3) foundm = false if alpha2 <= theta[2] @@ -2262,9 +2262,9 @@ function invsquaring(A0::UpperTriangular, theta) while !foundm more = false if s > s0 - d3 = cbrt(norm(AmI^3, 1)) + d3 = cbrt(opnorm(AmI^3, 1)) end - d4 = norm(AmI^4, 1)^(1/4) + d4 = opnorm(AmI^4, 1)^(1/4) alpha3 = max(d3, d4) if alpha3 <= theta[tmax] local j @@ -2287,7 +2287,7 @@ function invsquaring(A0::UpperTriangular, theta) end if !more - d5 = norm(AmI^5, 1)^(1/5) + d5 = opnorm(AmI^5, 1)^(1/5) alpha4 = max(d4, d5) eta = min(alpha3, alpha4) if eta <= theta[tmax] diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index d9861d373cdb4..01cba86974875 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -128,7 +128,7 @@ function (-)(J::UniformScaling, A::AbstractMatrix) end inv(J::UniformScaling) = UniformScaling(inv(J.λ)) -norm(J::UniformScaling, p::Real=2) = abs(J.λ) +opnorm(J::UniformScaling, p::Real=2) = opnorm(J.λ, p) function det(J::UniformScaling{T}) where T if isone(J.λ) @@ -194,11 +194,11 @@ end function isapprox(J::UniformScaling, A::AbstractMatrix; atol::Real = 0, rtol::Real = Base.rtoldefault(promote_leaf_eltypes(A), eltype(J), atol), - nans::Bool = false, norm::Function = vecnorm) + nans::Bool = false, norm::Function = norm) n = checksquare(A) - normJ = norm === LinearAlgebra.norm ? abs(J.λ) : - norm === vecnorm ? abs(J.λ) * sqrt(n) : - norm(Diagonal(fill(J.λ, n))) + normJ = norm === opnorm ? abs(J.λ) : + norm === norm ? abs(J.λ) * sqrt(n) : + norm(Diagonal(fill(J.λ, n))) return norm(A - J) <= max(atol, rtol * max(norm(A), normJ)) end isapprox(A::AbstractMatrix, J::UniformScaling; kwargs...) = isapprox(J, A; kwargs...) diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index f52b2bbb248d0..ae140d5f7a8c8 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -413,19 +413,30 @@ end @test (Transpose(complexvec) / Transpose(complexmat))::Transpose ≈ rowcomplexvec / copy(Transpose(complexmat)) end -@testset "norm of Adjoint/Transpose-wrapped vectors" begin +@testset "norm and opnorm of Adjoint/Transpose-wrapped vectors" begin # definitions are in base/linalg/generic.jl realvec, complexvec = [3, -4], [3im, -4im] - # one norm result should be maximum(abs.(realvec)) == 4 + # one norm result should be sum(abs.(realvec)) == 7 # two norm result should be sqrt(sum(abs.(realvec))) == 5 - # inf norm result should be sum(abs.(realvec)) == 7 + # inf norm result should be maximum(abs.(realvec)) == 4 for v in (realvec, complexvec) @test norm(Adjoint(v)) ≈ 5 - @test norm(Adjoint(v), 1) ≈ 4 - @test norm(Adjoint(v), Inf) ≈ 7 + @test norm(Adjoint(v), 1) ≈ 7 + @test norm(Adjoint(v), Inf) ≈ 4 @test norm(Transpose(v)) ≈ 5 - @test norm(Transpose(v), 1) ≈ 4 - @test norm(Transpose(v), Inf) ≈ 7 + @test norm(Transpose(v), 1) ≈ 7 + @test norm(Transpose(v), Inf) ≈ 4 + end + # one opnorm result should be maximum(abs.(realvec)) == 4 + # two opnorm result should be sqrt(sum(abs.(realvec))) == 5 + # inf opnorm result should be sum(abs.(realvec)) == 7 + for v in (realvec, complexvec) + @test opnorm(Adjoint(v)) ≈ 5 + @test opnorm(Adjoint(v), 1) ≈ 4 + @test opnorm(Adjoint(v), Inf) ≈ 7 + @test opnorm(Transpose(v)) ≈ 5 + @test opnorm(Transpose(v), 1) ≈ 4 + @test opnorm(Transpose(v), Inf) ≈ 7 end end diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index 4e4b2da53c301..b354d8b6fe2ae 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -251,7 +251,7 @@ srand(1) test_approx_eq_modphase(u1, u2) test_approx_eq_modphase(copy(v1), copy(v2)) end - @test 0 ≈ vecnorm(u2*Diagonal(d2)*v2'-Tfull) atol=n*max(n^2*eps(relty),vecnorm(u1*Diagonal(d1)*v1'-Tfull)) + @test 0 ≈ norm(u2*Diagonal(d2)*v2'-Tfull) atol=n*max(n^2*eps(relty),norm(u1*Diagonal(d1)*v1'-Tfull)) @inferred svdvals(T) @inferred svd(T) end diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index be35e774ac8d9..21aaff8f0325c 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -236,18 +236,18 @@ end end end - @testset "Matrix (Operator)" begin + @testset "Matrix (Operator) opnorm" begin A = fill(elty(1),10,10) As = view(A,1:5,1:5) - @test norm(A, 1) ≈ 10 - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A, 2) ≈ 10 - @test norm(A, Inf) ≈ 10 - @test norm(As, 1) ≈ 5 - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(As, 2) ≈ 5 - @test norm(As, Inf) ≈ 5 + @test opnorm(A, 1) ≈ 10 + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test opnorm(A, 2) ≈ 10 + @test opnorm(A, Inf) ≈ 10 + @test opnorm(As, 1) ≈ 5 + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test opnorm(As, 2) ≈ 5 + @test opnorm(As, Inf) ≈ 5 end - @testset "Absolute homogeneity, triangle inequality, & vecnorm" begin + @testset "Absolute homogeneity, triangle inequality, & norm" begin for i = 1:10 Ainit = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : elty <: Complex ? convert(Matrix{elty}, complex.(randn(mmat, nmat), randn(mmat, nmat))) : @@ -269,9 +269,9 @@ end 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) - # vecnorm: - for p = -2:3 - @test norm(reshape(A, length(A)), p) == vecnorm(A, p) + # norm + for p in (-Inf, Inf, (-2:3)...) + @test norm(A, p) == norm(vec(A), p) end end end @@ -577,7 +577,7 @@ end @test coth(acoth(coth(A))) ≈ coth(A) # Definition of principal values (Aprahamian & Higham, 2016, pp. 4-5) - abstol = sqrt(eps(real(elty))) * vecnorm(acosh(A)) + abstol = sqrt(eps(real(elty))) * norm(acosh(A)) @test all(z -> (0 < real(z) < π || abs(real(z)) < abstol && imag(z) >= 0 || abs(real(z) - π) < abstol && imag(z) <= 0), diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index f56a95e0b9211..58a85357054bc 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -197,14 +197,14 @@ end @test norm([2.4e-322, 4.4e-323]) ≈ 2.47e-322 @test norm([2.4e-322, 4.4e-323], 3) ≈ 2.4e-322 -@test_throws ArgumentError norm(Matrix{Float64}(undef,5,5),5) +@test_throws ArgumentError opnorm(Matrix{Float64}(undef,5,5),5) -@testset "generic vecnorm for arrays of arrays" begin +@testset "generic norm for arrays of arrays" begin x = Vector{Int}[[1,2], [3,4]] @test @inferred(norm(x)) ≈ sqrt(30) @test norm(x, 0) == length(x) - @test norm(x, 1) ≈ sqrt(5) + 5 - @test norm(x, 3) ≈ cbrt(sqrt(125)+125) + @test norm(x, 1) ≈ 10 + @test norm(x, 3) ≈ cbrt(100) end @testset "LinearAlgebra.axp(b)y! for element type without commutative multiplication" begin diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 171265dfbd04c..5146e29d72c9a 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -244,7 +244,7 @@ vecdot_(x,y) = invoke(vecdot, Tuple{Any,Any}, x,y) if n1 != n2 @test_throws DimensionMismatch d(1:n1, 1:n2) else - @test d(1:n1, 1:n2) ≈ vecnorm(1:n1)^2 + @test d(1:n1, 1:n2) ≈ norm(1:n1)^2 end end end diff --git a/stdlib/LinearAlgebra/test/pinv.jl b/stdlib/LinearAlgebra/test/pinv.jl index 42da2e1b7b791..489f7c0c8e421 100644 --- a/stdlib/LinearAlgebra/test/pinv.jl +++ b/stdlib/LinearAlgebra/test/pinv.jl @@ -89,14 +89,14 @@ end function test_pinv(a,m,n,tol1,tol2,tol3) apinv = @inferred pinv(a) - @test vecnorm(a*apinv*a-a)/vecnorm(a) ≈ 0 atol=tol1 + @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol1 x0 = randn(n); b = a*x0; x = apinv*b - @test vecnorm(a*x-b)/vecnorm(b) ≈ 0 atol=tol1 + @test norm(a*x-b)/norm(b) ≈ 0 atol=tol1 apinv = pinv(a,sqrt(eps(real(one(eltype(a)))))) - @test vecnorm(a*apinv*a-a)/vecnorm(a) ≈ 0 atol=tol2 + @test norm(a*apinv*a-a)/norm(a) ≈ 0 atol=tol2 x0 = randn(n); b = a*x0; x = apinv*b - @test vecnorm(a*x-b)/vecnorm(b) ≈ 0 atol=tol2 + @test norm(a*x-b)/norm(b) ≈ 0 atol=tol2 end @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64) diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index 2f9f0d763cbc3..91106676c52f6 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -174,7 +174,7 @@ end @testset "Issue 7304" begin A = [-√.5 -√.5; -√.5 √.5] Q = rectangularQ(qr(A).Q) - @test vecnorm(A-Q) < eps() + @test norm(A-Q) < eps() end @testset "qr on AbstractVector" begin diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 46aaeaf3140c0..13e2aeba68e7a 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -251,7 +251,7 @@ for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFlo if !(elty1 in (BigFloat, Complex{BigFloat})) # Not handled yet vals, vecs = eigen(A1) if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues. - @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(eps(float(real(one(vals[1])))))*(norm(A1,Inf)*n)^2 + @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(eps(float(real(one(vals[1])))))*(opnorm(A1,Inf)*n)^2 end end diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 03db4c90f6109..02d2ea36f16df 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -20,7 +20,7 @@ srand(123) @test -one(UniformScaling(2)) == UniformScaling(-1) @test sparse(3I,4,5) == sparse(1:4, 1:4, 3, 4, 5) @test sparse(3I,5,4) == sparse(1:4, 1:4, 3, 5, 4) - @test norm(UniformScaling(1+im)) ≈ sqrt(2) + @test opnorm(UniformScaling(1+im)) ≈ sqrt(2) end @testset "conjugation of UniformScaling" begin diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl index 33330420d385a..ece80669c4cc7 100644 --- a/stdlib/SparseArrays/src/SparseArrays.jl +++ b/stdlib/SparseArrays/src/SparseArrays.jl @@ -13,8 +13,8 @@ using LinearAlgebra import Base: +, -, *, \, /, &, |, xor, == import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, dot, eigen, - issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, vecdot, - vecnorm, cond, diagm, factorize, ishermitian, norm, lmul!, rmul!, tril, triu + issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, + cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu import Base: @get!, acos, acosd, acot, acotd, acsch, asech, asin, asind, asinh, atan, atand, atanh, broadcast!, conj!, cos, cosc, cosd, cosh, cospi, cot, diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index c70bc6ac12c17..8b81eacf42f54 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -560,15 +560,15 @@ end diff(a::SparseMatrixCSC; dims::Integer) = dims==1 ? sparse_diff1(a) : sparse_diff2(a) ## norm and rank -vecnorm(A::SparseMatrixCSC, p::Real=2) = vecnorm(view(A.nzval, 1:nnz(A)), p) +norm(A::SparseMatrixCSC, p::Real=2) = norm(view(A.nzval, 1:nnz(A)), p) -function norm(A::SparseMatrixCSC,p::Real=2) +function opnorm(A::SparseMatrixCSC, p::Real=2) m, n = size(A) if m == 0 || n == 0 || isempty(A) return float(real(zero(eltype(A)))) elseif m == 1 || n == 1 # TODO: compute more efficiently using A.nzval directly - return norm(Array(A), p) + return opnorm(Array(A), p) else Tnorm = typeof(float(real(zero(eltype(A))))) Tsum = promote_type(Float64,Tnorm) @@ -583,7 +583,7 @@ function norm(A::SparseMatrixCSC,p::Real=2) end return convert(Tnorm, nA) elseif p==2 - throw(ArgumentError("2-norm not yet implemented for sparse matrices. Try norm(Array(A)) or norm(A, p) where p=1 or Inf.")) + throw(ArgumentError("2-norm not yet implemented for sparse matrices. Try opnorm(Array(A)) or opnorm(A, p) where p=1 or Inf.")) elseif p==Inf rowSum = zeros(Tsum,m) for i=1:length(A.nzval) @@ -592,7 +592,7 @@ function norm(A::SparseMatrixCSC,p::Real=2) return convert(Tnorm, maximum(rowSum)) end end - throw(ArgumentError("invalid p-norm p=$p. Valid: 1, Inf")) + throw(ArgumentError("invalid operator p-norm p=$p. Valid: 1, Inf")) end # TODO rank @@ -600,12 +600,12 @@ end # cond function cond(A::SparseMatrixCSC, p::Real=2) if p == 1 - normAinv = normestinv(A) - normA = norm(A, 1) + normAinv = opnormestinv(A) + normA = opnorm(A, 1) return normA * normAinv elseif p == Inf - normAinv = normestinv(copy(A')) - normA = norm(A, Inf) + normAinv = opnormestinv(copy(A')) + normA = opnorm(A, Inf) return normA * normAinv elseif p == 2 throw(ArgumentError("2-norm condition number is not implemented for sparse matrices, try cond(Array(A), 2) instead")) @@ -614,7 +614,7 @@ function cond(A::SparseMatrixCSC, p::Real=2) end end -function normestinv(A::SparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T +function opnormestinv(A::SparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T maxiter = 5 # Check the input n = checksquare(A) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 313220b80b735..1fd076e737bd4 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -1320,7 +1320,7 @@ for f in [:sum, :maximum, :minimum], op in [:abs, :abs2] end end -vecnorm(x::SparseVectorUnion, p::Real=2) = vecnorm(nonzeros(x), p) +norm(x::SparseVectorUnion, p::Real=2) = norm(nonzeros(x), p) ### linalg.jl diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index 9684279b5873c..3e649bf227332 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -1568,8 +1568,8 @@ end @test norm(A) == zero(eltype(A)) A = sparse([1.0]) @test norm(A) == 1.0 - @test_throws ArgumentError norm(sprand(5,5,0.2),3) - @test_throws ArgumentError norm(sprand(5,5,0.2),2) + @test_throws ArgumentError opnorm(sprand(5,5,0.2),3) + @test_throws ArgumentError opnorm(sprand(5,5,0.2),2) end @testset "ishermitian/issymmetric" begin @@ -1691,30 +1691,30 @@ end Ac = sprandn(10,10,.1) + im* sprandn(10,10,.1) Ar = sprandn(10,10,.1) Ai = ceil.(Int,Ar*100) - @test norm(Ac,1) ≈ norm(Array(Ac),1) - @test norm(Ac,Inf) ≈ norm(Array(Ac),Inf) - @test vecnorm(Ac) ≈ vecnorm(Array(Ac)) - @test norm(Ar,1) ≈ norm(Array(Ar),1) - @test norm(Ar,Inf) ≈ norm(Array(Ar),Inf) - @test vecnorm(Ar) ≈ vecnorm(Array(Ar)) - @test norm(Ai,1) ≈ norm(Array(Ai),1) - @test norm(Ai,Inf) ≈ norm(Array(Ai),Inf) - @test vecnorm(Ai) ≈ vecnorm(Array(Ai)) + @test opnorm(Ac,1) ≈ opnorm(Array(Ac),1) + @test opnorm(Ac,Inf) ≈ opnorm(Array(Ac),Inf) + @test norm(Ac) ≈ norm(Array(Ac)) + @test opnorm(Ar,1) ≈ opnorm(Array(Ar),1) + @test opnorm(Ar,Inf) ≈ opnorm(Array(Ar),Inf) + @test norm(Ar) ≈ norm(Array(Ar)) + @test opnorm(Ai,1) ≈ opnorm(Array(Ai),1) + @test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf) + @test norm(Ai) ≈ norm(Array(Ai)) Ai = trunc.(Int, Ar*100) - @test norm(Ai,1) ≈ norm(Array(Ai),1) - @test norm(Ai,Inf) ≈ norm(Array(Ai),Inf) - @test vecnorm(Ai) ≈ vecnorm(Array(Ai)) + @test opnorm(Ai,1) ≈ opnorm(Array(Ai),1) + @test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf) + @test norm(Ai) ≈ norm(Array(Ai)) Ai = round.(Int, Ar*100) - @test norm(Ai,1) ≈ norm(Array(Ai),1) - @test norm(Ai,Inf) ≈ norm(Array(Ai),Inf) - @test vecnorm(Ai) ≈ vecnorm(Array(Ai)) + @test opnorm(Ai,1) ≈ opnorm(Array(Ai),1) + @test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf) + @test norm(Ai) ≈ norm(Array(Ai)) # make certain entries in nzval beyond # the range specified in colptr do not - # impact vecnorm of a sparse matrix + # impact norm of a sparse matrix foo = sparse(1.0I, 4, 4) resize!(foo.nzval, 5) setindex!(foo.nzval, NaN, 5) - @test vecnorm(foo) == 2.0 + @test norm(foo) == 2.0 end @testset "sparse matrix cond" begin @@ -1724,10 +1724,10 @@ end @test cond(A, 1) == 1.0 # For a discussion of the tolerance, see #14778 if Base.USE_GPL_LIBS - @test 0.99 <= cond(Ar, 1) \ norm(Ar, 1) * norm(inv(Array(Ar)), 1) < 3 - @test 0.99 <= cond(Ac, 1) \ norm(Ac, 1) * norm(inv(Array(Ac)), 1) < 3 - @test 0.99 <= cond(Ar, Inf) \ norm(Ar, Inf) * norm(inv(Array(Ar)), Inf) < 3 - @test 0.99 <= cond(Ac, Inf) \ norm(Ac, Inf) * norm(inv(Array(Ac)), Inf) < 3 + @test 0.99 <= cond(Ar, 1) \ opnorm(Ar, 1) * opnorm(inv(Array(Ar)), 1) < 3 + @test 0.99 <= cond(Ac, 1) \ opnorm(Ac, 1) * opnorm(inv(Array(Ac)), 1) < 3 + @test 0.99 <= cond(Ar, Inf) \ opnorm(Ar, Inf) * opnorm(inv(Array(Ar)), Inf) < 3 + @test 0.99 <= cond(Ac, Inf) \ opnorm(Ac, Inf) * opnorm(inv(Array(Ac)), Inf) < 3 end @test_throws ArgumentError cond(A,2) @test_throws ArgumentError cond(A,3) @@ -1737,21 +1737,21 @@ end @test_throws DimensionMismatch cond(Arect, Inf) end -@testset "sparse matrix normestinv" begin +@testset "sparse matrix opnormestinv" begin srand(1234) Ac = sprandn(20,20,.5) + im* sprandn(20,20,.5) Aci = ceil.(Int64, 100*sprand(20,20,.5)) + im*ceil.(Int64, sprand(20,20,.5)) Ar = sprandn(20,20,.5) Ari = ceil.(Int64, 100*Ar) if Base.USE_GPL_LIBS - # NOTE: normestinv is probabilistic, so requires a fixed seed (set above in srand(1234)) - @test SparseArrays.normestinv(Ac,3) ≈ norm(inv(Array(Ac)),1) atol=1e-4 - @test SparseArrays.normestinv(Aci,3) ≈ norm(inv(Array(Aci)),1) atol=1e-4 - @test SparseArrays.normestinv(Ar) ≈ norm(inv(Array(Ar)),1) atol=1e-4 - @test_throws ArgumentError SparseArrays.normestinv(Ac,0) - @test_throws ArgumentError SparseArrays.normestinv(Ac,21) - end - @test_throws DimensionMismatch SparseArrays.normestinv(sprand(3,5,.9)) + # NOTE: opnormestinv is probabilistic, so requires a fixed seed (set above in srand(1234)) + @test SparseArrays.opnormestinv(Ac,3) ≈ opnorm(inv(Array(Ac)),1) atol=1e-4 + @test SparseArrays.opnormestinv(Aci,3) ≈ opnorm(inv(Array(Aci)),1) atol=1e-4 + @test SparseArrays.opnormestinv(Ar) ≈ opnorm(inv(Array(Ar)),1) atol=1e-4 + @test_throws ArgumentError SparseArrays.opnormestinv(Ac,0) + @test_throws ArgumentError SparseArrays.opnormestinv(Ac,21) + end + @test_throws DimensionMismatch SparseArrays.opnormestinv(sprand(3,5,.9)) end @testset "issue #13008" begin diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index cf37d09489b88..e2f047a3eb0ee 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -708,16 +708,16 @@ end ### Reduction -@testset "sum, vecnorm" begin +@testset "sum, norm" begin x = spv_x1 @test sum(x) == 4.0 @test sum(abs, x) == 5.5 @test sum(abs2, x) == 14.375 - @test vecnorm(x) == sqrt(14.375) - @test vecnorm(x, 1) == 5.5 - @test vecnorm(x, 2) == sqrt(14.375) - @test vecnorm(x, Inf) == 3.5 + @test norm(x) == sqrt(14.375) + @test norm(x, 1) == 5.5 + @test norm(x, 2) == sqrt(14.375) + @test norm(x, Inf) == 3.5 end @testset "maximum, minimum" begin diff --git a/stdlib/SuiteSparse/test/cholmod.jl b/stdlib/SuiteSparse/test/cholmod.jl index 57a250bd4e8aa..bf255e145b0ab 100644 --- a/stdlib/SuiteSparse/test/cholmod.jl +++ b/stdlib/SuiteSparse/test/cholmod.jl @@ -290,8 +290,8 @@ end @test_throws BoundsError ADense[6, 1] @test_throws BoundsError ADense[1, 6] @test copy(ADense) == ADense - @test CHOLMOD.norm_dense(ADense, 1) ≈ norm(A, 1) - @test CHOLMOD.norm_dense(ADense, 0) ≈ norm(A, Inf) + @test CHOLMOD.norm_dense(ADense, 1) ≈ opnorm(A, 1) + @test CHOLMOD.norm_dense(ADense, 0) ≈ opnorm(A, Inf) @test_throws ArgumentError CHOLMOD.norm_dense(ADense, 2) @test_throws ArgumentError CHOLMOD.norm_dense(ADense, 3) From e90319a82f473ecbde8a18ad777daca79c33f943 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 3 Jun 2018 17:32:42 +0200 Subject: [PATCH 02/17] fixed vecnorm -> norm in tests for IterativeEigenolvers, complex, offsetarray --- test/complex.jl | 2 +- test/offsetarray.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/complex.jl b/test/complex.jl index 6dde71594ff53..67be8f2a5fa4d 100644 --- a/test/complex.jl +++ b/test/complex.jl @@ -929,7 +929,7 @@ end @inferred sin(x) @inferred cos(x) @inferred norm(x) - @inferred vecnorm(x) + @inferred opnorm(x) end end diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 3a97c7feda760..a566c90e46e22 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -376,8 +376,8 @@ I = findall(!iszero, z) @test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0)) @test sum(OffsetArray(fill(1,3000), -1000)) == 3000 -@test vecnorm(v) ≈ vecnorm(parent(v)) -@test vecnorm(A) ≈ vecnorm(parent(A)) +@test norm(v) ≈ norm(parent(v)) +@test norm(A) ≈ norm(parent(A)) @test vecdot(v, v) ≈ vecdot(v0, v0) # Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of From 7813808c07dfd3ec97fe058d8216538d9bd21157 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 3 Jun 2018 17:51:28 +0200 Subject: [PATCH 03/17] dot,vecdot -> inner; deprecate vecdot --- stdlib/LinearAlgebra/docs/src/index.md | 2 +- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 8 ++-- stdlib/LinearAlgebra/src/adjtrans.jl | 4 +- stdlib/LinearAlgebra/src/dense.jl | 2 +- stdlib/LinearAlgebra/src/deprecated.jl | 9 ++-- stdlib/LinearAlgebra/src/generic.jl | 53 +++++++++++------------ stdlib/LinearAlgebra/src/hessenberg.jl | 2 +- stdlib/LinearAlgebra/src/lapack.jl | 4 +- stdlib/LinearAlgebra/src/matmul.jl | 10 ++--- stdlib/LinearAlgebra/src/qr.jl | 2 +- stdlib/LinearAlgebra/test/adjtrans.jl | 8 ++-- stdlib/LinearAlgebra/test/blas.jl | 2 +- stdlib/LinearAlgebra/test/matmul.jl | 37 ++++++++-------- stdlib/SparseArrays/src/SparseArrays.jl | 2 +- stdlib/SparseArrays/src/sparsevector.jl | 24 +++++----- stdlib/SparseArrays/test/sparsevector.jl | 32 +++++++------- stdlib/SuiteSparse/src/spqr.jl | 4 +- 17 files changed, 104 insertions(+), 101 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index f8acc46722bef..dbabc6a98558d 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -298,8 +298,8 @@ Linear algebra functions in Julia are largely implemented by calling functions f ```@docs Base.:*(::AbstractMatrix, ::AbstractMatrix) Base.:\(::AbstractMatrix, ::AbstractVecOrMat) +LinearAlgebra.inner LinearAlgebra.dot -LinearAlgebra.vecdot LinearAlgebra.cross LinearAlgebra.factorize LinearAlgebra.Diagonal diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 6d4e082914350..0443de25e263e 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -79,7 +79,6 @@ export diag, diagind, diagm, - dot, eigen, eigen!, eigmax, @@ -91,6 +90,7 @@ export givens, hessenberg, hessenberg!, + inner, isdiag, ishermitian, isposdef, @@ -144,7 +144,6 @@ export triu, tril!, triu!, - vecdot, # Operators \, @@ -371,7 +370,10 @@ include("schur.jl") include("structuredbroadcast.jl") include("deprecated.jl") -const ⋅ = dot +#TODO: This yields: error during bootstrap. LoadError("sysimg.jl", 213, ...) +# const dot = inner +# export dot +const ⋅ = inner const × = cross export ⋅, × diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index d4b08a07502a6..70fed4c9439ae 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -194,8 +194,8 @@ broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs... ## multiplication * # Adjoint/Transpose-vector * vector -*(u::AdjointAbsVec, v::AbstractVector) = dot(u.parent, v) -*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v) +*(u::AdjointAbsVec, v::AbstractVector) = inner(u.parent, v) +*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = inner(u.parent, v) function *(u::TransposeAbsVec, v::AbstractVector) @boundscheck length(u) == length(v) || throw(DimensionMismatch()) return sum(@inbounds(u[k]*v[k]) for k in 1:length(u)) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 9c121abe26514..d7328c8c7f7d5 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -5,7 +5,7 @@ ## BLAS cutoff threshold constants const SCAL_CUTOFF = 2048 -const DOT_CUTOFF = 128 +const INNER_CUTOFF = 128 const ASUM_CUTOFF = 32 const NRM2_CUTOFF = 32 diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 4ed5e5f8a2734..3facd40ca9b7a 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -9,6 +9,9 @@ using Base: @deprecate, depwarn # deprecate vecnorm in favor of norm @deprecate vecnorm norm +# deprecate vecdot in favor of inner +@deprecate vecdot inner + # PR #22188 export cholfact, cholfact! @deprecate cholfact!(A::StridedMatrix, uplo::Symbol, ::Type{Val{false}}) cholesky!(Hermitian(A, uplo), Val(false)) @@ -541,9 +544,9 @@ IndexStyle(::Type{<:RowVector}) = IndexLinear() # Multiplication # # inner product -> dot product specializations -@inline *(rowvec::RowVector{T}, vec::AbstractVector{T}) where {T<:Real} = dot(parent(rowvec), vec) -@inline *(rowvec::ConjRowVector{T}, vec::AbstractVector{T}) where {T<:Real} = dot(rvadjoint(rowvec), vec) -@inline *(rowvec::ConjRowVector, vec::AbstractVector) = dot(rvadjoint(rowvec), vec) +@inline *(rowvec::RowVector{T}, vec::AbstractVector{T}) where {T<:Real} = inner(parent(rowvec), vec) +@inline *(rowvec::ConjRowVector{T}, vec::AbstractVector{T}) where {T<:Real} = inner(rvadjoint(rowvec), vec) +@inline *(rowvec::ConjRowVector, vec::AbstractVector) = inner(rvadjoint(rowvec), vec) # Generic behavior @inline function *(rowvec::RowVector, vec::AbstractVector) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 7f76a32ed9016..41be3f98e7495 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -615,56 +615,56 @@ opnorm(v::TransposeAbsVec) = norm(v.parent) norm(v::Union{TransposeAbsVec,AdjointAbsVec}, p::Real) = norm(v.parent, p) -function vecdot(x::AbstractArray, y::AbstractArray) +function inner(x::AbstractArray, y::AbstractArray) lx = _length(x) if lx != _length(y) 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 inner(zero(eltype(x)), zero(eltype(y))) end - s = zero(dot(first(x), first(y))) + s = zero(inner(first(x), first(y))) for (Ix, Iy) in zip(eachindex(x), eachindex(y)) - @inbounds s += dot(x[Ix], y[Iy]) + @inbounds s += inner(x[Ix], y[Iy]) end s end """ - vecdot(x, y) + inner(x, y) For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or -any element type for which `dot` is defined), compute the Euclidean dot product (the sum of -`dot(x[i],y[i])`) as if they were vectors. +any element type for which `inner` is defined), compute the inner product (or dot product +or scalar product, i.e. the sum of `inner(x[i],y[i])`) as if they were vectors. # Examples ```jldoctest -julia> vecdot(1:5, 2:6) +julia> inner(1:5, 2:6) 70 julia> x = fill(2., (5,5)); julia> y = fill(3., (5,5)); -julia> vecdot(x, y) +julia> inner(x, y) 150.0 ``` """ -function vecdot(x, y) # arbitrary iterables +function inner(x, y) # arbitrary iterables ix = iterate(x) iy = iterate(y) if ix === nothing if iy !== nothing throw(DimensionMismatch("x and y are of different lengths!")) end - return dot(zero(eltype(x)), zero(eltype(y))) + return inner(zero(eltype(x)), zero(eltype(y))) end if iy === nothing throw(DimensionMismatch("x and y are of different lengths!")) end (vx, xs) = ix (vy, ys) = iy - s = dot(vx, vy) + s = inner(vx, vy) while true ix = iterate(x, xs) iy = iterate(y, ys) @@ -672,7 +672,7 @@ function vecdot(x, y) # arbitrary iterables break end (vx, xs), (vy, ys) = ix, iy - s += dot(vx, vy) + s += inner(vx, vy) end if !(iy == nothing && ix == nothing) throw(DimensionMismatch("x and y are of different lengths!")) @@ -680,49 +680,46 @@ function vecdot(x, y) # arbitrary iterables return s end -vecdot(x::Number, y::Number) = conj(x) * y - -dot(x::Number, y::Number) = vecdot(x, y) +inner(x::Number, y::Number) = conj(x) * y """ + inner(x,y) dot(x, y) ⋅(x,y) -Compute the dot product between two vectors. For complex vectors, the first vector is conjugated. -When the vectors have equal lengths, calling `dot` is semantically equivalent to `sum(vx'vy for (vx,vy) in zip(x, y))`. +Compute the inner/dot product between two vectors. For complex vectors, the first +vector is conjugated. When the vectors have equal lengths, calling `inner` is +semantically equivalent to `sum(inner(vx,vy) for (vx,vy) in zip(x, y))`. # Examples ```jldoctest -julia> dot([1; 1], [2; 3]) +julia> inner([1; 1], [2; 3]) 5 -julia> dot([im; im], [1; 1]) +julia> inner([im; im], [1; 1]) 0 - 2im ``` """ -function dot(x::AbstractVector, y::AbstractVector) +function inner(x::AbstractVector, y::AbstractVector) if length(x) != length(y) - throw(DimensionMismatch("dot product arguments have lengths $(length(x)) and $(length(y))")) + throw(DimensionMismatch("inner product arguments have unequal lengths $(length(x)) and $(length(y))")) end ix = iterate(x) if ix === nothing # we only need to check the first vector, since equal lengths have been asserted - return zero(eltype(x))'zero(eltype(y)) + return inner(zero(eltype(x)), zero(eltype(y))) end iy = iterate(y) - s = ix[1]'iy[1] + s = inner(ix[1], iy[1]) ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) while ix != nothing - s += ix[1]'iy[1] + s += inner(ix[1], iy[1]) ix = iterate(x, ix[2]) iy = iterate(y, iy[2]) end return s end -# Call optimized BLAS methods for vectors of numbers -dot(x::AbstractVector{<:Number}, y::AbstractVector{<:Number}) = vecdot(x, y) - ########################################################################################### diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 050462c19d55f..19302710180cd 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -84,7 +84,7 @@ function getindex(A::HessenbergQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return dot(x, lmul!(A, y)) + return inner(x, lmul!(A, y)) end ## reconstruct the original matrix diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index 226c7aca5d222..4084fe116607d 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -12,7 +12,7 @@ import ..LinearAlgebra.BLAS.@blasfunc import ..LinearAlgebra: BlasFloat, BlasInt, LAPACKException, DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare -using ..LinearAlgebra: triu, tril, dot +using ..LinearAlgebra: triu, tril, inner using Base: iszero @@ -1444,7 +1444,7 @@ for (gglse, elty) in ((:dgglse_, :Float64), resize!(work, lwork) end end - X, dot(view(c, n - p + 1:m), view(c, n - p + 1:m)) + X, inner(view(c, n - p + 1:m), view(c, n - p + 1:m)) end end end diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 7b6e4a54437ea..0f787cafd7196 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -4,12 +4,12 @@ matprod(x, y) = x*y + x*y -# Dot products +# inner/dot products -vecdot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y) -vecdot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y) +inner(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y) +inner(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y) -function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasReal,TI<:Integer} +function inner(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasReal,TI<:Integer} if length(rx) != length(ry) throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))")) end @@ -22,7 +22,7 @@ function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector GC.@preserve x y BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry)) end -function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasComplex,TI<:Integer} +function inner(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasComplex,TI<:Integer} if length(rx) != length(ry) throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))")) end diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 602c143230cfb..97d9151428abc 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -469,7 +469,7 @@ function getindex(A::AbstractQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return dot(x, lmul!(A, y)) + return inner(x, lmul!(A, y)) end ## Multiplication by Q diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index ae140d5f7a8c8..35596e994a6fb 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -342,10 +342,10 @@ end realvec, realmat = [1, 2, 3], [1 2 3; 4 5 6; 7 8 9] complexvec, complexmat = [1im, 2, -3im], [1im 2 3; 4 5 -6im; 7im 8 9] # Adjoint/Transpose-vector * vector - @test Adjoint(realvec) * realvec == dot(realvec, realvec) - @test Transpose(realvec) * realvec == dot(realvec, realvec) - @test Adjoint(complexvec) * complexvec == dot(complexvec, complexvec) - @test Transpose(complexvec) * complexvec == dot(conj(complexvec), complexvec) + @test Adjoint(realvec) * realvec == inner(realvec, realvec) + @test Transpose(realvec) * realvec == inner(realvec, realvec) + @test Adjoint(complexvec) * complexvec == inner(complexvec, complexvec) + @test Transpose(complexvec) * complexvec == inner(conj(complexvec), complexvec) # vector * Adjoint/Transpose-vector @test realvec * Adjoint(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3))) @test realvec * Transpose(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3))) diff --git a/stdlib/LinearAlgebra/test/blas.jl b/stdlib/LinearAlgebra/test/blas.jl index d09f35899c358..bf5a578cd4f54 100644 --- a/stdlib/LinearAlgebra/test/blas.jl +++ b/stdlib/LinearAlgebra/test/blas.jl @@ -53,7 +53,7 @@ srand(100) v41 = convert(Vector{elty}, [4:-1:1;]) let n = 10 - @testset "dot products" begin + @testset "inner/dot products" begin if elty <: Real x1 = convert(Vector{elty}, randn(n)) x2 = convert(Vector{elty}, randn(n)) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 5146e29d72c9a..ec06f0964a2f6 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -215,32 +215,33 @@ end end # issue #6450 -@test dot(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 +@test inner(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 -@testset "dot" for elty in (Float32, Float64, ComplexF32, ComplexF64) +@testset "inner" for elty in (Float32, Float64, ComplexF32, ComplexF64) x = convert(Vector{elty},[1.0, 2.0, 3.0]) y = convert(Vector{elty},[3.5, 4.5, 5.5]) - @test_throws DimensionMismatch dot(x, 1:2, y, 1:3) - @test_throws BoundsError dot(x, 1:4, y, 1:4) - @test_throws BoundsError dot(x, 1:3, y, 2:4) - @test dot(x, 1:2,y, 1:2) == convert(elty, 12.5) + @test_throws DimensionMismatch inner(x, 1:2, y, 1:3) + @test_throws BoundsError inner(x, 1:4, y, 1:4) + @test_throws BoundsError inner(x, 1:3, y, 2:4) + @test inner(x, 1:2, y, 1:2) == convert(elty, 12.5) @test transpose(x)*y == convert(elty, 29.0) - @test_throws MethodError dot(rand(elty, 2, 2), randn(elty, 2, 2)) - X = convert(Vector{Matrix{elty}},[reshape(1:4, 2, 2), fill(1, 2, 2)]) - res = convert(Matrix{elty}, [7.0 13.0; 13.0 27.0]) - @test dot(X, X) == res + X = convert(Matrix{elty},[1.0 2.0; 3.0 4.0]) + Y = convert(Matrix{elty},[1.5 2.5; 3.5 4.5]) + @test inner(X, Y) == convert(elty, 35.0) + Z = convert(Vector{Matrix{elty}},[reshape(1:4, 2, 2), fill(1, 2, 2)]) + @test inner(Z, Z) == convert(elty, 34.0) end -vecdot_(x,y) = invoke(vecdot, Tuple{Any,Any}, x,y) -@testset "generic vecdot" begin +inner_(x,y) = invoke(inner, Tuple{Any,Any}, x,y) +@testset "generic inner" begin AA = [1+2im 3+4im; 5+6im 7+8im] BB = [2+7im 4+1im; 3+8im 6+5im] for A in (copy(AA), view(AA, 1:2, 1:2)), B in (copy(BB), view(BB, 1:2, 1:2)) - @test vecdot(A,B) == dot(vec(A),vec(B)) == vecdot_(A,B) == vecdot(float.(A),float.(B)) - @test vecdot(Int[], Int[]) == 0 == vecdot_(Int[], Int[]) - @test_throws MethodError vecdot(Any[], Any[]) - @test_throws MethodError vecdot_(Any[], Any[]) - for n1 = 0:2, n2 = 0:2, d in (vecdot, vecdot_) + @test inner(A,B) == inner(vec(A),vec(B)) == inner_(A,B) == inner(float.(A),float.(B)) + @test inner(Int[], Int[]) == 0 == inner_(Int[], Int[]) + @test_throws MethodError inner(Any[], Any[]) + @test_throws MethodError inner_(Any[], Any[]) + for n1 = 0:2, n2 = 0:2, d in (inner, inner_) if n1 != n2 @test_throws DimensionMismatch d(1:n1, 1:n2) else @@ -424,7 +425,7 @@ end @test transpose(Xv2)*Xv2 ≈ XtX @test (transpose(Xv3)*Xv3)[1] ≈ XtX @test Xv1'*Xv1 ≈ XcX - @test Xv2'*Xv2 ≈ XcX + @test_throws MethodError Xv2'*Xv2 ≈ XcX @test (Xv3'*Xv3)[1] ≈ XcX @test (Xv1*transpose(Xv1))[1] ≈ XXt @test Xv2*transpose(Xv2) ≈ XXt diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl index ece80669c4cc7..51a62a0417408 100644 --- a/stdlib/SparseArrays/src/SparseArrays.jl +++ b/stdlib/SparseArrays/src/SparseArrays.jl @@ -12,7 +12,7 @@ using Base.Sort: Forward using LinearAlgebra import Base: +, -, *, \, /, &, |, xor, == -import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, dot, eigen, +import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, eigen, inner, issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 1fd076e737bd4..4a2ac4bb99637 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -1385,27 +1385,27 @@ end (*)(a::Number, x::SparseVectorUnion) = SparseVector(length(x), copy(nonzeroinds(x)), a * nonzeros(x)) (/)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) / a) -# dot -function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} +# inner/dot +function inner(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} n = length(x) length(y) == n || throw(DimensionMismatch()) nzind = nonzeroinds(y) nzval = nonzeros(y) - s = zero(Tx) * zero(Ty) + s = inner(zero(Tx), zero(Ty)) for i = 1:length(nzind) - s += conj(x[nzind[i]]) * nzval[i] + s += inner(x[nzind[i]], nzval[i]) end return s end -function dot(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} +function inner(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} n = length(y) length(x) == n || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) - s = zero(Tx) * zero(Ty) + s = inner(zero(Tx), zero(Ty)) @inbounds for i = 1:length(nzind) - s += conj(nzval[i]) * y[nzind[i]] + s += inner(nzval[i], y[nzind[i]]) end return s end @@ -1431,7 +1431,7 @@ function _spdot(f::Function, s end -function dot(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) +function inner(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) x === y && return sum(abs2, x) n = length(x) length(y) == n || throw(DimensionMismatch()) @@ -1441,7 +1441,7 @@ function dot(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) xnzval = nonzeros(x) ynzval = nonzeros(y) - _spdot(dot, + _spdot(inner, 1, length(xnzind), xnzind, xnzval, 1, length(ynzind), ynzind, ynzval) end @@ -1612,7 +1612,7 @@ mul!(y::AbstractVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractS (A = adjA.parent; mul!(y, adjoint(A), x, one(Tx), zero(Ty))) mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, α::Number, β::Number) = - (A = adjA.parent; _At_or_Ac_mul_B!(dot, y, A, x, α, β)) + (A = adjA.parent; _At_or_Ac_mul_B!(inner, y, A, x, α, β)) function _At_or_Ac_mul_B!(tfun::Function, y::AbstractVector, A::SparseMatrixCSC, x::AbstractSparseVector, @@ -1633,7 +1633,7 @@ function _At_or_Ac_mul_B!(tfun::Function, mx = length(xnzind) for j = 1:n - # s <- dot(A[:,j], x) + # s <- inner(A[:,j], x) s = _spdot(tfun, Acolptr[j], Acolptr[j+1]-1, Arowval, Anzval, 1, mx, xnzind, xnzval) @inbounds y[j] += s * α @@ -1654,7 +1654,7 @@ end (A = transA.parent; _At_or_Ac_mul_B(*, A, x)) *(adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector) = - (A = adjA.parent; _At_or_Ac_mul_B(dot, A, x)) + (A = adjA.parent; _At_or_Ac_mul_B(inner, A, x)) function _At_or_Ac_mul_B(tfun::Function, A::SparseMatrixCSC{TvA,TiA}, x::AbstractSparseVector{TvX,TiX}) where {TvA,TiA,TvX,TiX} m, n = size(A) diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index e2f047a3eb0ee..58dfd24dfb5c0 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -795,14 +795,14 @@ end @test exact_equal(xc, sx) end - @testset "dot" begin - dv = dot(xf, xf2) - @test dot(x, x) == sum(abs2, x) - @test dot(x2, x2) == sum(abs2, x2) - @test dot(x, x2) ≈ dv - @test dot(x2, x) ≈ dv - @test dot(Array(x), x2) ≈ dv - @test dot(x, Array(x2)) ≈ dv + @testset "inner/dot" begin + dv = inner(xf, xf2) + @test inner(x, x) == sum(abs2, x) + @test inner(x2, x2) == sum(abs2, x2) + @test inner(x, x2) ≈ dv + @test inner(x2, x) ≈ dv + @test inner(Array(x), x2) ≈ dv + @test inner(x, Array(x2)) ≈ dv end end @@ -810,8 +810,8 @@ end y = complex.(sprand(32, 0.6), sprand(32, 0.6)) xf = Array(x)::Vector{ComplexF64} yf = Array(y)::Vector{ComplexF64} - @test dot(x, x) ≈ dot(xf, xf) - @test dot(x, y) ≈ dot(xf, yf) + @test inner(x, x) ≈ inner(xf, xf) + @test inner(x, y) ≈ inner(xf, yf) end end @@ -1247,12 +1247,12 @@ end A = sprandn(n, n, 0.01) for j in 1:5:n Aj, Ajview = A[:, j], view(A, :, j) - @test norm(Aj) == norm(Ajview) - @test dot(Aj, copy(Aj)) == dot(Ajview, Aj) # don't alias since it takes a different code path - @test rmul!(Aj, 0.1) == rmul!(Ajview, 0.1) - @test Aj*0.1 == Ajview*0.1 - @test 0.1*Aj == 0.1*Ajview - @test Aj/0.1 == Ajview/0.1 + @test norm(Aj) == norm(Ajview) + @test inner(Aj, copy(Aj)) == inner(Ajview, Aj) # don't alias since it takes a different code path + @test rmul!(Aj, 0.1) == rmul!(Ajview, 0.1) + @test Aj*0.1 == Ajview*0.1 + @test 0.1*Aj == 0.1*Ajview + @test Aj/0.1 == Ajview/0.1 @test LinearAlgebra.axpy!(1.0, Aj, sparse(fill(1., n))) == LinearAlgebra.axpy!(1.0, Ajview, sparse(fill(1., n))) @test LinearAlgebra.lowrankupdate!(Matrix(1.0*I, n, n), fill(1.0, n), Aj) == diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl index 5059f1d75a587..b66faa07c3b58 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -206,7 +206,7 @@ function LinearAlgebra.lmul!(Q::QRSparseQ, A::StridedVecOrMat) h = view(Q.factors, :, l) for j in 1:size(A, 2) a = view(A, :, j) - LinearAlgebra.axpy!(τl*dot(h, a), h, a) + LinearAlgebra.axpy!(τl*inner(h, a), h, a) end end return A @@ -236,7 +236,7 @@ function LinearAlgebra.lmul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMa h = view(Q.factors, :, l) for j in 1:size(A, 2) a = view(A, :, j) - LinearAlgebra.axpy!(τl'*dot(h, a), h, a) + LinearAlgebra.axpy!(τl'*inner(h, a), h, a) end end return A From e238d14b265df56158239f8d37850e61fe0ad5c9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 3 Jun 2018 20:11:35 +0200 Subject: [PATCH 04/17] add dot(x,y)=inner(x,y) --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 4 +--- stdlib/LinearAlgebra/src/generic.jl | 9 +++++++++ 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 0443de25e263e..2f27a038f13e3 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -79,6 +79,7 @@ export diag, diagind, diagm, + dot, eigen, eigen!, eigmax, @@ -370,9 +371,6 @@ include("schur.jl") include("structuredbroadcast.jl") include("deprecated.jl") -#TODO: This yields: error during bootstrap. LoadError("sysimg.jl", 213, ...) -# const dot = inner -# export dot const ⋅ = inner const × = cross export ⋅, × diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 41be3f98e7495..163c466e60941 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -720,6 +720,15 @@ function inner(x::AbstractVector, y::AbstractVector) return s end +""" + inner(x,y) + dot(x, y) + ⋅(x,y) + +Compute the inner/dot product between `x` and `y`, see [`inner`](@ref). +""" +dot(x, y) = inner(x, y) + ########################################################################################### From 5ff39c3e8f5342c0aa917b6bc98c7ab0d9867f8f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 4 Jun 2018 05:38:33 +0200 Subject: [PATCH 05/17] use dot/inner instead of vecdot in test/offsetarray.jl --- test/offsetarray.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index a566c90e46e22..22b7c16ae3bba 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -378,7 +378,7 @@ I = findall(!iszero, z) @test norm(v) ≈ norm(parent(v)) @test norm(A) ≈ norm(parent(A)) -@test vecdot(v, v) ≈ vecdot(v0, v0) +@test inner(v, v) ≈ inner(v0, v0) # Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of # accuracy in the tests, we need to use BigFloats with enlarged precision. From 1b4d1fa5d7baf77eea626e8e012a7c2e5c994b95 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 4 Jun 2018 20:18:34 +0200 Subject: [PATCH 06/17] replace length with length(LinearIndices(... --- stdlib/LinearAlgebra/src/generic.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 163c466e60941..245d8d5ac36e0 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -701,8 +701,8 @@ julia> inner([im; im], [1; 1]) ``` """ function inner(x::AbstractVector, y::AbstractVector) - if length(x) != length(y) - throw(DimensionMismatch("inner product arguments have unequal lengths $(length(x)) and $(length(y))")) + if length(LinearIndices(x)) != length(LinearIndices(y)) + throw(DimensionMismatch("inner product arguments have unequal lengths $(length(LinearIndices(x))) and $(length(LinearIndices(y)))")) end ix = iterate(x) if ix === nothing @@ -726,6 +726,15 @@ end ⋅(x,y) Compute the inner/dot product between `x` and `y`, see [`inner`](@ref). + +# Examples +```jldoctest +julia> dot([1; 1], [2; 3]) +5 + +julia> dot([im; im], [1; 1]) +0 - 2im +``` """ dot(x, y) = inner(x, y) From 1258f0a21ce54cc950e5ac016400d3df721e4ff6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Jun 2018 05:53:11 +0200 Subject: [PATCH 07/17] const dot = inner --- stdlib/LinearAlgebra/src/generic.jl | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 245d8d5ac36e0..54a8ea388355d 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -720,23 +720,7 @@ function inner(x::AbstractVector, y::AbstractVector) return s end -""" - inner(x,y) - dot(x, y) - ⋅(x,y) - -Compute the inner/dot product between `x` and `y`, see [`inner`](@ref). - -# Examples -```jldoctest -julia> dot([1; 1], [2; 3]) -5 - -julia> dot([im; im], [1; 1]) -0 - 2im -``` -""" -dot(x, y) = inner(x, y) +const dot = inner ########################################################################################### From 682bbc7d109f9336f2894140cca6b2350692b9d0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Jun 2018 12:45:57 +0200 Subject: [PATCH 08/17] norm(x,p) uses norm(xi) instead of norm(xi,p) --- stdlib/LinearAlgebra/src/generic.jl | 33 ++++++++++++++-------------- stdlib/LinearAlgebra/test/generic.jl | 4 ++-- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 54a8ea388355d..55f3038116341 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -252,12 +252,12 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons # special cases of norm; note that they don't need to handle isempty(x) function generic_normMinusInf(x) (v, s) = iterate(x)::Tuple - minabs = norm(v, -Inf) + minabs = norm(v) while true y = iterate(x, s) y === nothing && break (v, s) = y - vnorm = norm(v, -Inf) + vnorm = norm(v) minabs = ifelse(isnan(minabs) | (minabs < vnorm), minabs, vnorm) end return float(minabs) @@ -265,12 +265,12 @@ end function generic_normInf(x) (v, s) = iterate(x)::Tuple - maxabs = norm(v, Inf) + maxabs = norm(v) while true y = iterate(x, s) y === nothing && break (v, s) = y - vnorm = norm(v, Inf) + vnorm = norm(v) maxabs = ifelse(isnan(maxabs) | (maxabs > vnorm), maxabs, vnorm) end return float(maxabs) @@ -278,14 +278,14 @@ end function generic_norm1(x) (v, s) = iterate(x)::Tuple - av = float(norm(v, 1)) + av = float(norm(v)) T = typeof(av) sum::promote_type(Float64, T) = av while true y = iterate(x, s) y === nothing && break (v, s) = y - sum += norm(v, 1) + sum += norm(v) end return convert(T, sum) end @@ -330,25 +330,25 @@ function generic_normp(x, p) (maxabs == 0 || isinf(maxabs)) && return maxabs T = typeof(maxabs) else - T = typeof(float(norm(v,p))) + 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,p)^spp + sum::promote_type(Float64, T) = norm(v)^spp while true y = iterate(x, s) y === nothing && break (v, s) = y - sum += norm(v,p)^spp + sum += norm(v)^spp end return convert(T, sum^inv(spp)) else # rescaling - sum = (norm(v,p)/maxabs)^spp + sum = (norm(v)/maxabs)^spp while true y = iterate(x, s) y == nothing && break (v, s) = y - sum += (norm(v,p)/maxabs)^spp + sum += (norm(v)/maxabs)^spp end return convert(T, maxabs*sum^inv(spp)) end @@ -368,11 +368,12 @@ For any iterable container `A` (including arrays of any dimension) of numbers (o element type for which `norm` is defined), compute the `p`-norm (defaulting to `p=2`) as if `A` were a vector of the corresponding length. -The `p`-norm is defined as: +The `p`-norm is defined as ```math \\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p} ``` -with ``a_i`` the entries of ``A`` and ``n`` its length. +with ``a_i`` the entries of ``A``, ``| a_i |`` are their [`norm`](@ref)s, and +``n`` the length of ``A``. `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 @@ -407,13 +408,13 @@ julia> norm([1 2 3 4 5 6 7 8 9]) julia> norm([1 2 3 4 5 6 7 8 9]') 16.881943016134134 -julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) == norm([v,v], 1) +julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) != norm([v,v], 1) true julia> norm(hcat(v,v), 2) == norm(vcat(v,v), 2) == norm([v,v], 2) true -julia> norm(hcat(v,v), Inf) == norm(vcat(v,v), Inf) == norm([v,v], Inf) +julia> norm(hcat(v,v), Inf) == norm(vcat(v,v), Inf) != norm([v,v], Inf) true ``` """ @@ -430,7 +431,7 @@ function norm(itr, p::Real=2) elseif p == -Inf return normMinusInf(itr) else - normp(itr,p) + normp(itr, p) end end diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index 58a85357054bc..6675e025a7912 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -203,8 +203,8 @@ end x = Vector{Int}[[1,2], [3,4]] @test @inferred(norm(x)) ≈ sqrt(30) @test norm(x, 0) == length(x) - @test norm(x, 1) ≈ 10 - @test norm(x, 3) ≈ cbrt(100) + @test norm(x, 1) ≈ 5+sqrt(5) + @test norm(x, 3) ≈ cbrt(5^3 +sqrt(5)^3) end @testset "LinearAlgebra.axp(b)y! for element type without commutative multiplication" begin From 614776b6c3be2b9922380757dfdcc405cf6f8e2d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 7 Jun 2018 05:50:52 +0200 Subject: [PATCH 09/17] update docstrings of norm and inner --- stdlib/LinearAlgebra/src/generic.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 55f3038116341..49e6710cad2d2 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -373,13 +373,18 @@ The `p`-norm is defined as \\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p} ``` with ``a_i`` the entries of ``A``, ``| a_i |`` are their [`norm`](@ref)s, and -``n`` the length of ``A``. +``n`` the length of ``A``. Since the `p`-norm is computed using the [`norm`](@ref)s +of the entries of `A`, the `p`-norm of a vector of vectors is not compatible with +the interpretation of it as a block vector in general if `p != 2`. `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. If `A` is a matrix and `p=2`, then this is equivalent to the Frobenius norm. +The second argument `p` is not necessarily a part of the interface for `norm`, i.e. a custom +type may only implement `norm(A)` without second argument. + Use [`opnorm`](@ref) to compute the operator norm of a matrix. # Examples @@ -636,7 +641,10 @@ end For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or any element type for which `inner` is defined), compute the inner product (or dot product -or scalar product, i.e. the sum of `inner(x[i],y[i])`) as if they were vectors. +or scalar product), i.e. the sum of `inner(x[i],y[i])`, as if they were vectors. + +`dot(x, y)`` and `x ⋅ y` (where `⋅` can be typed by tab-completing `\cdot` in the REPL) are +synonyms for `inner(x, y)`. # Examples ```jldoctest From 83098b1175a2cbf413516f2f207ecb089e1bfda4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 9 Jun 2018 06:10:29 +0200 Subject: [PATCH 10/17] fixed isapprox for UniformScaling --- stdlib/LinearAlgebra/src/uniformscaling.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 01cba86974875..67d2456a5c316 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -196,9 +196,9 @@ function isapprox(J::UniformScaling, A::AbstractMatrix; rtol::Real = Base.rtoldefault(promote_leaf_eltypes(A), eltype(J), atol), nans::Bool = false, norm::Function = norm) n = checksquare(A) - normJ = norm === opnorm ? abs(J.λ) : - norm === norm ? abs(J.λ) * sqrt(n) : - norm(Diagonal(fill(J.λ, n))) + normJ = norm === opnorm ? abs(J.λ) : + norm === LinearAlgebra.norm ? abs(J.λ) * sqrt(n) : + norm(Diagonal(fill(J.λ, n))) return norm(A - J) <= max(atol, rtol * max(norm(A), normJ)) end isapprox(A::AbstractMatrix, J::UniformScaling; kwargs...) = isapprox(J, A; kwargs...) From dc5c660c88f7f3e49c6d7d16532d0c9f6ae0ca60 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 9 Jun 2018 21:44:50 +0200 Subject: [PATCH 11/17] rename inner to dot --- stdlib/LinearAlgebra/docs/src/index.md | 1 - stdlib/LinearAlgebra/src/LinearAlgebra.jl | 3 +- stdlib/LinearAlgebra/src/adjtrans.jl | 4 +- stdlib/LinearAlgebra/src/dense.jl | 2 +- stdlib/LinearAlgebra/src/deprecated.jl | 10 ++-- stdlib/LinearAlgebra/src/generic.jl | 64 +++++++++++------------ stdlib/LinearAlgebra/src/hessenberg.jl | 2 +- stdlib/LinearAlgebra/src/lapack.jl | 4 +- stdlib/LinearAlgebra/src/matmul.jl | 10 ++-- stdlib/LinearAlgebra/src/qr.jl | 2 +- stdlib/LinearAlgebra/src/rowvector.jl | 2 +- stdlib/LinearAlgebra/test/adjtrans.jl | 8 +-- stdlib/LinearAlgebra/test/blas.jl | 2 +- stdlib/LinearAlgebra/test/dense.jl | 2 +- stdlib/LinearAlgebra/test/matmul.jl | 30 +++++------ stdlib/SparseArrays/src/SparseArrays.jl | 2 +- stdlib/SparseArrays/src/sparsevector.jl | 24 ++++----- stdlib/SparseArrays/test/sparsevector.jl | 32 ++++++------ stdlib/SuiteSparse/src/spqr.jl | 4 +- test/offsetarray.jl | 2 +- 20 files changed, 103 insertions(+), 107 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index dbabc6a98558d..6e7d122da83cc 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -298,7 +298,6 @@ Linear algebra functions in Julia are largely implemented by calling functions f ```@docs Base.:*(::AbstractMatrix, ::AbstractMatrix) Base.:\(::AbstractMatrix, ::AbstractVecOrMat) -LinearAlgebra.inner LinearAlgebra.dot LinearAlgebra.cross LinearAlgebra.factorize diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 2f27a038f13e3..d841f2fc8f2e5 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -91,7 +91,6 @@ export givens, hessenberg, hessenberg!, - inner, isdiag, ishermitian, isposdef, @@ -371,7 +370,7 @@ include("schur.jl") include("structuredbroadcast.jl") include("deprecated.jl") -const ⋅ = inner +const ⋅ = dot const × = cross export ⋅, × diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index 70fed4c9439ae..d4b08a07502a6 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -194,8 +194,8 @@ broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs... ## multiplication * # Adjoint/Transpose-vector * vector -*(u::AdjointAbsVec, v::AbstractVector) = inner(u.parent, v) -*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = inner(u.parent, v) +*(u::AdjointAbsVec, v::AbstractVector) = dot(u.parent, v) +*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v) function *(u::TransposeAbsVec, v::AbstractVector) @boundscheck length(u) == length(v) || throw(DimensionMismatch()) return sum(@inbounds(u[k]*v[k]) for k in 1:length(u)) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index d7328c8c7f7d5..0f422b70d44cb 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -5,7 +5,7 @@ ## BLAS cutoff threshold constants const SCAL_CUTOFF = 2048 -const INNER_CUTOFF = 128 +#TODO const DOT_CUTOFF = 128 const ASUM_CUTOFF = 32 const NRM2_CUTOFF = 32 diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 3facd40ca9b7a..55a40e175fce2 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -9,8 +9,8 @@ using Base: @deprecate, depwarn # deprecate vecnorm in favor of norm @deprecate vecnorm norm -# deprecate vecdot in favor of inner -@deprecate vecdot inner +# deprecate vecdot in favor of dot +@deprecate vecdot dot # PR #22188 export cholfact, cholfact! @@ -544,9 +544,9 @@ IndexStyle(::Type{<:RowVector}) = IndexLinear() # Multiplication # # inner product -> dot product specializations -@inline *(rowvec::RowVector{T}, vec::AbstractVector{T}) where {T<:Real} = inner(parent(rowvec), vec) -@inline *(rowvec::ConjRowVector{T}, vec::AbstractVector{T}) where {T<:Real} = inner(rvadjoint(rowvec), vec) -@inline *(rowvec::ConjRowVector, vec::AbstractVector) = inner(rvadjoint(rowvec), vec) +@inline *(rowvec::RowVector{T}, vec::AbstractVector{T}) where {T<:Real} = dot(parent(rowvec), vec) +@inline *(rowvec::ConjRowVector{T}, vec::AbstractVector{T}) where {T<:Real} = dot(rvadjoint(rowvec), vec) +@inline *(rowvec::ConjRowVector, vec::AbstractVector) = dot(rvadjoint(rowvec), vec) # Generic behavior @inline function *(rowvec::RowVector, vec::AbstractVector) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 49e6710cad2d2..d627e973a5a4c 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -247,7 +247,7 @@ tril!(M::AbstractMatrix) = tril!(M,0) diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to construct a diagonal matrix")) ########################################################################################### -# Inner products and norms +# Dot products and norms # special cases of norm; note that they don't need to handle isempty(x) function generic_normMinusInf(x) @@ -577,8 +577,8 @@ equivalent to the `p`-norm with value `p = q/(q-1)`. They coincide at `p = q = 2 Use [`norm`](@ref) to compute the `p` norm of `A` as a vector. The difference in norm between a vector space and its dual arises to preserve -the relationship between duality and the inner product, and the result is -consistent with the operator `p`-norm of `1 × n` matrix. +the relationship between duality and the dot product, and the result is +consistent with the operator `p`-norm of a `1 × n` matrix. # Examples ```jldoctest @@ -621,59 +621,60 @@ opnorm(v::TransposeAbsVec) = norm(v.parent) norm(v::Union{TransposeAbsVec,AdjointAbsVec}, p::Real) = norm(v.parent, p) -function inner(x::AbstractArray, y::AbstractArray) +function dot(x::AbstractArray, y::AbstractArray) lx = _length(x) if lx != _length(y) throw(DimensionMismatch("first array has length $(lx) which does not match the length of the second, $(_length(y)).")) end if lx == 0 - return inner(zero(eltype(x)), zero(eltype(y))) + return dot(zero(eltype(x)), zero(eltype(y))) end - s = zero(inner(first(x), first(y))) + s = zero(dot(first(x), first(y))) for (Ix, Iy) in zip(eachindex(x), eachindex(y)) - @inbounds s += inner(x[Ix], y[Iy]) + @inbounds s += dot(x[Ix], y[Iy]) end s end """ - inner(x, y) + dot(x, y) + x ⋅ y For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or -any element type for which `inner` is defined), compute the inner product (or dot product -or scalar product), i.e. the sum of `inner(x[i],y[i])`, as if they were vectors. +any element type for which `dot` is defined), compute the dot product (or inner product +or scalar product), i.e. the sum of `dot(x[i],y[i])`, as if they were vectors. -`dot(x, y)`` and `x ⋅ y` (where `⋅` can be typed by tab-completing `\cdot` in the REPL) are -synonyms for `inner(x, y)`. +`x ⋅ y` (where `⋅` can be typed by tab-completing `\\cdot` in the REPL) is a synonym for +`dot(x, y)`. # Examples ```jldoctest -julia> inner(1:5, 2:6) +julia> dot(1:5, 2:6) 70 julia> x = fill(2., (5,5)); julia> y = fill(3., (5,5)); -julia> inner(x, y) +julia> dot(x, y) 150.0 ``` """ -function inner(x, y) # arbitrary iterables +function dot(x, y) # arbitrary iterables ix = iterate(x) iy = iterate(y) if ix === nothing if iy !== nothing throw(DimensionMismatch("x and y are of different lengths!")) end - return inner(zero(eltype(x)), zero(eltype(y))) + return dot(zero(eltype(x)), zero(eltype(y))) end if iy === nothing throw(DimensionMismatch("x and y are of different lengths!")) end (vx, xs) = ix (vy, ys) = iy - s = inner(vx, vy) + s = dot(vx, vy) while true ix = iterate(x, xs) iy = iterate(y, ys) @@ -681,7 +682,7 @@ function inner(x, y) # arbitrary iterables break end (vx, xs), (vy, ys) = ix, iy - s += inner(vx, vy) + s += dot(vx, vy) end if !(iy == nothing && ix == nothing) throw(DimensionMismatch("x and y are of different lengths!")) @@ -689,48 +690,45 @@ function inner(x, y) # arbitrary iterables return s end -inner(x::Number, y::Number) = conj(x) * y +dot(x::Number, y::Number) = conj(x) * y """ - inner(x,y) dot(x, y) - ⋅(x,y) + x ⋅ y -Compute the inner/dot product between two vectors. For complex vectors, the first -vector is conjugated. When the vectors have equal lengths, calling `inner` is -semantically equivalent to `sum(inner(vx,vy) for (vx,vy) in zip(x, y))`. +Compute the dot product between two vectors. For complex vectors, the first +vector is conjugated. When the vectors have equal lengths, calling `dot` is +semantically equivalent to `sum(dot(vx,vy) for (vx,vy) in zip(x, y))`. # Examples ```jldoctest -julia> inner([1; 1], [2; 3]) +julia> dot([1; 1], [2; 3]) 5 -julia> inner([im; im], [1; 1]) +julia> dot([im; im], [1; 1]) 0 - 2im ``` """ -function inner(x::AbstractVector, y::AbstractVector) +function dot(x::AbstractVector, y::AbstractVector) if length(LinearIndices(x)) != length(LinearIndices(y)) - throw(DimensionMismatch("inner product arguments have unequal lengths $(length(LinearIndices(x))) and $(length(LinearIndices(y)))")) + throw(DimensionMismatch("dot product arguments have unequal lengths $(length(LinearIndices(x))) and $(length(LinearIndices(y)))")) end ix = iterate(x) if ix === nothing # we only need to check the first vector, since equal lengths have been asserted - return inner(zero(eltype(x)), zero(eltype(y))) + return dot(zero(eltype(x)), zero(eltype(y))) end iy = iterate(y) - s = inner(ix[1], iy[1]) + s = dot(ix[1], iy[1]) ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) while ix != nothing - s += inner(ix[1], iy[1]) + s += dot(ix[1], iy[1]) ix = iterate(x, ix[2]) iy = iterate(y, iy[2]) end return s end -const dot = inner - ########################################################################################### diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 19302710180cd..050462c19d55f 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -84,7 +84,7 @@ function getindex(A::HessenbergQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return inner(x, lmul!(A, y)) + return dot(x, lmul!(A, y)) end ## reconstruct the original matrix diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index 4084fe116607d..226c7aca5d222 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -12,7 +12,7 @@ import ..LinearAlgebra.BLAS.@blasfunc import ..LinearAlgebra: BlasFloat, BlasInt, LAPACKException, DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare -using ..LinearAlgebra: triu, tril, inner +using ..LinearAlgebra: triu, tril, dot using Base: iszero @@ -1444,7 +1444,7 @@ for (gglse, elty) in ((:dgglse_, :Float64), resize!(work, lwork) end end - X, inner(view(c, n - p + 1:m), view(c, n - p + 1:m)) + X, dot(view(c, n - p + 1:m), view(c, n - p + 1:m)) end end end diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 0f787cafd7196..8609ef2dc1ba3 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -4,12 +4,12 @@ matprod(x, y) = x*y + x*y -# inner/dot products +# dot products -inner(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y) -inner(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y) +dot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y) +dot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y) -function inner(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasReal,TI<:Integer} +function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasReal,TI<:Integer} if length(rx) != length(ry) throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))")) end @@ -22,7 +22,7 @@ function inner(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vect GC.@preserve x y BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry)) end -function inner(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasComplex,TI<:Integer} +function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector{T}, ry::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasComplex,TI<:Integer} if length(rx) != length(ry) throw(DimensionMismatch("length of rx, $(length(rx)), does not equal length of ry, $(length(ry))")) end diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 97d9151428abc..602c143230cfb 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -469,7 +469,7 @@ function getindex(A::AbstractQ, i::Integer, j::Integer) x[i] = 1 y = zeros(eltype(A), size(A, 2)) y[j] = 1 - return inner(x, lmul!(A, y)) + return dot(x, lmul!(A, y)) end ## Multiplication by Q diff --git a/stdlib/LinearAlgebra/src/rowvector.jl b/stdlib/LinearAlgebra/src/rowvector.jl index 1b63a35d70625..01ea9fe8c35fa 100644 --- a/stdlib/LinearAlgebra/src/rowvector.jl +++ b/stdlib/LinearAlgebra/src/rowvector.jl @@ -12,7 +12,7 @@ recursively). By convention, a vector can be multiplied by a matrix on its left (`A * v`) whereas a row vector can be multiplied by a matrix on its right (such that `RowVector(v) * A = RowVector(transpose(A) * v)`). It differs from a `1×n`-sized matrix by the facts that its transpose returns a vector and the -inner product `RowVector(v1) * v2` returns a scalar, but will otherwise behave similarly. +dot product `RowVector(v1) * v2` returns a scalar, but will otherwise behave similarly. # Examples ```jldoctest diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index 35596e994a6fb..ae140d5f7a8c8 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -342,10 +342,10 @@ end realvec, realmat = [1, 2, 3], [1 2 3; 4 5 6; 7 8 9] complexvec, complexmat = [1im, 2, -3im], [1im 2 3; 4 5 -6im; 7im 8 9] # Adjoint/Transpose-vector * vector - @test Adjoint(realvec) * realvec == inner(realvec, realvec) - @test Transpose(realvec) * realvec == inner(realvec, realvec) - @test Adjoint(complexvec) * complexvec == inner(complexvec, complexvec) - @test Transpose(complexvec) * complexvec == inner(conj(complexvec), complexvec) + @test Adjoint(realvec) * realvec == dot(realvec, realvec) + @test Transpose(realvec) * realvec == dot(realvec, realvec) + @test Adjoint(complexvec) * complexvec == dot(complexvec, complexvec) + @test Transpose(complexvec) * complexvec == dot(conj(complexvec), complexvec) # vector * Adjoint/Transpose-vector @test realvec * Adjoint(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3))) @test realvec * Transpose(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3))) diff --git a/stdlib/LinearAlgebra/test/blas.jl b/stdlib/LinearAlgebra/test/blas.jl index bf5a578cd4f54..d09f35899c358 100644 --- a/stdlib/LinearAlgebra/test/blas.jl +++ b/stdlib/LinearAlgebra/test/blas.jl @@ -53,7 +53,7 @@ srand(100) v41 = convert(Vector{elty}, [4:-1:1;]) let n = 10 - @testset "inner/dot products" begin + @testset "dot products" begin if elty <: Real x1 = convert(Vector{elty}, randn(n)) x2 = convert(Vector{elty}, randn(n)) diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index 21aaff8f0325c..ed192251a958a 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -796,7 +796,7 @@ end r = (elty <: Complex ? adjoint : transpose)(rand(elty, 5)) cm = rand(elty, 5, 1) rm = rand(elty, 1, 5) - @testset "inner products" begin + @testset "dot products" begin test_div_pinv_consistency(r, c) test_div_pinv_consistency(rm, c) test_div_pinv_consistency(r, cm) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index ec06f0964a2f6..f6b47fa9afc6a 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -215,33 +215,33 @@ end end # issue #6450 -@test inner(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 +@test dot(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 -@testset "inner" for elty in (Float32, Float64, ComplexF32, ComplexF64) +@testset "dot" for elty in (Float32, Float64, ComplexF32, ComplexF64) x = convert(Vector{elty},[1.0, 2.0, 3.0]) y = convert(Vector{elty},[3.5, 4.5, 5.5]) - @test_throws DimensionMismatch inner(x, 1:2, y, 1:3) - @test_throws BoundsError inner(x, 1:4, y, 1:4) - @test_throws BoundsError inner(x, 1:3, y, 2:4) - @test inner(x, 1:2, y, 1:2) == convert(elty, 12.5) + @test_throws DimensionMismatch dot(x, 1:2, y, 1:3) + @test_throws BoundsError dot(x, 1:4, y, 1:4) + @test_throws BoundsError dot(x, 1:3, y, 2:4) + @test dot(x, 1:2, y, 1:2) == convert(elty, 12.5) @test transpose(x)*y == convert(elty, 29.0) X = convert(Matrix{elty},[1.0 2.0; 3.0 4.0]) Y = convert(Matrix{elty},[1.5 2.5; 3.5 4.5]) - @test inner(X, Y) == convert(elty, 35.0) + @test dot(X, Y) == convert(elty, 35.0) Z = convert(Vector{Matrix{elty}},[reshape(1:4, 2, 2), fill(1, 2, 2)]) - @test inner(Z, Z) == convert(elty, 34.0) + @test dot(Z, Z) == convert(elty, 34.0) end -inner_(x,y) = invoke(inner, Tuple{Any,Any}, x,y) -@testset "generic inner" begin +dot_(x,y) = invoke(dot, Tuple{Any,Any}, x,y) +@testset "generic dot" begin AA = [1+2im 3+4im; 5+6im 7+8im] BB = [2+7im 4+1im; 3+8im 6+5im] for A in (copy(AA), view(AA, 1:2, 1:2)), B in (copy(BB), view(BB, 1:2, 1:2)) - @test inner(A,B) == inner(vec(A),vec(B)) == inner_(A,B) == inner(float.(A),float.(B)) - @test inner(Int[], Int[]) == 0 == inner_(Int[], Int[]) - @test_throws MethodError inner(Any[], Any[]) - @test_throws MethodError inner_(Any[], Any[]) - for n1 = 0:2, n2 = 0:2, d in (inner, inner_) + @test dot(A,B) == dot(vec(A),vec(B)) == dot_(A,B) == dot(float.(A),float.(B)) + @test dot(Int[], Int[]) == 0 == dot_(Int[], Int[]) + @test_throws MethodError dot(Any[], Any[]) + @test_throws MethodError dot_(Any[], Any[]) + for n1 = 0:2, n2 = 0:2, d in (dot, dot_) if n1 != n2 @test_throws DimensionMismatch d(1:n1, 1:n2) else diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl index 51a62a0417408..af2aa64f6db3d 100644 --- a/stdlib/SparseArrays/src/SparseArrays.jl +++ b/stdlib/SparseArrays/src/SparseArrays.jl @@ -12,7 +12,7 @@ using Base.Sort: Forward using LinearAlgebra import Base: +, -, *, \, /, &, |, xor, == -import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, eigen, inner, +import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, eigen, dot, issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 4a2ac4bb99637..f30ea679b88d2 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -1385,27 +1385,27 @@ end (*)(a::Number, x::SparseVectorUnion) = SparseVector(length(x), copy(nonzeroinds(x)), a * nonzeros(x)) (/)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) / a) -# inner/dot -function inner(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} +# dot +function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} n = length(x) length(y) == n || throw(DimensionMismatch()) nzind = nonzeroinds(y) nzval = nonzeros(y) - s = inner(zero(Tx), zero(Ty)) + s = dot(zero(Tx), zero(Ty)) for i = 1:length(nzind) - s += inner(x[nzind[i]], nzval[i]) + s += dot(x[nzind[i]], nzval[i]) end return s end -function inner(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} +function dot(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} n = length(y) length(x) == n || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) - s = inner(zero(Tx), zero(Ty)) + s = dot(zero(Tx), zero(Ty)) @inbounds for i = 1:length(nzind) - s += inner(nzval[i], y[nzind[i]]) + s += dot(nzval[i], y[nzind[i]]) end return s end @@ -1431,7 +1431,7 @@ function _spdot(f::Function, s end -function inner(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) +function dot(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) x === y && return sum(abs2, x) n = length(x) length(y) == n || throw(DimensionMismatch()) @@ -1441,7 +1441,7 @@ function inner(x::SparseVectorUnion{<:Number}, y::SparseVectorUnion{<:Number}) xnzval = nonzeros(x) ynzval = nonzeros(y) - _spdot(inner, + _spdot(dot, 1, length(xnzind), xnzind, xnzval, 1, length(ynzind), ynzind, ynzval) end @@ -1612,7 +1612,7 @@ mul!(y::AbstractVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractS (A = adjA.parent; mul!(y, adjoint(A), x, one(Tx), zero(Ty))) mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, α::Number, β::Number) = - (A = adjA.parent; _At_or_Ac_mul_B!(inner, y, A, x, α, β)) + (A = adjA.parent; _At_or_Ac_mul_B!(dot, y, A, x, α, β)) function _At_or_Ac_mul_B!(tfun::Function, y::AbstractVector, A::SparseMatrixCSC, x::AbstractSparseVector, @@ -1633,7 +1633,7 @@ function _At_or_Ac_mul_B!(tfun::Function, mx = length(xnzind) for j = 1:n - # s <- inner(A[:,j], x) + # s <- dot(A[:,j], x) s = _spdot(tfun, Acolptr[j], Acolptr[j+1]-1, Arowval, Anzval, 1, mx, xnzind, xnzval) @inbounds y[j] += s * α @@ -1654,7 +1654,7 @@ end (A = transA.parent; _At_or_Ac_mul_B(*, A, x)) *(adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector) = - (A = adjA.parent; _At_or_Ac_mul_B(inner, A, x)) + (A = adjA.parent; _At_or_Ac_mul_B(dot, A, x)) function _At_or_Ac_mul_B(tfun::Function, A::SparseMatrixCSC{TvA,TiA}, x::AbstractSparseVector{TvX,TiX}) where {TvA,TiA,TvX,TiX} m, n = size(A) diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 58dfd24dfb5c0..e2f047a3eb0ee 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -795,14 +795,14 @@ end @test exact_equal(xc, sx) end - @testset "inner/dot" begin - dv = inner(xf, xf2) - @test inner(x, x) == sum(abs2, x) - @test inner(x2, x2) == sum(abs2, x2) - @test inner(x, x2) ≈ dv - @test inner(x2, x) ≈ dv - @test inner(Array(x), x2) ≈ dv - @test inner(x, Array(x2)) ≈ dv + @testset "dot" begin + dv = dot(xf, xf2) + @test dot(x, x) == sum(abs2, x) + @test dot(x2, x2) == sum(abs2, x2) + @test dot(x, x2) ≈ dv + @test dot(x2, x) ≈ dv + @test dot(Array(x), x2) ≈ dv + @test dot(x, Array(x2)) ≈ dv end end @@ -810,8 +810,8 @@ end y = complex.(sprand(32, 0.6), sprand(32, 0.6)) xf = Array(x)::Vector{ComplexF64} yf = Array(y)::Vector{ComplexF64} - @test inner(x, x) ≈ inner(xf, xf) - @test inner(x, y) ≈ inner(xf, yf) + @test dot(x, x) ≈ dot(xf, xf) + @test dot(x, y) ≈ dot(xf, yf) end end @@ -1247,12 +1247,12 @@ end A = sprandn(n, n, 0.01) for j in 1:5:n Aj, Ajview = A[:, j], view(A, :, j) - @test norm(Aj) == norm(Ajview) - @test inner(Aj, copy(Aj)) == inner(Ajview, Aj) # don't alias since it takes a different code path - @test rmul!(Aj, 0.1) == rmul!(Ajview, 0.1) - @test Aj*0.1 == Ajview*0.1 - @test 0.1*Aj == 0.1*Ajview - @test Aj/0.1 == Ajview/0.1 + @test norm(Aj) == norm(Ajview) + @test dot(Aj, copy(Aj)) == dot(Ajview, Aj) # don't alias since it takes a different code path + @test rmul!(Aj, 0.1) == rmul!(Ajview, 0.1) + @test Aj*0.1 == Ajview*0.1 + @test 0.1*Aj == 0.1*Ajview + @test Aj/0.1 == Ajview/0.1 @test LinearAlgebra.axpy!(1.0, Aj, sparse(fill(1., n))) == LinearAlgebra.axpy!(1.0, Ajview, sparse(fill(1., n))) @test LinearAlgebra.lowrankupdate!(Matrix(1.0*I, n, n), fill(1.0, n), Aj) == diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl index b66faa07c3b58..5059f1d75a587 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -206,7 +206,7 @@ function LinearAlgebra.lmul!(Q::QRSparseQ, A::StridedVecOrMat) h = view(Q.factors, :, l) for j in 1:size(A, 2) a = view(A, :, j) - LinearAlgebra.axpy!(τl*inner(h, a), h, a) + LinearAlgebra.axpy!(τl*dot(h, a), h, a) end end return A @@ -236,7 +236,7 @@ function LinearAlgebra.lmul!(adjQ::Adjoint{<:Any,<:QRSparseQ}, A::StridedVecOrMa h = view(Q.factors, :, l) for j in 1:size(A, 2) a = view(A, :, j) - LinearAlgebra.axpy!(τl'*inner(h, a), h, a) + LinearAlgebra.axpy!(τl'*dot(h, a), h, a) end end return A diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 22b7c16ae3bba..2d5c73744e724 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -378,7 +378,7 @@ I = findall(!iszero, z) @test norm(v) ≈ norm(parent(v)) @test norm(A) ≈ norm(parent(A)) -@test inner(v, v) ≈ inner(v0, v0) +@test dot(v, v) ≈ dot(v0, v0) # Prior to its removal from Base, cumsum_kbn was used here. To achieve the same level of # accuracy in the tests, we need to use BigFloats with enlarged precision. From 5b6ea61db8ee55f1c12e6aff5a6b14f6574e7cc8 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 11 Jun 2018 08:17:40 +0200 Subject: [PATCH 12/17] fix norm of RowVector and some docstrings as suggested by @Sacha0 --- stdlib/LinearAlgebra/src/deprecated.jl | 4 ++-- stdlib/LinearAlgebra/src/generic.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 55a40e175fce2..6842c21a895ca 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -1207,8 +1207,8 @@ _mat_ldiv_rowvec_error() = throw(DimensionMismatch("Cannot left-divide matrix by *(A::RowVector, B::Adjoint{<:Any,<:AbstractRotation}) = A * adjoint(B.parent) # methods involving RowVector from base/linalg/generic.jl, to deprecate -norm(tv::RowVector, q::Real) = q == Inf ? opnorm(rvtranspose(tv), 1) : opnorm(rvtranspose(tv), q/(q-1)) -norm(tv::RowVector) = opnorm(rvtranspose(tv)) +norm(tv::RowVector, q::Real) = q == Inf ? norm(rvtranspose(tv), 1) : norm(rvtranspose(tv), q/(q-1)) +norm(tv::RowVector) = norm(rvtranspose(tv)) # methods involving RowVector from base/linalg/factorization.jl, to deprecate \(A::Adjoint{<:Any,<:Factorization}, B::RowVector) = adjoint(A.parent) \ B diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index d627e973a5a4c..bed2605cb75a2 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -372,7 +372,7 @@ The `p`-norm is defined as ```math \\|A\\|_p = \\left( \\sum_{i=1}^n | a_i | ^p \\right)^{1/p} ``` -with ``a_i`` the entries of ``A``, ``| a_i |`` are their [`norm`](@ref)s, and +with ``a_i`` the entries of ``A``, ``| a_i |`` the [`norm`](@ref) of ``a_i``, and ``n`` the length of ``A``. Since the `p`-norm is computed using the [`norm`](@ref)s of the entries of `A`, the `p`-norm of a vector of vectors is not compatible with the interpretation of it as a block vector in general if `p != 2`. @@ -410,7 +410,7 @@ julia> norm([1 2 3; 4 5 6; 7 8 9]) julia> norm([1 2 3 4 5 6 7 8 9]) 16.881943016134134 -julia> norm([1 2 3 4 5 6 7 8 9]') +julia> norm(1:9) 16.881943016134134 julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) != norm([v,v], 1) @@ -515,7 +515,7 @@ end """ opnorm(A::AbstractMatrix, p::Real=2) -Computes the operator norm (or matrix norm) induced by the vector `p`-norm, +Compute the operator norm (or matrix norm) induced by the vector `p`-norm, where valid values of `p` are `1`, `2`, or `Inf`. (Note that for sparse matrices, `p=2` is currently not implemented.) Use [`norm`](@ref) to compute the Frobenius norm. From d6fbe755637d80aec26132135e93a5995e735fb2 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 12 Jun 2018 06:46:06 +0200 Subject: [PATCH 13/17] additional test in 'VecOrMat of Vectors' --- stdlib/LinearAlgebra/test/matmul.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index f6b47fa9afc6a..c96bc0c2061c1 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -425,7 +425,8 @@ end @test transpose(Xv2)*Xv2 ≈ XtX @test (transpose(Xv3)*Xv3)[1] ≈ XtX @test Xv1'*Xv1 ≈ XcX - @test_throws MethodError Xv2'*Xv2 ≈ XcX + @test_throws MethodError Xv2'*Xv2 ≈ XcX # Xv2'*Xv2 is a scalar, XcX a matrix + @test Xv2'*Xv2 ≈ norm(Xv2)^2 @test (Xv3'*Xv3)[1] ≈ XcX @test (Xv1*transpose(Xv1))[1] ≈ XXt @test Xv2*transpose(Xv2) ≈ XXt From 9b24aee638e78351dabe8357b2a52b24ead2245d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 12 Jun 2018 07:31:30 +0200 Subject: [PATCH 14/17] removed a test in 'VecOrMat of Vectors' --- stdlib/LinearAlgebra/test/matmul.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index c96bc0c2061c1..db906303c89ac 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -425,7 +425,6 @@ end @test transpose(Xv2)*Xv2 ≈ XtX @test (transpose(Xv3)*Xv3)[1] ≈ XtX @test Xv1'*Xv1 ≈ XcX - @test_throws MethodError Xv2'*Xv2 ≈ XcX # Xv2'*Xv2 is a scalar, XcX a matrix @test Xv2'*Xv2 ≈ norm(Xv2)^2 @test (Xv3'*Xv3)[1] ≈ XcX @test (Xv1*transpose(Xv1))[1] ≈ XXt From e5127f85f444a65ec1e38f0ab949dcadf7abe5a6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Jun 2018 07:27:44 +0200 Subject: [PATCH 15/17] change vecdot for sparse matrices to dot --- stdlib/SparseArrays/src/linalg.jl | 8 ++++---- stdlib/SparseArrays/test/sparse.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index 8b81eacf42f54..3d6abf879f60f 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -203,11 +203,11 @@ function spmatmul(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}; return C end -# Frobenius inner product: trace(A'B) -function vecdot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T2,S1,S2} +# Frobenius dot/inner product: trace(A'B) +function dot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T2,S1,S2} m, n = size(A) size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions")) - r = vecdot(zero(T1), zero(T2)) + r = dot(zero(T1), zero(T2)) @inbounds for j = 1:n ia = A.colptr[j]; ia_nxt = A.colptr[j+1] ib = B.colptr[j]; ib_nxt = B.colptr[j+1] @@ -223,7 +223,7 @@ function vecdot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T ib < ib_nxt || break rb = B.rowval[ib] else # ra == rb - r += vecdot(A.nzval[ia], B.nzval[ib]) + r += dot(A.nzval[ia], B.nzval[ib]) ia += oneunit(S1); ib += oneunit(S2) ia < ia_nxt && ib < ib_nxt || break ra = A.rowval[ia]; rb = B.rowval[ib] diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index 3e649bf227332..6309efc11b8a0 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -351,13 +351,13 @@ end end end -@testset "sparse Frobenius inner product" begin +@testset "sparse Frobenius dot/inner product" begin for i = 1:5 A = sprand(ComplexF64,10,15,0.4) B = sprand(ComplexF64,10,15,0.5) - @test vecdot(A,B) ≈ vecdot(Matrix(A),Matrix(B)) + @test dot(A,B) ≈ dot(Matrix(A),Matrix(B)) end - @test_throws DimensionMismatch vecdot(sprand(5,5,0.2),sprand(5,6,0.2)) + @test_throws DimensionMismatch dot(sprand(5,5,0.2),sprand(5,6,0.2)) end sA = sprandn(3, 7, 0.5) From 274d59732cb465b036f3c19a64cefb1429bb4e06 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 15 Jun 2018 07:19:29 +0200 Subject: [PATCH 16/17] added news --- NEWS.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/NEWS.md b/NEWS.md index b1592f2d63e9c..6ad4e28c44c86 100644 --- a/NEWS.md +++ b/NEWS.md @@ -511,6 +511,13 @@ This section lists changes that do not have deprecation warnings. This change makes `@schedule` redundant with `@async`, so `@schedule` has been deprecated ([#27164]). + * `norm(A::AbstractMatrix, p=2)` computes no longer the operator/matrix norm but the `norm` of `A` + as for other iterables, i.e. as if it were a vector. Especially, `norm(A::AbstractMatrix)` is the + Frobenius norm. To compute the operator/matrix norm, use the new function `opnorm` ([#27401]). + + * `dot(u, v)` now acts recursively. Instead of `sum(u[i]' * v[i] for i in ...)`, it computes + `sum(dot(u[i], v[i]) for i in ...)`, similarly to `vecdot` before ([#27401]). + Library improvements -------------------- @@ -1274,6 +1281,8 @@ Deprecated or removed * The functions `eigs` and `svds` have been moved to the `Arpack.jl` package ([#27616]). + * `vecdot` and `norm` are deprecated in favor of `dot` and `norm`, respectively ([#27401]). + Command-line option changes --------------------------- @@ -1610,3 +1619,4 @@ Command-line option changes [#27189]: https://github.com/JuliaLang/julia/issues/27189 [#27212]: https://github.com/JuliaLang/julia/issues/27212 [#27248]: https://github.com/JuliaLang/julia/issues/27248 +[#27401]: https://github.com/JuliaLang/julia/issues/27401 From 4282c4f7cd40f641bcb22d2f4cc27931237b4064 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 15 Jun 2018 19:19:01 +0200 Subject: [PATCH 17/17] fixed typo in NEWS and test in givens.jl --- NEWS.md | 2 +- stdlib/LinearAlgebra/test/givens.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6ad4e28c44c86..5e8f35fd77543 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1281,7 +1281,7 @@ Deprecated or removed * The functions `eigs` and `svds` have been moved to the `Arpack.jl` package ([#27616]). - * `vecdot` and `norm` are deprecated in favor of `dot` and `norm`, respectively ([#27401]). + * `vecdot` and `vecnorm` are deprecated in favor of `dot` and `norm`, respectively ([#27401]). Command-line option changes --------------------------- diff --git a/stdlib/LinearAlgebra/test/givens.jl b/stdlib/LinearAlgebra/test/givens.jl index 25a11959812cc..c2fc783547917 100644 --- a/stdlib/LinearAlgebra/test/givens.jl +++ b/stdlib/LinearAlgebra/test/givens.jl @@ -38,7 +38,7 @@ using LinearAlgebra: rmul!, lmul! @test_throws DimensionMismatch lmul!(G, A) @test_throws DimensionMismatch rmul!(A, adjoint(G)) @test abs.(A) ≈ abs.(hessenberg(Ac).H) - @test norm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty) + @test opnorm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty) I10 = Matrix{elty}(I, 10, 10) G, _ = givens(one(elty),zero(elty),9,10)