Add unconjugated dot product dotu

stdlib/LinearAlgebra/docs/src/index.md

@@ -299,6 +299,7 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 Base.:*(::AbstractMatrix, ::AbstractMatrix)
 Base.:\(::AbstractMatrix, ::AbstractVecOrMat)
 LinearAlgebra.dot
+LinearAlgebra.dotu
 LinearAlgebra.cross
 LinearAlgebra.factorize
 LinearAlgebra.Diagonal

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -80,6 +80,7 @@ export
     diagind,
     diagm,
     dot,
+    dotu,
     eigen,
     eigen!,
     eigmax,

stdlib/LinearAlgebra/src/generic.jl

@@ -729,6 +729,83 @@ function dot(x::AbstractVector, y::AbstractVector)
     return s
 end
 
+"""
+    dotu(x, y)
+
+For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or
+any element type for which `*` is defined), compute the unconjugated dot product, i.e. the
+sum of `x[i]*y[i]`, as if they were vectors.
+
+# Examples
+```jldoctest
+julia> dotu(1:5, 2:6)
+70
+
+julia> v = [1, im]; dotu(v, v)
+0 + 0im
+
+julia> σ = [[0 1; 1 0], [0 -im; im 0], [1 0; 0 -1]]; n = [1, 2, 3]; dotu(σ, n)
+2×2 Array{Complex{Int64},2}:
+ 3+0im   1-2im
+ 1+2im  -3+0im
+
+julia> dotu(σ[1:1], σ[1:1])
+2×2 Array{Complex{Int64},2}:
+ 1+0im  0+0im
+ 0+0im  1+0im
+
+julia> dotu(σ, σ)
+ 2×2 Array{Complex{Int64},2}:
+  3+0im  0+0im
+  0+0im  3+0im
+```
+"""
+function dotu(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 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 = vx * vy
+    while true
+        ix = iterate(x, xs)
+        iy = iterate(y, ys)
+        ix === nothing && break
+        iy === nothing && break
+        (vx, xs), (vy, ys) = ix, iy
+        s += vx * vy
+    end
+    if !(iy === nothing && ix === nothing)
+        throw(DimensionMismatch("x and y are of different lengths!"))
+    end
+    return s
+end
+
+dotu(x::Number, y::Number) = x * y
+
+function dotu(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 zero(eltype(x)) * zero(eltype(y))
+    end
+    s = zero(first(x) * first(y))
+    @inbounds for (Ix, Iy) in zip(eachindex(x), eachindex(y))
+        s += x[Ix] * y[Iy]
+    end
+    s
+end
+
 
 ###########################################################################################
 

stdlib/LinearAlgebra/src/matmul.jl

@@ -35,6 +35,35 @@ function dot(x::Vector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}, y::Vector
     GC.@preserve x y BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
 end
 
+dotu(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y)
+dotu(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotu(x, y)
+
+function dotu(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
+    if minimum(rx) < 1 || maximum(rx) > length(x)
+        throw(BoundsError(x, rx))
+    end
+    if minimum(ry) < 1 || maximum(ry) > length(y)
+        throw(BoundsError(y, ry))
+    end
+    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 dotu(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
+    if minimum(rx) < 1 || maximum(rx) > length(x)
+        throw(BoundsError(x, rx))
+    end
+    if minimum(ry) < 1 || maximum(ry) > length(y)
+        throw(BoundsError(y, ry))
+    end
+    GC.@preserve x y BLAS.dotu(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
+end
+
 function *(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) where {T<:BlasComplex}
     x = transx.parent
     return BLAS.dotu(x, y)

stdlib/LinearAlgebra/test/matmul.jl

@@ -251,6 +251,46 @@ dot_(x,y) = invoke(dot, Tuple{Any,Any}, x,y)
     end
 end
 
+@testset "dotu" 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 dotu(x, 1:2, y, 1:3)
+    @test_throws BoundsError dotu(x, 1:4, y, 1:4)
+    @test_throws BoundsError dotu(x, 1:3, y, 2:4)
+    @test dotu(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 dotu(X, Y) == convert(elty, 35.0)
+    Z = convert(Vector{Matrix{elty}},[reshape(1:4, 2, 2), fill(1, 2, 2)])
+    @test dotu(Z, Z) == convert(Matrix{elty},[9 17; 12 24])
+    @test dotu(one(elty), one(elty)) == one(elty) == dotu(ones(elty, 1), ones(elty, 1))
+    @test dotu(im*one(elty), one(elty)) == im*one(elty) == dotu(im*ones(elty, 1), ones(elty, 1))
+    @test dotu(one(elty), im*one(elty)) == im*one(elty) == dotu(ones(elty, 1), im*ones(elty, 1))
+end
+@test dotu(Any[1.0,2.0], Any[3.5,4.5]) === 12.5
+
+dotu1(x,y) = invoke(dotu, Tuple{Any,Any}, x,y)
+dotu2(x,y) = invoke(dotu, Tuple{AbstractArray,AbstractArray}, x,y)
+@testset "generic dotu" 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 dotu(A,B) == dotu(vec(A),vec(B)) == dotu1(A,B) == dotu2(A,B) == dotu(float.(A),float.(B)) == sum(A .* B)
+        @test dotu(Int[], Int[]) == 0 == dotu1(Int[], Int[]) == dotu2(Int[], Int[])
+        @test_throws MethodError dotu(Any[], Any[])
+        @test_throws MethodError dotu1(Any[], Any[])
+        @test_throws MethodError dotu2(Any[], Any[])
+        for n1 = 0:2, n2 = 0:2, d in (dotu, dotu1, dotu2)
+            if n1 != n2
+                @test_throws DimensionMismatch d(1:n1, 1:n2)
+            else
+                @test d(1:n1, 1:n2) ≈ norm(1:n1)^2
+            end
+        end
+    end
+end
+
 @testset "Issue 11978" begin
     A = Matrix{Matrix{Float64}}(undef, 2, 2)
     A[1,1] = Matrix(1.0I, 3, 3)

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, dot,
+import LinearAlgebra: mul!, ldiv!, rdiv!, chol, adjoint!, diag, eigen, dot, dotu,
     issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!,
     cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu
 

stdlib/SparseArrays/src/linalg.jl

@@ -234,6 +234,36 @@ function dot(A::SparseMatrixCSC{T1,S1},B::SparseMatrixCSC{T2,S2}) where {T1,T2,S
     return r
 end
 
+function dotu(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 = 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]
+        if ia < ia_nxt && ib < ib_nxt
+            ra = A.rowval[ia]; rb = B.rowval[ib]
+            while true
+                if ra < rb
+                    ia += oneunit(S1)
+                    ia < ia_nxt || break
+                    ra = A.rowval[ia]
+                elseif ra > rb
+                    ib += oneunit(S2)
+                    ib < ib_nxt || break
+                    rb = B.rowval[ib]
+                else # ra == rb
+                    r += 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]
+                end
+            end
+        end
+    end
+    return r
+end
+
 ## solvers
 function fwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat)
 # forward substitution for CSC matrices

stdlib/SparseArrays/test/sparse.jl

@@ -356,8 +356,10 @@ end
         A = sprand(ComplexF64,10,15,0.4)
         B = sprand(ComplexF64,10,15,0.5)
         @test dot(A,B) ≈ dot(Matrix(A),Matrix(B))
+        @test dotu(A,B) ≈ dotu(Matrix(A),Matrix(B))
     end
     @test_throws DimensionMismatch dot(sprand(5,5,0.2),sprand(5,6,0.2))
+    @test_throws DimensionMismatch dotu(sprand(5,5,0.2),sprand(5,6,0.2))
 end
 
 sA = sprandn(3, 7, 0.5)

test/offsetarray.jl

@@ -379,6 +379,7 @@ I = findall(!iszero, z)
 @test norm(v) ≈ norm(parent(v))
 @test norm(A) ≈ norm(parent(A))
 @test dot(v, v) ≈ dot(v0, v0)
+@test dotu(v, v) ≈ dotu(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.