Merge branch 'master' into jishnub/issymmetric_number

Author: dkarrasch

.github/workflows/ci.yml

@@ -15,7 +15,7 @@ jobs:
   docs:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: 'nightly'

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}
@@ -325,7 +326,7 @@ StridedMatrixStride1{T} = StridedArrayStride1{T,2}
 """
     LinearAlgebra.checksquare(A)
 
-Check that a matrix is square, then return its common dimension.
+Checks whether a matrix is square, returning its common dimension if it is the case, or throwing a DimensionMismatch error otherwise.
 For multiple arguments, return a vector.
 
 # Examples
@@ -339,19 +340,13 @@ julia> LinearAlgebra.checksquare(A, B)
 ```
 """
 function checksquare(A)
-    m,n = size(A)
-    m == n || throw(DimensionMismatch(lazy"matrix is not square: dimensions are $(size(A))"))
-    m
+    sizeA = size(A)
+    length(sizeA) == 2 || throw(DimensionMismatch(lazy"input is not a matrix: dimensions are $sizeA"))
+    sizeA[1] == sizeA[2] || throw(DimensionMismatch(lazy"matrix is not square: dimensions are $sizeA"))
+    return sizeA[1]
 end
 
-function checksquare(A...)
-    sizes = Int[]
-    for a in A
-        size(a,1)==size(a,2) || throw(DimensionMismatch(lazy"matrix is not square: dimensions are $(size(a))"))
-        push!(sizes, size(a,1))
-    end
-    return sizes
-end
+checksquare(A...) = [checksquare(a) for a in A]
 
 function char_uplo(uplo::Symbol)
     if uplo === :U
@@ -363,7 +358,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 +374,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/abstractq.jl

