Merge branch 'master' into dk/diag_svd_eigen

Author: dkarrasch

.ci/Manifest.toml

@@ -52,7 +52,7 @@ version = "1.11.0"
 
 [[deps.JuliaSyntaxHighlighting]]
 deps = ["StyledStrings"]
-uuid = "dc6e5ff7-fb65-4e79-a425-ec3bc9c03011"
+uuid = "ac6e5ff7-fb65-4e79-a425-ec3bc9c03011"
 version = "1.12.0"
 
 [[deps.LazyArtifacts]]
@@ -93,7 +93,7 @@ version = "1.11.0"
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
 path = ".."
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-version = "1.11.0"
+version = "1.12.0"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

src/LinearAlgebra.jl

@@ -16,6 +16,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     permutedims, permuterows!, power_by_squaring, promote_rule, real, sec, sech, setindex!,
     show, similar, sin, sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc,
     typed_hcat, vec, view, zero
+import Base: AbstractArray, AbstractMatrix, Array, Matrix
 using Base: IndexLinear, promote_eltype, promote_op, print_matrix,
     @propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing,
     splat, BitInteger

src/blas.jl

@@ -84,7 +84,8 @@ export
     trsm!,
     trsm
 
-using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, chkstride1
+using ..LinearAlgebra: libblastrampoline, BlasReal, BlasComplex, BlasFloat, BlasInt,
+    DimensionMismatch, checksquare, chkstride1, SingularException
 
 include("lbt.jl")
 
@@ -128,6 +129,10 @@ Set the number of threads the BLAS library should use equal to `n::Integer`.
 
 Also accepts `nothing`, in which case julia tries to guess the default number of threads.
 Passing `nothing` is discouraged and mainly exists for historical reasons.
+
+!!! note
+    Some BLAS libraries, such as Apple Accelerate, cannot be configured to use a fixed number of threads.
+    For these backends, `set_num_threads()` is a no-op. See also [`get_num_threads`](@ref).
 """
 set_num_threads(nt::Integer)::Nothing = lbt_set_num_threads(Int32(nt))
 function set_num_threads(::Nothing)
@@ -147,6 +152,10 @@ Get the number of threads the BLAS library is using.
 
 !!! compat "Julia 1.6"
     `get_num_threads` requires at least Julia 1.6.
+
+!!! note
+    Some BLAS libraries, such as Apple Accelerate, cannot be configured to use a fixed number of threads.
+    For these backends, `get_num_threads()` always returns `1`. See also [`set_num_threads`](@ref).
 """
 get_num_threads()::Int = lbt_get_num_threads()
 
