Merge remote-tracking branch 'upstream/master'

Author: GianmarcoCuppari

docs/src/index.md

@@ -604,6 +604,7 @@ LinearAlgebra.hermitianpart
 LinearAlgebra.hermitianpart!
 LinearAlgebra.copy_adjoint!
 LinearAlgebra.copy_transpose!
+LinearAlgebra.uplo
 ```
 
 ## Low-level matrix operations

src/LinearAlgebra.jl

@@ -181,7 +181,8 @@ public AbstractTriangular,
         zeroslike,
         matprod_dest,
         fillstored!,
-        fillband!
+        fillband!,
+        uplo
 
 const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
 const BlasReal = Union{Float64,Float32}
@@ -363,7 +364,14 @@ function char_uplo(uplo::Symbol)
     end
 end
 
+"""
+    sym_uplo(uplo::Char)
+
+Return the `Symbol` corresponding the `uplo` by checking for validity.
+"""
 function sym_uplo(uplo::Char)
+    # This method is called by other packages, and isn't used within LinearAlgebra
+    # It's retained here for backward compatibility.
     if uplo == 'U'
         return :U
     elseif uplo == 'L'
@@ -372,6 +380,13 @@ function sym_uplo(uplo::Char)
         throw_uplo()
     end
 end
+"""
+    _sym_uplo(uplo::Char)
+
+Return the `Symbol` corresponding to `uplo` without checking for validity.
+See also `sym_uplo`, which checks for validity.
+"""
+_sym_uplo(uplo::Char) = uplo == 'U' ? (:U) : (:L)
 
 @noinline throw_uplo() = throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
 

src/bidiag.jl

@@ -115,6 +115,13 @@ Bidiagonal(A::Bidiagonal) = A
 Bidiagonal{T}(A::Bidiagonal{T}) where {T} = A
 Bidiagonal{T}(A::Bidiagonal) where {T} = Bidiagonal{T}(A.dv, A.ev, A.uplo)
 
+"""
+    LinearAlgebra.uplo(S::Bidiagonal)::Symbol
+
+Return a `Symbol` corresponding to whether the upper (`:U`) or lower (`:L`) off-diagonal band is stored.
+"""
+uplo(B::Bidiagonal) = sym_uplo(B.uplo)
+
 _offdiagind(uplo) = uplo == 'U' ? 1 : -1
 
 @inline function Base.isassigned(A::Bidiagonal, i::Int, j::Int)
@@ -295,7 +302,7 @@ function show(io::IO, M::Bidiagonal)
     print(io, ", ")
     show(io, M.ev)
     print(io, ", ")
-    show(io, sym_uplo(M.uplo))
+    show(io, _sym_uplo(M.uplo))
     print(io, ")")
 end
 
@@ -910,12 +917,12 @@ function _bidimul!(C::AbstractMatrix, A::BiTriSym, B::Diagonal, _add::MulAddMul)
     @inbounds begin
         # first row of C
         for j in 1:min(2, n)
-            C[1,j] += _add(A[1,j]*B[j,j])
+            C[1,j] += _add(A[1,j]*Bd[j])
         end
         # second row of C
         if n > 1
             for j in 1:min(3, n)
-                C[2,j] += _add(A[2,j]*B[j,j])
+                C[2,j] += _add(A[2,j]*Bd[j])
             end
         end
         for j in 3:n-2
@@ -926,13 +933,13 @@ function _bidimul!(C::AbstractMatrix, A::BiTriSym, B::Diagonal, _add::MulAddMul)
         if n > 3
             # row before last of C
             for j in n-2:n
-                C[n-1,j] += _add(A[n-1,j]*B[j,j])
+                C[n-1,j] += _add(A[n-1,j]*Bd[j])
             end
         end
         # last row of C
         if n > 2
             for j in n-1:n
-                C[n,j] += _add(A[n,j]*B[j,j])
+                C[n,j] += _add(A[n,j]*Bd[j])
             end
         end
     end # inbounds

src/bunchkaufman.jl

@@ -130,7 +130,9 @@ function bunchkaufman!(A::StridedMatrix{<:BlasFloat}, rook::Bool = false; check:
 end
 
 bkcopy_oftype(A, S) = eigencopy_oftype(A, S)
-bkcopy_oftype(A::Symmetric{<:Complex}, S) = Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), sym_uplo(A.uplo))
+function bkcopy_oftype(A::Symmetric{<:Complex}, S)
+    Symmetric(copytrito!(similar(parent(A), S, size(A)), A.data, A.uplo), _sym_uplo(A.uplo))
+end
 
 """
     bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman

src/generic.jl

@@ -579,20 +579,14 @@ function generic_norm2(x)
     T = typeof(maxabs)
     if isfinite(length(x)*maxabs*maxabs) && !iszero(maxabs*maxabs) # Scaling not necessary
         sum::promote_type(Float64, T) = norm_sqr(v)
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             sum += norm_sqr(v)
         end
         ismissing(sum) && return missing
         return convert(T, sqrt(sum))
     else
         sum = abs2(norm(v)/maxabs)
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             sum += (norm(v)/maxabs)^2
         end
         ismissing(sum) && return missing
@@ -614,21 +608,15 @@ function generic_normp(x, p)
     spp::promote_type(Float64, T) = p
     if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && !iszero(maxabs^spp)) # scaling not necessary
         sum::promote_type(Float64, T) = norm(v)^spp
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             ismissing(v) && return missing
             sum += norm(v)^spp
         end
         return convert(T, sum^inv(spp))
     else # rescaling
         sum = (norm(v)/maxabs)^spp
         ismissing(sum) && return missing
-        while true
-            y = iterate(x, s)
-            y === nothing && break
-            (v, s) = y
+        for v in Iterators.rest(x, s)
             ismissing(v) && return missing
             sum += (norm(v)/maxabs)^spp
         end
@@ -996,7 +984,8 @@ function dot(x, y) # arbitrary iterables
     return s
 end
 
-dot(x::Number, y::Number) = conj(x) * y
+# the unary + is for type promotion in the Boolean case, mimicking the reduction in usual dot
+dot(x::Number, y::Number) = +(conj(x) * y)
 
 function dot(x::AbstractArray, y::AbstractArray)
     lx = length(x)

src/special.jl

@@ -292,19 +292,19 @@ end
 
 for f in (:+, :-)
     @eval function $f(D::Diagonal{<:Number}, S::Symmetric)
-        uplo = sym_uplo(S.uplo)
+        uplo = _sym_uplo(S.uplo)
         return Symmetric(parentof_applytri($f, Symmetric(D, uplo), S), uplo)
     end
     @eval function $f(S::Symmetric, D::Diagonal{<:Number})
-        uplo = sym_uplo(S.uplo)
+        uplo = _sym_uplo(S.uplo)
         return Symmetric(parentof_applytri($f, S, Symmetric(D, uplo)), uplo)
     end
     @eval function $f(D::Diagonal{<:Real}, H::Hermitian)
-        uplo = sym_uplo(H.uplo)
+        uplo = _sym_uplo(H.uplo)
         return Hermitian(parentof_applytri($f, Hermitian(D, uplo), H), uplo)
     end
     @eval function $f(H::Hermitian, D::Diagonal{<:Real})
-        uplo = sym_uplo(H.uplo)
+        uplo = _sym_uplo(H.uplo)
         return Hermitian(parentof_applytri($f, H, Hermitian(D, uplo)), uplo)
     end
 end
@@ -608,8 +608,8 @@ end
 # tridiagonal cholesky factorization
 function cholesky(S::RealSymHermitian{<:BiTriSym}, ::NoPivot = NoPivot(); check::Bool = true)
     T = choltype(S)
-    B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), sym_uplo(S.uplo))
-    cholesky!(Hermitian(B, sym_uplo(S.uplo)), NoPivot(); check = check)
+    B = Bidiagonal{T}(diag(S, 0), diag(S, S.uplo == 'U' ? 1 : -1), _sym_uplo(S.uplo))
+    cholesky!(Hermitian(B, _sym_uplo(S.uplo)), NoPivot(); check = check)
 end
 
 # istriu/istril for triangular wrappers of structured matrices

src/structuredbroadcast.jl

@@ -230,6 +230,30 @@ end
 # All structured matrices are square, and therefore they only broadcast out if they are size (1, 1)
 Broadcast.newindex(D::StructuredMatrix, I::CartesianIndex{2}) = size(D) == (1,1) ? CartesianIndex(1,1) : I
 
+# Recursively replace wrapped matrices by their parents to improve broadcasting performance
+# We may do this because the indexing within `copyto!` is restricted to the stored indices
+preprocess_broadcasted(::Type{T}, A) where {T} = _preprocess_broadcasted(T, A)
+function preprocess_broadcasted(::Type{T}, bc::Broadcasted) where {T}
+    args = map(x -> preprocess_broadcasted(T, x), bc.args)
+    Broadcast.broadcasted(bc.f, args...)
+end
+# fallback case that doesn't unwrap at all
+_preprocess_broadcasted(::Type, x) = x
+
+_preprocess_broadcasted(::Type{Diagonal}, d::Diagonal) = d.diag
+# fallback for types that might opt into Diagonal-like structured broadcasting, e.g. wrappers
+_preprocess_broadcasted(::Type{Diagonal}, d::AbstractMatrix) = diagview(d)
+
+function copy(bc::Broadcasted{StructuredMatrixStyle{Diagonal}})
+    if isstructurepreserving(bc) || fzeropreserving(bc)
+        # forward the broadcasting operation to the diagonal
+        bc2 = preprocess_broadcasted(Diagonal, bc)
+        return Diagonal(copy(bc2))
+    else
+        @invoke copy(bc::Broadcasted)
+    end
+end
+
 function copyto!(dest::Diagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
@@ -291,13 +315,6 @@ function copyto!(dest::Tridiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     return dest
 end
 
-# Recursively replace wrapped matrices by their parents to improve broadcasting performance
-# We may do this because the indexing within `copyto!` is restricted to the stored indices
-preprocess_broadcasted(::Type{T}, A) where {T} = _preprocess_broadcasted(T, A)
-function preprocess_broadcasted(::Type{T}, bc::Broadcasted) where {T}
-    args = map(x -> preprocess_broadcasted(T, x), bc.args)
-    Broadcast.Broadcasted(bc.f, args, bc.axes)
-end
 _preprocess_broadcasted(::Type{LowerTriangular}, A) = lowertridata(A)
 _preprocess_broadcasted(::Type{UpperTriangular}, A) = uppertridata(A)
 _preprocess_broadcasted(::Type{UpperHessenberg}, A) = upperhessenbergdata(A)