Merge branch 'master' into sgj/evalpoly

Author: LilithHafner

src/LinearAlgebra.jl

@@ -10,7 +10,7 @@ module LinearAlgebra
 import Base: \, /, //, *, ^, +, -, ==
 import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, asec, asech,
     asin, asinh, atan, atanh, axes, big, broadcast, cbrt, ceil, cis, collect, conj, convert,
-    copy, copyto!, copymutable, cos, cosh, cot, coth, csc, csch, eltype, exp, fill!, floor,
+    copy, copy!, copyto!, copymutable, cos, cosh, cot, coth, csc, csch, eltype, exp, fill!, floor,
     getindex, hcat, getproperty, imag, inv, invpermuterows!, isapprox, isequal, isone, iszero,
     IndexStyle, kron, kron!, length, log, map, ndims, one, oneunit, parent, permutecols!,
     permutedims, permuterows!, power_by_squaring, promote_rule, real, isreal, sec, sech, setindex!,
@@ -708,9 +708,9 @@ function ldiv(F::Factorization, B::AbstractVecOrMat)
 
     if n > size(B, 1)
         # Underdetermined
-        copyto!(view(BB, 1:m, :), B)
+        copy!(view(BB, axes(B,1), ntuple(_->:, ndims(B)-1)...), B)
     else
-        copyto!(BB, B)
+        copy!(BB, B)
     end
 
     ldiv!(FF, BB)
@@ -790,31 +790,22 @@ function versioninfo(io::IO=stdout)
     println(io, indent, "LinearAlgebra.BLAS.get_num_threads() = ", BLAS.get_num_threads())
     println(io, "Relevant environment variables:")
     env_var_names = [
-        "JULIA_NUM_THREADS",
-        "MKL_DYNAMIC",
-        "MKL_NUM_THREADS",
+        ["JULIA_NUM_THREADS"],
+        ["MKL_DYNAMIC"],
+        ["MKL_NUM_THREADS"],
          # OpenBLAS has a hierarchy of environment variables for setting the
          # number of threads, see
          # https://github.com/xianyi/OpenBLAS/blob/c43ec53bdd00d9423fc609d7b7ecb35e7bf41b85/README.md#setting-the-number-of-threads-using-environment-variables
-        ("OPENBLAS_NUM_THREADS", "GOTO_NUM_THREADS", "OMP_NUM_THREADS"),
+        ["OPENBLAS_NUM_THREADS", "GOTO_NUM_THREADS", "OMP_NUM_THREADS"],
     ]
     printed_at_least_one_env_var = false
     print_var(io, indent, name) = println(io, indent, name, " = ", ENV[name])
     for name in env_var_names
-        if name isa Tuple
-            # If `name` is a Tuple, then find the first environment which is
-            # defined, and disregard the following ones.
-            for nm in name
-                if haskey(ENV, nm)
-                    print_var(io, indent, nm)
-                    printed_at_least_one_env_var = true
-                    break
-                end
-            end
-        else
-            if haskey(ENV, name)
-                print_var(io, indent, name)
+        for nm in name
+            if haskey(ENV, nm)
+                print_var(io, indent, nm)
                 printed_at_least_one_env_var = true
+                break
             end
         end
     end

src/bidiag.jl

@@ -538,8 +538,9 @@ function lmul!(D::Diagonal, B::Bidiagonal)
     matmul_size_check(size(D), size(B))
     (; dv, ev) = B
     isL = B.uplo == 'L'
-    dv[1] = D.diag[1] * dv[1]
-    for i in axes(ev,1)
+    iszero(size(D,1)) && return B
+    @inbounds dv[1] = D.diag[1] * dv[1]
+    @inbounds for i in axes(ev,1)
         ev[i] = D.diag[i + isL] * ev[i]
         dv[i+1] = D.diag[i+1] * dv[i+1]
     end
@@ -575,8 +576,9 @@ function rmul!(B::Bidiagonal, D::Diagonal)
     matmul_size_check(size(B), size(D))
     (; dv, ev) = B
     isU = B.uplo == 'U'
-    dv[1] *= D.diag[1]
-    for i in axes(ev,1)
+    iszero(size(D,1)) && return B
+    @inbounds dv[1] *= D.diag[1]
+    @inbounds for i in axes(ev,1)
         ev[i] *= D.diag[i + isU]
         dv[i+1] *= D.diag[i+1]
     end

src/dense.jl

