From aeb2a38e8e258e7d3139a90c74a83f0cfeaaf6c0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 20 Jun 2018 07:14:50 +0200 Subject: [PATCH 1/2] add unconjugated dot product: dotu --- stdlib/LinearAlgebra/docs/src/index.md | 1 + stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 + stdlib/LinearAlgebra/src/generic.jl | 72 +++++++++++++++++++++++ stdlib/LinearAlgebra/src/matmul.jl | 29 +++++++++ stdlib/LinearAlgebra/test/matmul.jl | 40 +++++++++++++ stdlib/SparseArrays/src/SparseArrays.jl | 2 +- stdlib/SparseArrays/src/linalg.jl | 30 ++++++++++ stdlib/SparseArrays/test/sparse.jl | 2 + test/offsetarray.jl | 1 + 9 files changed, 177 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 6e7d122da83cc..9e04a474a1660 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/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 diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index d841f2fc8f2e5..d51d732164c0a 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -80,6 +80,7 @@ export diagind, diagm, dot, + dotu, eigen, eigen!, eigmax, diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index bed2605cb75a2..3c476c03f2df1 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -729,6 +729,78 @@ 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 + s = ix[1] * iy[1] + ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) + while ix != nothing && iy != nothing + s += ix[1] * iy[1] + ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) + 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 + ########################################################################################### diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 8609ef2dc1ba3..a5493483f5f23 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/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) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index db906303c89ac..7853eaef540b2 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/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) diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl index af2aa64f6db3d..7dabb35363a05 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, 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 diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index 3d6abf879f60f..6f88439e03771 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/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 diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index 6309efc11b8a0..89536a5203fbd 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/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) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 2d5c73744e724..6e87482e3f22b 100644 --- a/test/offsetarray.jl +++ b/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. From 11078ccf63af6d3dfbd02400e443d7e04e42f528 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 21 Jun 2018 09:38:20 +0200 Subject: [PATCH 2/2] improve performance of generic dotu similar to dot, cf. #27678 --- stdlib/LinearAlgebra/src/generic.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 3c476c03f2df1..caee8a92d80cb 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -763,22 +763,27 @@ julia> dotu(σ, σ) function dotu(x, y) # arbitrary iterables ix = iterate(x) iy = iterate(y) - if ix == nothing - if iy != nothing + 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 + if iy === nothing throw(DimensionMismatch("x and y are of different lengths!")) end - s = ix[1] * iy[1] - ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) - while ix != nothing && iy != nothing - s += ix[1] * iy[1] - ix, iy = iterate(x, ix[2]), iterate(y, iy[2]) + (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) + if !(iy === nothing && ix === nothing) throw(DimensionMismatch("x and y are of different lengths!")) end return s