add vecdot, in analogy with vecnorm, as discussed in #11064

NEWS.md

@@ -156,6 +156,8 @@ Library improvements
 
     * OpenBLAS 64-bit (ILP64) interface is now compiled with a `64_` suffix ([#8734]) to avoid conflicts with external libraries using a 32-bit BLAS ([#4923]).
 
+    * New `vecdot` function, analogous to `vecnorm`, for Euclidean inner products over any iterable container ([#11067]).
+
   * Strings
 
     * NUL-terminated strings should now be passed to C via the new `Cstring` type, not `Ptr{UInt8}` or `Ptr{Cchar}`,
@@ -1393,3 +1395,4 @@ Too numerous to mention.
 [#10893]: https://github.com/JuliaLang/julia/issues/10893
 [#10914]: https://github.com/JuliaLang/julia/issues/10914
 [#10994]: https://github.com/JuliaLang/julia/issues/10994
+[#11067]: https://github.com/JuliaLang/julia/issues/11067

base/exports.jl

@@ -685,6 +685,7 @@ export
     tril,
     triu!,
     triu,
+    vecdot,
     vecnorm,
     ⋅,
     ×,

base/linalg.jl

@@ -115,6 +115,7 @@ export
     triu,
     tril!,
     triu!,
+    vecdot,
     vecnorm,
 
 # Operators

base/linalg/generic.jl

@@ -62,6 +62,9 @@ diag(A::AbstractVector) = error("use diagm instead of diag to construct a diagon
 
 #diagm{T}(v::AbstractVecOrMat{T})
 
+###########################################################################################
+# Inner products and norms
+
 # special cases of vecnorm; note that they don't need to handle isempty(x)
 function generic_vecnormMinusInf(x)
     s = start(x)
@@ -215,6 +218,47 @@ end
 norm(x::Number, p::Real=2) =
     p == 0 ? convert(typeof(real(x)), ifelse(x != 0, 1, 0)) : abs(x)
 
+function vecdot(x::AbstractVector, y::AbstractVector)
+    lx = length(x)
+    lx==length(y) || throw(DimensionMismatch())
+    if lx == 0
+        return dot(zero(eltype(x)), zero(eltype(y)))
+    end
+    s = dot(x[1], y[1])
+    @inbounds for i = 2:lx
+        s += dot(x[i], y[i])
+    end
+    s
+end
+
+function vecdot(x, y) # arbitrary iterables
+    ix = start(x)
+    if done(x, ix)
+        isempty(y) || throw(DimensionMismatch())
+        return dot(zero(eltype(x)), zero(eltype(y)))
+    end
+    iy = start(y)
+    done(y, iy) && throw(DimensionMismatch())
+    (vx, ix) = next(x, ix)
+    (vy, iy) = next(y, iy)
+    s = dot(vx, vy)
+    while !done(x, ix)
+        done(y, iy) && throw(DimensionMismatch())
+        (vx, ix) = next(x, ix)
+        (vy, iy) = next(y, iy)
+        s += dot(vx, vy)
+    end
+    !done(y, iy) && throw(DimensionMismatch())
+    return s
+end
+
+vecdot(x::Number, y::Number) = conj(x) * y
+
+dot(x::Number, y::Number) = vecdot(x, y)
+dot(x::AbstractVector, y::AbstractVector) = vecdot(x, y)
+
+###########################################################################################
+
 rank(A::AbstractMatrix, tol::Real) = sum(svdvals(A) .> tol)
 function rank(A::AbstractMatrix)
     m,n = size(A)

base/linalg/matmul.jl

@@ -33,8 +33,8 @@ scale(b::Vector, A::Matrix) = scale!(similar(b, promote_type(eltype(A),eltype(b)
 
 # Dot products
 
-dot{T<:BlasReal}(x::StridedVector{T}, y::StridedVector{T}) = BLAS.dot(x, y)
-dot{T<:BlasComplex}(x::StridedVector{T}, y::StridedVector{T}) = BLAS.dotc(x, y)
+vecdot{T<:BlasReal}(x::Union(DenseArray{T},StridedVector{T}), y::Union(DenseArray{T},StridedVector{T})) = BLAS.dot(x, y)
+vecdot{T<:BlasComplex}(x::Union(DenseArray{T},StridedVector{T}), y::Union(DenseArray{T},StridedVector{T})) = BLAS.dotc(x, y)
 function dot{T<:BlasReal, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},Range{TI}), y::Vector{T}, ry::Union(UnitRange{TI},Range{TI}))
     length(rx)==length(ry) || throw(DimensionMismatch())
     if minimum(rx) < 1 || maximum(rx) > length(x) || minimum(ry) < 1 || maximum(ry) > length(y)
@@ -49,19 +49,7 @@ function dot{T<:BlasComplex, TI<:Integer}(x::Vector{T}, rx::Union(UnitRange{TI},
     end
     BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
 end
-function dot(x::AbstractVector, y::AbstractVector)
-    lx = length(x)
-    lx==length(y) || throw(DimensionMismatch())
-    if lx == 0
-        return zero(eltype(x))*zero(eltype(y))
-    end
-    s = conj(x[1])*y[1]
-    @inbounds for i = 2:lx
-        s += conj(x[i])*y[i]
-    end
-    s
-end
-dot(x::Number, y::Number) = conj(x) * y
+
 Ac_mul_B(x::AbstractVector, y::AbstractVector) = [dot(x, y)]
 At_mul_B{T<:Real}(x::AbstractVector{T}, y::AbstractVector{T}) = [dot(x, y)]
 At_mul_B{T<:BlasComplex}(x::StridedVector{T}, y::StridedVector{T}) = [BLAS.dotu(x, y)]

doc/stdlib/linalg.rst

@@ -30,6 +30,13 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    Compute the dot product. For complex vectors, the first vector is conjugated.
 
+.. function:: vecdot(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.
+
 .. function:: cross(x, y)
               ×(x,y)
 
@@ -493,9 +500,10 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
 .. function:: vecnorm(A, [p])
 
-   For any iterable container ``A`` (including arrays of any dimension)
-   of numbers, compute the ``p``-norm (defaulting to ``p=2``) as if ``A``
-   were a vector of the corresponding length.
+   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
+   ``A`` were a vector of the corresponding length.
 
    For example, if ``A`` is a matrix and ``p=2``, then this is equivalent
    to the Frobenius norm.

test/linalg3.jl

@@ -275,3 +275,18 @@ A = [ 1  5  9
 @test diag(zeros(0,1),1) == []
 @test_throws BoundsError diag(zeros(0,1),-1)
 @test_throws BoundsError diag(zeros(0,1),2)
+
+vecdot_(x,y) = invoke(vecdot, (Any,Any), x,y) # generic vecdot
+let A = [1+2im 3+4im; 5+6im 7+8im], B = [2+7im 4+1im; 3+8im 6+5im]
+    @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_)
+        if n1 != n2
+            @test_throws DimensionMismatch d(1:n1, 1:n2)
+        else
+            @test_approx_eq d(1:n1, 1:n2) vecnorm(1:n1)^2
+        end
+    end
+end