@@ -42,10 +42,13 @@ convert(::Type{AbstractQ{T}}, adjQ::AdjointQ) where {T} = convert(AbstractQ{T},
 
 # ... to matrix
 Matrix{T}(Q::AbstractQ) where {T} = convert(Matrix{T}, Q*I) # generic fallback, yields square matrix
-Matrix{T}(adjQ::AdjointQ{S}) where {T,S} = convert(Matrix{T}, lmul!(adjQ, Matrix{S}(I, size(adjQ))))
 Matrix(Q::AbstractQ{T}) where {T} = Matrix{T}(Q)
 Array{T}(Q::AbstractQ) where {T} = Matrix{T}(Q)
 Array(Q::AbstractQ) = Matrix(Q)
+AbstractMatrix(Q::AbstractQ) = Q*I
+AbstractArray(Q::AbstractQ) = AbstractMatrix(Q)
+AbstractMatrix{T}(Q::AbstractQ) where {T} = Matrix{T}(Q)
+AbstractArray{T}(Q::AbstractQ) where {T} = AbstractMatrix{T}(Q)
 convert(::Type{T}, Q::AbstractQ) where {T<:AbstractArray} = T(Q)
 # legacy
 @deprecate(convert(::Type{AbstractMatrix{T}}, Q::AbstractQ) where {T},
@@ -172,10 +175,10 @@ qsize_check(Q::AbstractQ, P::AbstractQ) =
 # mimic the AbstractArray fallback
 *(Q::AbstractQ{<:Number}) = Q
 
-(*)(Q::AbstractQ, J::UniformScaling) = Q*J.λ
-function (*)(Q::AbstractQ, b::Number)
-    T = promote_type(eltype(Q), typeof(b))
-    lmul!(convert(AbstractQ{T}, Q), Matrix{T}(b*I, size(Q)))
+(*)(Q::AbstractQ, b::Number) = Q*(b*I)
+function (*)(Q::AbstractQ, J::UniformScaling)
+    T = promote_type(eltype(Q), eltype(J))
+    lmul!(convert(AbstractQ{T}, Q), Matrix{T}(J, size(Q)))
 end
 function (*)(Q::AbstractQ, B::AbstractVector)
     T = promote_type(eltype(Q), eltype(B))
@@ -188,10 +191,10 @@ function (*)(Q::AbstractQ, B::AbstractMatrix)
     mul!(similar(B, T, (size(Q, 1), size(B, 2))), convert(AbstractQ{T}, Q), B)
 end
 
-(*)(J::UniformScaling, Q::AbstractQ) = J.λ*Q
-function (*)(a::Number, Q::AbstractQ)
-    T = promote_type(typeof(a), eltype(Q))
-    rmul!(Matrix{T}(a*I, size(Q)), convert(AbstractQ{T}, Q))
+(*)(a::Number, Q::AbstractQ) = (a*I)*Q
+function (*)(J::UniformScaling, Q::AbstractQ)
+    T = promote_type(eltype(J), eltype(Q))
+    rmul!(Matrix{T}(J, size(Q)), convert(AbstractQ{T}, Q))
 end
 function (*)(A::AbstractVector, Q::AbstractQ)
     T = promote_type(eltype(A), eltype(Q))

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
@@ -1285,7 +1292,7 @@ function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
     nx, ny = length(x), length(y)
     (nx == size(B, 1) == ny) || throw(DimensionMismatch())
     if nx ≤ 1
-        nx == 0 && return dot(zero(eltype(x)), zero(eltype(B)), zero(eltype(y)))
+        nx == 0 && return zero(dot(zero(eltype(x)), zero(eltype(B)), zero(eltype(y))))
         return dot(x[1], B.dv[1], y[1])
     end
     ev, dv = B.ev, B.dv

src/blas.jl

@@ -2102,7 +2102,7 @@ end
     her2k!(uplo, trans, alpha, A, B, beta, C)
 
 Rank-2k update of the Hermitian matrix `C` as
-`alpha*A*B' + alpha*B*A' + beta*C` or `alpha*A'*B + alpha*B'*A + beta*C`
+`alpha*A*B' + alpha'*B*A' + beta*C` or `alpha*A'*B + alpha'*B'*A + beta*C`
 according to [`trans`](@ref stdlib-blas-trans). The scalar `beta` has to be real.
 Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `C` is used. Return `C`.
 """
@@ -2111,8 +2111,8 @@ function her2k! end
 """
     her2k(uplo, trans, alpha, A, B)
 
-Return the [`uplo`](@ref stdlib-blas-uplo) triangle of `alpha*A*B' + alpha*B*A'`
-or `alpha*A'*B + alpha*B'*A`, according to [`trans`](@ref stdlib-blas-trans).
+Return the [`uplo`](@ref stdlib-blas-uplo) triangle of `alpha*A*B' + alpha'*B*A'`
+or `alpha*A'*B + alpha'*B'*A`, according to [`trans`](@ref stdlib-blas-trans).
 """
 her2k(uplo, trans, alpha, A, B)
 

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
@@ -215,6 +217,12 @@ BunchKaufman{T}(B::BunchKaufman) where {T} =
     BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)
 Factorization{T}(B::BunchKaufman) where {T} = BunchKaufman{T}(B)
 
+AbstractMatrix(B::BunchKaufman) = B.uplo == 'U' ? B.P'B.U*B.D*B.U'B.P : B.P'B.L*B.D*B.L'B.P
+AbstractArray(B::BunchKaufman) = AbstractMatrix(B)
+Matrix(B::BunchKaufman) = convert(Array, AbstractArray(B))
+Array(B::BunchKaufman) = Matrix(B)
+
+
 size(B::BunchKaufman) = size(getfield(B, :LD))
 size(B::BunchKaufman, d::Integer) = size(getfield(B, :LD), d)
 issymmetric(B::BunchKaufman) = B.symmetric

src/cholesky.jl

@@ -646,7 +646,7 @@ Factorization{T}(C::CholeskyPivoted) where {T} = CholeskyPivoted{T}(C)
 
 AbstractMatrix(C::Cholesky) = C.uplo == 'U' ? C.U'C.U : C.L*C.L'
 AbstractArray(C::Cholesky) = AbstractMatrix(C)
-Matrix(C::Cholesky) = Array(AbstractArray(C))
+Matrix(C::Cholesky) = convert(Array, AbstractArray(C))
 Array(C::Cholesky) = Matrix(C)
 
 function AbstractMatrix(F::CholeskyPivoted)
@@ -655,7 +655,7 @@ function AbstractMatrix(F::CholeskyPivoted)
     U'U
 end
 AbstractArray(F::CholeskyPivoted) = AbstractMatrix(F)
-Matrix(F::CholeskyPivoted) = Array(AbstractArray(F))
+Matrix(F::CholeskyPivoted) = convert(Array, AbstractArray(F))
 Array(F::CholeskyPivoted) = Matrix(F)
 
 copy(C::Cholesky) = Cholesky(copy(C.factors), C.uplo, C.info)

src/dense.jl

@@ -237,22 +237,14 @@ end
 
 fillstored!(A::AbstractMatrix, v) = fill!(A, v)
 
-diagind(m::Integer, n::Integer, k::Integer=0) = diagind(IndexLinear(), m, n, k)
-diagind(::IndexLinear, m::Integer, n::Integer, k::Integer=0) =
-    k <= 0 ? range(1-k, step=m+1, length=min(m+k, n)) : range(k*m+1, step=m+1, length=min(m, n-k))
-
-function diagind(::IndexCartesian, m::Integer, n::Integer, k::Integer=0)
-    Cstart = CartesianIndex(1 + max(0,-k), 1 + max(0,k))
-    Cstep = CartesianIndex(1, 1)
-    length = max(0, k <= 0 ? min(m+k, n) : min(m, n-k))
-    StepRangeLen(Cstart, Cstep, length)
-end
-
 """
     diagind(M::AbstractMatrix, k::Integer = 0, indstyle::IndexStyle = IndexLinear())
     diagind(M::AbstractMatrix, indstyle::IndexStyle = IndexLinear())
+    diagind(::IndexStyle, m::Integer, n::Integer, k::Integer = 0)
+    diagind(m::Integer, n::Integer, k::Integer = 0)
 
-An `AbstractRange` giving the indices of the `k`th diagonal of the matrix `M`.
+An `AbstractRange` giving the indices of the `k`th diagonal of a matrix,
+specified either by the matrix `M` itself or by its dimensions `m` and `n`.
 Optionally, an index style may be specified which determines the type of the range returned.
 If `indstyle isa IndexLinear` (default), this returns an `AbstractRange{Integer}`.
 On the other hand, if `indstyle isa IndexCartesian`, this returns an `AbstractRange{CartesianIndex{2}}`.
@@ -262,6 +254,7 @@ If `k` is not provided, it is assumed to be `0` (corresponding to the main diago
 See also: [`diag`](@ref), [`diagm`](@ref), [`Diagonal`](@ref).
 
 # Examples
+The matrix itself may be passed to `diagind`:
 ```jldoctest
 julia> A = [1 2 3; 4 5 6; 7 8 9]
 3×3 Matrix{Int64}:
@@ -276,6 +269,18 @@ julia> diagind(A, IndexCartesian())
 StepRangeLen(CartesianIndex(1, 1), CartesianIndex(1, 1), 3)
 ```
 
+Alternatively, dimensions `m` and `n` may be passed to get the diagonal of an `m×n` matrix:
+```jldoctest
+julia> m, n = 5, 7
+(5, 7)
+
+julia> diagind(m, n, 2)
+11:6:35
+
+julia> diagind(IndexCartesian(), m, n)
+StepRangeLen(CartesianIndex(1, 1), CartesianIndex(1, 1), 5)
+```
+
 !!! compat "Julia 1.11"
      Specifying an `IndexStyle` requires at least Julia 1.11.
 """
@@ -286,6 +291,17 @@ end
 
 diagind(A::AbstractMatrix, indexstyle::IndexStyle) = diagind(A, 0, indexstyle)
 
+function diagind(::IndexCartesian, m::Integer, n::Integer, k::Integer=0)
+    Cstart = CartesianIndex(1 + max(0,-k), 1 + max(0,k))
+    Cstep = CartesianIndex(1, 1)
+    length = max(0, k <= 0 ? min(m+k, n) : min(m, n-k))
+    StepRangeLen(Cstart, Cstep, length)
+end
+
+diagind(::IndexLinear, m::Integer, n::Integer, k::Integer=0) =
+    k <= 0 ? range(1-k, step=m+1, length=min(m+k, n)) : range(k*m+1, step=m+1, length=min(m, n-k))
+diagind(m::Integer, n::Integer, k::Integer=0) = diagind(IndexLinear(), m, n, k)
+
 """
     diag(M, k::Integer=0)
 

src/diagonal.jl

@@ -970,6 +970,10 @@ function inv(D::Diagonal{T}) where T
     Diagonal(Di)
 end
 
+# Ensure doubly wrapped matrices use efficient diagonal methods and return a Symmetric/Hermitian type
+inv(A::Symmetric{<:Number,<:Diagonal}) = Symmetric(inv(A.data), sym_uplo(A.uplo))
+inv(A::Hermitian{<:Number,<:Diagonal}) = Hermitian(inv(real(A.data)), sym_uplo(A.uplo))
+
 function pinv(D::Diagonal{T}) where T
     Di = similar(D.diag, typeof(inv(oneunit(T))))
     for i = 1:length(D.diag)
@@ -1111,6 +1115,23 @@ end
 /(u::AdjointAbsVec, D::Diagonal) = (D' \ u')'
 /(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u))
 
+# norm
+function generic_normMinusInf(D::Diagonal)
+    norm_diag = norm(D.diag, -Inf)
+    return size(D,1) > 1 ? min(norm_diag, zero(norm_diag)) : norm_diag
+end
+generic_normInf(D::Diagonal) = norm(D.diag, Inf)
+generic_norm1(D::Diagonal) = norm(D.diag, 1)
+generic_norm2(D::Diagonal) = norm(D.diag)
+function generic_normp(D::Diagonal, p)
+    v = norm(D.diag, p)
+    if size(D,1) > 1 && p < 0
+        v = norm(zero(v), p)
+    end
+    return v
+end
+norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag)
+
 _opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag)
 _opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag)
 _opnorm12Inf(A::Diagonal, p) = maximum(opnorm(x, p) for x in A.diag)