Merge branch 'master' into jishnub/genmatvecmul_onepass

Author: jishnub

src/LinearAlgebra.jl

@@ -840,7 +840,10 @@ function lbt_openblas_onload_callback()
     end
 end
 
-# If users want to lazily load a different BLAS, they'd need to either change this call, or
-# clear the datastructures modified by this call and call it again with their own.
-libblastrampoline_jll.add_dependency!(OpenBLAS_jll, libopenblas, lbt_openblas_onload_callback)
+function __init__()
+    # If users want to lazily load a different BLAS, they'd need to either change this call, or
+    # clear the datastructures modified by this call and call it again with their own.
+    libblastrampoline_jll.add_dependency!(OpenBLAS_jll, libopenblas, lbt_openblas_onload_callback)
+end
+
 end # module LinearAlgebra

src/bunchkaufman.jl

@@ -182,10 +182,8 @@ julia> d, u, p = S; # destructuring via iteration
 julia> d == S.D && u == S.U && p == S.p
 true
 
-julia> S.U*S.D*S.U' - S.P*A*S.P'
-2×2 Matrix{Float64}:
- 0.0  0.0
- 0.0  0.0
+julia> S.U * S.D * S.U' ≈ S.P * A * S.P'
+true
 
 julia> S = bunchkaufman(Symmetric(A, :L))
 BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}
@@ -202,10 +200,8 @@ permutation:
  2
  1
 
-julia> S.L*S.D*S.L' - A[S.p, S.p]
-2×2 Matrix{Float64}:
- 0.0  0.0
- 0.0  0.0
+julia> S.L * S.D * S.L' ≈ A[S.p, S.p]
+true
 ```
 """
 bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =

src/matmul.jl

@@ -692,7 +692,7 @@ end
 # the aggressive constprop pushes tA and tB into gemm_wrapper!, which is needed for wrap calls within it
 # to be concretely inferred
 Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T},
-        alpha::Number, beta::Number) where {T<:BlasFloat}
+        α::Number, β::Number) where {T<:BlasFloat}
     nC = checksquare(C)
     tA_uc = uppercase(tA) # potentially convert a WrapperChar to a Char
     if tA_uc == 'T'
@@ -708,16 +708,16 @@ Base.@constprop :aggressive function syrk_wrapper!(C::StridedMatrix{T}, tA::Abst
 
     # BLAS.syrk! only updates symmetric C
     # alternatively, make non-zero β a show-stopper for BLAS.syrk!
-    if iszero(beta) || issymmetric(C)
-        α, β = promote(alpha, beta, zero(T))
+    if iszero(β) || issymmetric(C)
+        alpha, beta = promote(α, β, zero(T))
         if (alpha isa Union{Bool,T} &&
                 beta isa Union{Bool,T} &&
                 stride(A, 1) == stride(C, 1) == 1 &&
                 _fullstride2(A) && _fullstride2(C))
             return copytri!(BLAS.syrk!('U', tA, alpha, A, beta, C), 'U')
         end
     end
-    return gemm_wrapper!(C, tA, tAt, A, A, alpha, beta)
+    return gemm_wrapper!(C, tA, tAt, A, A, α, β)
 end
 # legacy method
 syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T}, _add::MulAddMul = MulAddMul()) where {T<:BlasFloat} =
@@ -743,7 +743,7 @@ Base.@constprop :aggressive function herk_wrapper!(C::Union{StridedMatrix{T}, St
     # Result array does not need to be initialized as long as beta==0
     #    C = Matrix{T}(undef, mA, mA)
 
-    if iszero(β) || issymmetric(C)
+    if iszero(β) || ishermitian(C)
         alpha, beta = promote(α, β, zero(T))
         if (alpha isa Union{Bool,T} &&
                 beta isa Union{Bool,T} &&
@@ -932,7 +932,7 @@ function __generic_matvecmul!(f::F, C::AbstractVector, A::AbstractVecOrMat, B::A
     @inbounds begin
         if length(B) == 0
             for k = eachindex(C)
-                @stable_muladdmul _modify!(MulAddMul(alpha,beta), false, C, k)
+                @stable_muladdmul _modify!(MulAddMul(alpha,beta), zero(eltype(C)), C, k)
             end
         else
             for k = eachindex(C)
@@ -953,7 +953,7 @@ function __generic_matvecmul!(::typeof(identity), C::AbstractVector, A::Abstract
     @inbounds begin
         for i in eachindex(C)
             if length(B) == 0
-                C[i] = false
+                C[i] = zero(eltype(C))
             else
                 @stable_muladdmul _modify!(MulAddMul(alpha, beta), A[i]*B[1], C, i)
             end

src/symmetric.jl

@@ -352,7 +352,7 @@ function copyto!(dest::Symmetric, src::Symmetric)
     elseif src.uplo == dest.uplo
         copytrito!(dest.data, src.data, src.uplo)
     else
-        transpose!(dest.data, Base.unalias(dest.data, src.data))
+        copytrito!(dest.data, transpose(Base.unalias(dest.data, src.data)), dest.uplo)
     end
     return dest
 end
@@ -363,7 +363,7 @@ function copyto!(dest::Hermitian, src::Hermitian)
     elseif src.uplo == dest.uplo
         copytrito!(dest.data, src.data, src.uplo)
     else
-        adjoint!(dest.data, Base.unalias(dest.data, src.data))
+        copytrito!(dest.data, adjoint(Base.unalias(dest.data, src.data)), dest.uplo)
     end
     return dest
 end
@@ -445,8 +445,6 @@ adjoint(A::Hermitian) = A
 transpose(A::Symmetric) = A
 adjoint(A::Symmetric{<:Real}) = A
 transpose(A::Hermitian{<:Real}) = A
-adjoint(A::Symmetric) = Adjoint(A)
-transpose(A::Hermitian) = Transpose(A)
 
 real(A::Symmetric{<:Real}) = A
 real(A::Hermitian{<:Real}) = A

test/matmul.jl

@@ -6,6 +6,11 @@ using Base: rtoldefault
 using Test, LinearAlgebra, Random
 using LinearAlgebra: mul!, Symmetric, Hermitian
 
+const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+
+isdefined(Main, :SizedArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "SizedArrays.jl"))
+using .Main.SizedArrays
+
 ## Test Julia fallbacks to BLAS routines
 
 mul_wrappers = [
@@ -1176,4 +1181,16 @@ end
     end
 end
 
+@testset "zero-length generic matvec" begin
+    m = SizedArrays.SizedArray{(2,2)}(ones(2,2))
+    A = fill(m, 2, 0)
+    v = fill(m, size(A,2))
+    w = similar(v, size(A,1))
+    mul!(w, A, v)
+    @test all(iszero, w)
+    A = fill(m, 0, 2)
+    mul!(w, A', v)
+    @test all(iszero, w)
+end
+
 end # module TestMatmul

test/symmetric.jl

@@ -1163,4 +1163,19 @@ end
     end
 end
 
+@testset "copyto! with mismatched uplo" begin
+    for (S,tf) in ((Symmetric, transpose), (Hermitian, adjoint))
+        for (uplo1,uplo2) in [(:U,:L), (:L, :U)]
+            M = Matrix{Complex{BigFloat}}(undef, 2, 2)
+            M[1,1] = M[2,2] = 3
+            isupper = uplo1 == :U
+            M[1+!isupper, 1+isupper] = 4+3im
+            H1 = S(M, uplo1)
+            H2 = 2 * S(tf(M), uplo2)
+            copyto!(H2, H1)
+            @test H2 == H1
+        end
+    end
+end
+
 end # module TestSymmetric