@@ -720,7 +720,7 @@ julia> ℯ^[1 2; 0 3]
  0.0      20.0855
 ```
 """
-Base.:^(b::Number, A::AbstractMatrix) = exp!(log(b)*A)
+Base.:^(b::Number, A::AbstractMatrix) = exp_maybe_inplace(log(b)*A)
 # method for ℯ to explicitly elide the log(b) multiplication
 Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A)
 

src/diagonal.jl

@@ -191,7 +191,7 @@ end
 Return the appropriate zero element `A[i, j]` corresponding to a banded matrix `A`.
 """
 diagzero(A::AbstractMatrix, i, j) = zero(eltype(A))
-diagzero(A::AbstractMatrix{M}, i, j) where {M<:AbstractMatrix} =
+@propagate_inbounds diagzero(A::AbstractMatrix{M}, i, j) where {M<:AbstractMatrix} =
     zeroslike(M, axes(A[i,i], 1), axes(A[j,j], 2))
 diagzero(A::AbstractMatrix, inds...) = diagzero(A, to_indices(A, inds)...)
 # dispatching on the axes permits specializing on the axis types to return something other than an Array

src/matmul.jl

@@ -725,8 +725,9 @@ syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T}, _add
 
 # 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 herk_wrapper!(C::Union{StridedMatrix{T}, StridedMatrix{Complex{T}}}, tA::AbstractChar, A::Union{StridedVecOrMat{T}, StridedVecOrMat{Complex{T}}},
-        α::Number, β::Number) where {T<:BlasReal}
+Base.@constprop :aggressive function herk_wrapper!(C::StridedMatrix{TC}, tA::AbstractChar, A::StridedVecOrMat{TC},
+        α::Number, β::Number) where {TC<:BlasComplex}
+    T = real(TC)
     nC = checksquare(C)
     tA_uc = uppercase(tA) # potentially convert a WrapperChar to a Char
     if tA_uc == 'C'
@@ -740,13 +741,10 @@ Base.@constprop :aggressive function herk_wrapper!(C::Union{StridedMatrix{T}, St
         throw(DimensionMismatch(lazy"output matrix has size: $(size(C)), but should have size $((mA, mA))"))
     end
 
-    # Result array does not need to be initialized as long as beta==0
-    #    C = Matrix{T}(undef, mA, mA)
-
+    # BLAS.herk! only updates hermitian C, alpha and beta need to be real
     if iszero(β) || ishermitian(C)
         alpha, beta = promote(α, β, zero(T))
-        if (alpha isa Union{Bool,T} &&
-                beta isa Union{Bool,T} &&
+        if (alpha isa T && beta isa T &&
                 stride(A, 1) == stride(C, 1) == 1 &&
                 _fullstride2(A) && _fullstride2(C))
             return copytri!(BLAS.herk!('U', tA, alpha, A, beta, C), 'U', true)

src/triangular.jl

@@ -178,7 +178,7 @@ end
 parent(A::UpperOrLowerTriangular) = A.data
 
 # For strided matrices, we may only loop over the filled triangle
-copy(A::UpperOrLowerTriangular{<:Any, <:StridedMaybeAdjOrTransMat}) = copyto!(similar(A), A)
+copy(A::UpperOrLowerTriangular{<:Any, <:StridedMaybeAdjOrTransMat}) = copy!(similar(A), A)
 
 # then handle all methods that requires specific handling of upper/lower and unit diagonal
 
@@ -532,30 +532,34 @@ for T in (:UpperOrUnitUpperTriangular, :LowerOrUnitLowerTriangular)
         if axes(dest) != axes(U)
             @invoke copyto!(dest::AbstractArray, U::AbstractArray)
         else
-            _copyto!(dest, U)
+            copy!(dest, U)
         end
         return dest
     end
+    @eval function copy!(dest::$T, U::$T)
+        axes(dest) == axes(U) || throw(ArgumentError(
+            "arrays must have the same axes for copy! (consider using `copyto!`)"))
+        _copy!(dest, U)
+    end
 end
 
 # copy and scale
 for (T, UT) in ((:UpperTriangular, :UnitUpperTriangular), (:LowerTriangular, :UnitLowerTriangular))
-    @eval @inline function _copyto!(A::$T, B::$T)
-        @boundscheck checkbounds(A, axes(B)...)
+    @eval @inline function _copy!(A::$T, B::$T)
         copytrito!(parent(A), parent(B), uplo_char(A))
         return A
     end
-    @eval @inline function _copyto!(A::$UT, B::$T)
+    @eval @inline function _copy!(A::$UT, B::$T)
         for dind in diagind(A, IndexStyle(A))
             if A[dind] != B[dind]
                 throw_nononeerror(typeof(A), B[dind], Tuple(dind)...)
             end
         end
-        _copyto!($T(parent(A)), B)
+        _copy!($T(parent(A)), B)
         return A
     end
 end
-@inline function _copyto!(A::UpperOrUnitUpperTriangular, B::UnitUpperTriangular)
+@inline function _copy!(A::UpperOrUnitUpperTriangular, B::UnitUpperTriangular)
     @boundscheck checkbounds(A, axes(B)...)
     B2 = Base.unalias(A, B)
     Ap = parent(A)
@@ -570,7 +574,7 @@ end
     end
     return A
 end
-@inline function _copyto!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular)
+@inline function _copy!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular)
     @boundscheck checkbounds(A, axes(B)...)
     B2 = Base.unalias(A, B)
     Ap = parent(A)
@@ -595,23 +599,30 @@ _triangularize!(::LowerOrUnitLowerTriangular) = tril!
     if axes(dest) != axes(U)
         @invoke copyto!(dest::StridedMatrix, U::AbstractArray)
     else
-        _copyto!(dest, U)
+        copy!(dest, U)
     end
     return dest
 end
-@propagate_inbounds function _copyto!(dest::StridedMatrix, U::UpperOrLowerTriangular)
+
+function copy!(dest::StridedMatrix, U::UpperOrLowerTriangular)
+    axes(dest) == axes(U) || throw(ArgumentError(
+            "arrays must have the same axes for copy! (consider using `copyto!`)"))
+    _copy!(dest, U)
+end
+
+@propagate_inbounds function _copy!(dest::StridedMatrix, U::UpperOrLowerTriangular)
     copytrito!(dest, parent(U), U isa UpperOrUnitUpperTriangular ? 'U' : 'L')
     copytrito!(dest, U, U isa UpperOrUnitUpperTriangular ? 'L' : 'U')
     return dest
 end
-@propagate_inbounds function _copyto!(dest::StridedMatrix, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix})
+@propagate_inbounds function _copy!(dest::StridedMatrix, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix})
     U2 = Base.unalias(dest, U)
-    copyto_unaliased!(dest, U2)
+    copy_unaliased!(dest, U2)
     return dest
 end
 # for strided matrices, we explicitly loop over the arrays to improve cache locality
 # This fuses the copytrito! for the two halves
-@inline function copyto_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix})
+@inline function copy_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix})
     @boundscheck checkbounds(dest, axes(U)...)
     isunit = U isa UnitUpperTriangular
     for col in axes(dest,2)
@@ -624,7 +635,7 @@ end
     end
     return dest
 end
-@inline function copyto_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix})
+@inline function copy_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix})
     @boundscheck checkbounds(dest, axes(L)...)
     isunit = L isa UnitLowerTriangular
     for col in axes(dest,2)
@@ -645,7 +656,7 @@ Base.@constprop :aggressive function copytrito_triangular!(Bdata, Adata, uplo, u
         BLAS.chkuplo(uplo)
         LAPACK.lacpy_size_check(size(Bdata), sz)
         # only the diagonal is copied in this case
-        copyto!(diagview(Bdata), diagview(Adata))
+        copy!(diagview(Bdata), diagview(Adata))
     end
     return Bdata
 end
@@ -1055,15 +1066,17 @@ isunit_char(::UnitUpperTriangular) = 'U'
 isunit_char(::LowerTriangular) = 'N'
 isunit_char(::UnitLowerTriangular) = 'U'
 
+_copy_or_copyto!(dest, src) = ndims(dest) == ndims(src) ? copy!(dest, src) : copyto!(dest, src)
+
 # generic fallback for AbstractTriangular matrices outside of the four subtypes provided here
 _trimul!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVector) =
-    lmul!(A, copyto!(C, B))
+    lmul!(A, _copy_or_copyto!(C, B))
 _trimul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractMatrix) =
-    lmul!(A, copyto!(C, B))
+    lmul!(A, copy!(C, B))
 _trimul!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) =
-    rmul!(copyto!(C, A), B)
+    rmul!(copy!(C, A), B)
 _trimul!(C::AbstractMatrix, A::AbstractTriangular, B::AbstractTriangular) =
-    lmul!(A, copyto!(C, B))
+    lmul!(A, copy!(C, B))
 # redirect for UpperOrLowerTriangular
 _trimul!(C::AbstractVecOrMat, A::UpperOrLowerTriangular, B::AbstractVector) =
     generic_trimatmul!(C, uplo_char(A), isunit_char(A), wrapperop(parent(A)), _unwrap_at(parent(A)), B)
@@ -1124,9 +1137,9 @@ end
 ldiv!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = _ldiv!(C, A, B)
 # generic fallback for AbstractTriangular, directs to 2-arg [l/r]div!
 _ldiv!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) =
-    ldiv!(A, copyto!(C, B))
+    ldiv!(A, _copy_or_copyto!(C, B))
 _rdiv!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractTriangular) =
-    rdiv!(copyto!(C, A), B)
+    rdiv!(copy!(C, A), B)
 # redirect for UpperOrLowerTriangular to generic_*div!
 _ldiv!(C::AbstractVecOrMat, A::UpperOrLowerTriangular, B::AbstractVecOrMat) =
     generic_trimatdiv!(C, uplo_char(A), isunit_char(A), wrapperop(parent(A)), _unwrap_at(parent(A)), B)
@@ -1204,40 +1217,40 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
             elseif p == Inf
                 return inv(LAPACK.trcon!('I', $uploc, $isunitc, A.data))
             else # use fallback
-                return cond(copyto!(similar(parent(A)), A), p)
+                return cond(copy!(similar(parent(A)), A), p)
             end
         end
     end
 end
 
 # multiplication
 generic_trimatmul!(c::StridedVector{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, b::AbstractVector{T}) where {T<:BlasFloat} =
-    BLAS.trmv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, c === b ? c : copyto!(c, b))
+    BLAS.trmv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, c === b ? c : copy!(c, b))
 function generic_trimatmul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractMatrix{T}) where {T<:BlasFloat}
     if stride(C,1) == stride(A,1) == 1
-        BLAS.trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C === B ? C : copyto!(C, B))
+        BLAS.trmm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C === B ? C : copy!(C, B))
     else # incompatible with BLAS
         @invoke generic_trimatmul!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix)
     end
 end
 function generic_mattrimul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::AbstractMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat}
     if stride(C,1) == stride(B,1) == 1
-        BLAS.trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C === A ? C : copyto!(C, A))
+        BLAS.trmm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C === A ? C : copy!(C, A))
     else # incompatible with BLAS
         @invoke generic_mattrimul!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix)
     end
 end
 # division
 function generic_trimatdiv!(C::StridedVecOrMat{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractVecOrMat{T}) where {T<:BlasFloat}
     if stride(C,1) == stride(A,1) == 1
-        LAPACK.trtrs!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : copyto!(C, B))
+        LAPACK.trtrs!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : _copy_or_copyto!(C, B))
     else # incompatible with LAPACK
         @invoke generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractVecOrMat)
     end
 end
 function generic_mattridiv!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::AbstractMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat}
     if stride(C,1) == stride(B,1) == 1
-        BLAS.trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C === A ? C : copyto!(C, A))
+        BLAS.trsm!('R', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), B, C === A ? C : copy!(C, A))
     else # incompatible with BLAS
         @invoke generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix)
     end
@@ -1962,22 +1975,22 @@ function powm!(A0::UpperTriangular, p::Real)
         for i in axes(S,1)
             @inbounds S[i, i] = S[i, i] + 1
         end
-        copyto!(Stmp, S)
+        copy!(Stmp, S)
         mul!(S, A, c)
         ldiv!(Stmp, S)
 
         c = (p - j) / (j4 - 2)
         for i in axes(S,1)
             @inbounds S[i, i] = S[i, i] + 1
         end
-        copyto!(Stmp, S)
+        copy!(Stmp, S)
         mul!(S, A, c)
         ldiv!(Stmp, S)
     end
     for i in axes(S,1)
         S[i, i] = S[i, i] + 1
     end
-    copyto!(Stmp, S)
+    copy!(Stmp, S)
     mul!(S, A, -p)
     ldiv!(Stmp, S)
     for i in axes(S,1)
@@ -1987,7 +2000,7 @@ function powm!(A0::UpperTriangular, p::Real)
     blockpower!(A0, S, p/(2^s))
     for m = 1:s
         mul!(Stmp.data, S, S)
-        copyto!(S, Stmp)
+        copy!(S, Stmp)
         blockpower!(A0, S, p/(2^(s-m)))
     end
     rmul!(S, normA0^p)
@@ -2174,7 +2187,7 @@ function _find_params_log_quasitriu!(A)
             break
         end
         _sqrt_quasitriu!(A isa UpperTriangular ? parent(A) : A, A)
-        copyto!(AmI, A)
+        copy!(AmI, A)
         for i in axes(AmI,1)
             @inbounds AmI[i,i] -= 1
         end

test/bidiag.jl

@@ -1192,4 +1192,13 @@ end
     @test_throws msg ldiv!(C, B, zeros(2,1))
 end
 
+@testset "l/rmul with 0-sized matrices" begin
+    n = 0
+    B = Bidiagonal(ones(n), ones(max(n-1,0)), :U)
+    B2 = copy(B)
+    D = Diagonal(ones(n))
+    @test lmul!(D, B) == B2
+    @test rmul!(B, D) == B2
+end
+
 end # module TestBidiagonal