@@ -1369,6 +1378,11 @@ for (fname, elty) in ((:dtrsv_,:Float64),
                 throw(DimensionMismatch(lazy"size of A is $n != length(x) = $(length(x))"))
             end
             chkstride1(A)
+            if diag == 'N'
+                for i in 1:n
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             px, stx = vec_pointer_stride(x, ArgumentError("input vector with 0 stride is not allowed"))
             GC.@preserve x ccall((@blasfunc($fname), libblastrampoline), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
@@ -2217,6 +2231,11 @@ for (mmname, smname, elty) in
             end
             chkstride1(A)
             chkstride1(B)
+            if diag == 'N'
+                for i in 1:k
+                    iszero(A[i,i]) && throw(SingularException(i))
+                end
+            end
             ccall((@blasfunc($smname), libblastrampoline), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                     Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},

src/cholesky.jl

@@ -65,7 +65,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -305,7 +305,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{UpperTriangular}, tol::Real, ch
     rTA = real(eltype(A))
     # checks
     Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
+    n = checksquare(A)
     # initialization
     piv = collect(1:n)
     dots = zeros(rTA, n)
@@ -354,7 +354,7 @@ function _cholpivoted!(A::AbstractMatrix, ::Type{LowerTriangular}, tol::Real, ch
     rTA = real(eltype(A))
     # checks
     Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
+    n = checksquare(A)
     # initialization
     piv = collect(1:n)
     dots = zeros(rTA, n)
@@ -530,7 +530,7 @@ julia> C.U
   ⋅    ⋅    3.0
 
 julia> C.L
-3×3 LowerTriangular{Float64, Matrix{Float64}}:
+3×3 LowerTriangular{Float64, Adjoint{Matrix{Float64}}}:
   2.0   ⋅    ⋅
   6.0  1.0   ⋅
  -8.0  5.0  3.0
@@ -664,30 +664,15 @@ copy(C::CholeskyPivoted) = CholeskyPivoted(copy(C.factors), C.uplo, C.piv, C.ran
 size(C::Union{Cholesky, CholeskyPivoted}) = size(C.factors)
 size(C::Union{Cholesky, CholeskyPivoted}, d::Integer) = size(C.factors, d)
 
-function _choleskyUfactor(Cfactors, Cuplo)
-    if Cuplo === 'U'
-        return UpperTriangular(Cfactors)
-    else
-        return copy(LowerTriangular(Cfactors)')
-    end
-end
-function _choleskyLfactor(Cfactors, Cuplo)
-    if Cuplo === 'L'
-        return LowerTriangular(Cfactors)
-    else
-        return copy(UpperTriangular(Cfactors)')
-    end
-end
-
 function getproperty(C::Cholesky, d::Symbol)
     Cfactors = getfield(C, :factors)
     Cuplo    = getfield(C, :uplo)
     if d === :U
-        _choleskyUfactor(Cfactors, Cuplo)
+        UpperTriangular(Cuplo == 'U' ? Cfactors : Cfactors')
     elseif d === :L
-        _choleskyLfactor(Cfactors, Cuplo)
+        LowerTriangular(Cuplo == 'L' ? Cfactors : Cfactors')
     elseif d === :UL
-        return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
+        return (Cuplo == 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors))
     else
         return getfield(C, d)
     end
@@ -704,9 +689,9 @@ function getproperty(C::CholeskyPivoted{T}, d::Symbol) where {T}
     Cfactors = getfield(C, :factors)
     Cuplo    = getfield(C, :uplo)
     if d === :U
-        _choleskyUfactor(Cfactors, Cuplo)
+        UpperTriangular(Cuplo == 'U' ? Cfactors : Cfactors')
     elseif d === :L
-        _choleskyLfactor(Cfactors, Cuplo)
+        LowerTriangular(Cuplo == 'L' ? Cfactors : Cfactors')
     elseif d === :p
         return getfield(C, :piv)
     elseif d === :P
@@ -813,7 +798,7 @@ function rdiv!(B::AbstractMatrix, C::Cholesky)
     end
 end
 
-function LinearAlgebra.rdiv!(B::AbstractMatrix, C::CholeskyPivoted)
+function rdiv!(B::AbstractMatrix, C::CholeskyPivoted)
     n = size(C, 2)
     for i in 1:size(B, 1)
         permute!(view(B, i, 1:n), C.piv)
@@ -965,7 +950,7 @@ function lowrankdowndate!(C::Cholesky, v::AbstractVector)
         s = conj(v[i]/Aii)
         s2 = abs2(s)
         if s2 > 1
-            throw(LinearAlgebra.PosDefException(i))
+            throw(PosDefException(i))
         end
         c = sqrt(1 - abs2(s))
 

src/dense.jl

@@ -972,7 +972,8 @@ sqrt(::AbstractMatrix)
 function sqrt(A::AbstractMatrix{T}) where {T<:Union{Real,Complex}}
     if checksquare(A) == 0
         return copy(float(A))
-    elseif isdiag(A)
+    elseif isdiag(A) && (T <: Complex || all(x -> x ≥ zero(x), diagview(A)))
+        # Real Diagonal sqrt requires each diagonal element to be positive
         return applydiagonal(sqrt, A)
     elseif ishermitian(A)
         sqrtHermA = sqrt(Hermitian(A))
@@ -1587,24 +1588,25 @@ factorize(A::Transpose) = transpose(factorize(parent(A)))
 factorize(a::Number)    = a # same as how factorize behaves on Diagonal types
 
 function getstructure(A::StridedMatrix)
+    require_one_based_indexing(A)
     m, n = size(A)
     if m == 1 return A[1] end
     utri    = true
     utri1   = true
     herm    = true
     sym     = true
-    for j = 1:n-1, i = j+1:m
-        if utri1
+    for j = 1:n, i = j:m
+        if (j < n) && (i > j) && utri1 # indices are off-diagonal
             if A[i,j] != 0
                 utri1 = i == j + 1
                 utri = false
             end
         end
         if sym
-            sym &= A[i,j] == A[j,i]
+            sym &= A[i,j] == transpose(A[j,i])
         end
         if herm
-            herm &= A[i,j] == conj(A[j,i])
+            herm &= A[i,j] == adjoint(A[j,i])
         end
         if !(utri1|herm|sym) break end
     end
@@ -1617,10 +1619,12 @@ function getstructure(A::StridedMatrix)
     if ltri1
         for i = 1:n-1
             if A[i,i+1] != 0
-                ltri &= false
+                ltri = false
                 break
             end
         end
+    else
+        ltri = false
     end
     return (utri, utri1, ltri, ltri1, sym, herm)
 end
@@ -1779,6 +1783,10 @@ Condition number of the matrix `M`, computed using the operator `p`-norm. Valid
 """
 function cond(A::AbstractMatrix, p::Real=2)
     if p == 2
+        if isempty(A)
+            checksquare(A)
+            return zero(real(eigtype(eltype(A))))
+        end
         v = svdvals(A)
         maxv = maximum(v)
         return iszero(maxv) ? oftype(real(maxv), Inf) : maxv / minimum(v)

src/exceptions.jl

@@ -28,6 +28,10 @@ struct SingularException <: Exception
     info::BlasInt
 end
 
+function Base.showerror(io::IO, ex::SingularException)
+    print(io, "SingularException: matrix is singular; factorization failed. Zero pivot found at index ", ex.info)
+end
+
 """
     PosDefException
 

src/lapack.jl

@@ -2975,10 +2975,10 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in
             mA  = size(A, 1)
             k   = length(tau)
             if side == 'L' && m != mA
-                throw(DimensionMismatch(lazy"for a left-sided multiplication, the first dimension of C, $m, must equal the second dimension of A, $mA"))
+                throw(DimensionMismatch(lazy"for a left-sided multiplication, the first dimension of C, $m, must equal the first dimension of A, $mA"))
             end
             if side == 'R' && n != mA
-                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $m, must equal the second dimension of A, $mA"))
+                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $n, must equal the first dimension of A, $mA"))
             end
             if side == 'L' && k > m
                 throw(DimensionMismatch(lazy"invalid number of reflectors: k = $k should be <= m = $m"))
@@ -3025,10 +3025,10 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in
             mA  = size(A, 1)
             k   = length(tau)
             if side == 'L' && m != mA
-                throw(DimensionMismatch(lazy"for a left-sided multiplication, the first dimension of C, $m, must equal the second dimension of A, $mA"))
+                throw(DimensionMismatch(lazy"for a left-sided multiplication, the first dimension of C, $m, must equal the first dimension of A, $mA"))
             end
             if side == 'R' && n != mA
-                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $m, must equal the second dimension of A, $mA"))
+                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $n, must equal the first dimension of A, $mA"))
             end
             if side == 'L' && k > m
                 throw(DimensionMismatch(lazy"invalid number of reflectors: k = $k should be <= m = $m"))
@@ -3078,7 +3078,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in
                 throw(DimensionMismatch(lazy"for a left-sided multiplication, the first dimension of C, $m, must equal the second dimension of A, $nA"))
             end
             if side == 'R' && n != nA
-                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $m, must equal the second dimension of A, $nA"))
+                throw(DimensionMismatch(lazy"for a right-sided multiplication, the second dimension of C, $n, must equal the second dimension of A, $nA"))
             end
             if side == 'L' && k > m
                 throw(DimensionMismatch(lazy"invalid number of reflectors: k = $k should be <= m = $m"))
@@ -5391,9 +5391,9 @@ for (syev, syevr, syevd, sygvd, elty) in
             W = similar(A, $elty, n)
             ldz = n
             if jobz == 'N'
-                Z = similar(A, $elty, ldz, 0)
+                Z = similar(A, $elty, 0)
             elseif jobz == 'V'
-                Z = similar(A, $elty, ldz, n)
+                Z = similar(A, $elty, ldz * n)
             end
             isuppz = similar(A, BlasInt, 2*n)
             work   = Vector{$elty}(undef, 1)
@@ -5423,7 +5423,8 @@ for (syev, syevr, syevd, sygvd, elty) in
                     resize!(iwork, liwork)
                 end
             end
-            W[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
+            zm = jobz == 'V' ? m[] : 0
+            resize!(W, m[]), reshape(resize!(Z, ldz * zm), ldz, zm)
         end
         syevr!(jobz::AbstractChar, A::AbstractMatrix{$elty}) =
             syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0)
@@ -5593,10 +5594,10 @@ for (syev, syevr, syevd, sygvd, elty, relty) in
             W = similar(A, $relty, n)
             if jobz == 'N'
                 ldz = 1
-                Z = similar(A, $elty, ldz, 0)
+                Z = similar(A, $elty, 0)
             elseif jobz == 'V'
                 ldz = n
-                Z = similar(A, $elty, ldz, n)
+                Z = similar(A, $elty, ldz * n)
             end
             isuppz = similar(A, BlasInt, 2*n)
             work   = Vector{$elty}(undef, 1)
@@ -5632,7 +5633,8 @@ for (syev, syevr, syevd, sygvd, elty, relty) in
                     resize!(iwork, liwork)
                 end
             end
-            W[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)]
+            zm = jobz == 'V' ? m[] : 0
+            resize!(W, m[]), reshape(resize!(Z, ldz * zm), ldz, zm)
         end
         syevr!(jobz::AbstractChar, A::AbstractMatrix{$elty}) =
             syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0)

src/matmul.jl

@@ -1021,6 +1021,7 @@ function _generic_matmatmul_nonadjtrans!(C, A, B, alpha, beta)
     @inbounds for n in axes(B, 2), k in axes(B, 1)
         # Balpha = B[k,n] * alpha, but we skip the multiplication in case isone(alpha)
         Balpha = @stable_muladdmul MulAddMul(alpha, false)(B[k,n])
+        !ismissing(Balpha) && iszero(Balpha) && continue
         @simd for m in axes(A, 1)
             C[m,n] = muladd(A[m,k], Balpha, C[m,n])
         end

src/symmetric.jl

@@ -703,6 +703,19 @@ for f in (:+, :-)
 end
 
 *(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
+# catch a few potential BLAS-cases
+function *(A::HermOrSym{<:BlasFloat,<:StridedMatrix}, B::AdjOrTrans{<:BlasFloat,<:StridedMatrix})
+    T = promote_type(eltype(A), eltype(B))
+    mul!(similar(B, T, (size(A, 1), size(B, 2))),
+            convert(AbstractMatrix{T}, A),
+            copy_oftype(B, T)) # make sure the AdjOrTrans wrapper is resolved
+end
+function *(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloat,<:StridedMatrix})
+    T = promote_type(eltype(A), eltype(B))
+    mul!(similar(B, T, (size(A, 1), size(B, 2))),
+            copy_oftype(A, T), # make sure the AdjOrTrans wrapper is resolved
+            convert(AbstractMatrix{T}, B))
+end
 
 function dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
     require_one_based_indexing(x, y)

src/triangular.jl

@@ -1223,11 +1223,13 @@ function generic_mattrimul!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function,
     end
 end
 # division
-function generic_trimatdiv!(C::StridedVecOrMat{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractVecOrMat{T}) where {T<:BlasFloat}
+generic_trimatdiv!(C::StridedVector{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractVector{T}) where {T<:BlasFloat} =
+    BLAS.trsv!(uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, A, C === B ? C : copyto!(C, B))
+function generic_trimatdiv!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::StridedMatrix{T}, B::AbstractMatrix{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))
+        BLAS.trsm!('L', uploc, tfun === identity ? 'N' : tfun === transpose ? 'T' : 'C', isunitc, one(T), A, C === B ? C : copyto!(C, B))
     else # incompatible with LAPACK
-        @invoke generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractVecOrMat)
+        @invoke generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix)
     end
 end
 function generic_mattridiv!(C::StridedMatrix{T}, uploc, isunitc, tfun::Function, A::AbstractMatrix{T}, B::StridedMatrix{T}) where {T<:BlasFloat}