Merge branch 'master' of https://github.com/JuliaLang/LinearAlgebra.jl

Author: GianmarcoCuppari

src/LinearAlgebra.jl

@@ -326,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
@@ -340,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

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/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

@@ -217,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

@@ -1111,6 +1111,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)

src/eigen.jl

@@ -276,6 +276,8 @@ eigvecs(A::Union{Number, AbstractMatrix}; kws...) =
 eigvecs(F::Union{Eigen, GeneralizedEigen}) = F.vectors
 
 eigvals(F::Union{Eigen, GeneralizedEigen}) = F.values
+eigmin(F::Union{Eigen, GeneralizedEigen}) = minimum(eigvals(F))
+eigmax(F::Union{Eigen, GeneralizedEigen}) = maximum(eigvals(F))
 
 """
     eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values
@@ -357,13 +359,8 @@ eigvals(x::Number; kwargs...) = imag(x) == 0 ? real(x) : x
 """
     eigmax(A; permute::Bool=true, scale::Bool=true)
 
-Return the largest eigenvalue of `A`.
-The option `permute=true` permutes the matrix to become
-closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
-make rows and columns more equal in norm.
-Note that if the eigenvalues of `A` are complex,
-this method will fail, since complex numbers cannot
-be sorted.
+Return the largest eigenvalue of `A`, assuming they are all real, and throwing an error otherwise.
+The `permute` and `scale` keyword arguments are the same as for [`eigen`](@ref).
 
 # Examples
 ```jldoctest
@@ -388,23 +385,18 @@ Stacktrace:
 ```
 """
 function eigmax(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true)
-    v = eigvals(A, permute = permute, scale = scale)
+    v = eigvals(A; permute, scale)
     if eltype(v)<:Complex
         throw(DomainError(A, "`A` cannot have complex eigenvalues."))
     end
-    maximum(v)
+    return maximum(v)
 end
 
 """
     eigmin(A; permute::Bool=true, scale::Bool=true)
 
-Return the smallest eigenvalue of `A`.
-The option `permute=true` permutes the matrix to become
-closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
-make rows and columns more equal in norm.
-Note that if the eigenvalues of `A` are complex,
-this method will fail, since complex numbers cannot
-be sorted.
+Return the smallest eigenvalue of `A`, assuming they are all real, and throwing an error otherwise.
+The `permute` and `scale` keyword arguments are the same as for [`eigen`](@ref).
 
 # Examples
 ```jldoctest
@@ -428,13 +420,12 @@ Stacktrace:
 [...]
 ```
 """
-function eigmin(A::Union{Number, AbstractMatrix};
-                permute::Bool=true, scale::Bool=true)
-    v = eigvals(A, permute = permute, scale = scale)
+function eigmin(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true)
+    v = eigvals(A; permute, scale)
     if eltype(v)<:Complex
         throw(DomainError(A, "`A` cannot have complex eigenvalues."))
     end
-    minimum(v)
+    return minimum(v)
 end
 
 inv(A::Eigen) = A.vectors * inv(Diagonal(A.values)) / A.vectors

test/abstractq.jl

@@ -68,6 +68,7 @@ n = 5
         end
         @test convert(Matrix, Q) ≈ Matrix(Q) ≈ Q[:,:] ≈ copyto!(zeros(T, size(Q)), Q) ≈ Q.Q*I
         @test convert(Matrix, Q') ≈ Matrix(Q') ≈ (Q')[:,:] ≈ copyto!(zeros(T, size(Q)), Q') ≈ Q.Q'*I
+        @test AbstractMatrix(Q) ≈ AbstractArray(Q) ≈ AbstractMatrix{T}(Q) ≈ AbstractArray{T}(Q)
         @test Q[1,:] == Q.Q[1,:] == view(Q, 1, :)
         @test Q[:,1] == Q.Q[:,1] == view(Q, :, 1)
         @test Q[1,1] == Q.Q[1,1]

test/bunchkaufman.jl

@@ -259,4 +259,13 @@ end
     @test B.U * B.D * B.U' ≈ S
 end
 
+@testset "BunchKaufman array constructors #1461" begin
+    a = randn(5,5)
+    A = a'a
+    for ul in (:U, :L)
+        B = bunchkaufman(Symmetric(A, ul))
+        @test A ≈ Array(B) ≈ Matrix(B) ≈ AbstractArray(B) ≈ AbstractMatrix(B)
+    end
+end
+
 end # module TestBunchKaufman