From f927d25f21fd336b1b459369d7f59f1c7cb8d4d6 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 7 Oct 2022 12:13:06 +0200 Subject: [PATCH] audit `copymutable_oftype` usage (#47063) --- stdlib/LinearAlgebra/src/cholesky.jl | 6 ++--- stdlib/LinearAlgebra/src/diagonal.jl | 31 +++++-------------------- stdlib/LinearAlgebra/src/hessenberg.jl | 2 +- stdlib/LinearAlgebra/src/lq.jl | 6 ++--- stdlib/LinearAlgebra/src/lu.jl | 6 ++--- stdlib/LinearAlgebra/src/special.jl | 4 ++-- stdlib/LinearAlgebra/src/svd.jl | 21 +++++++---------- stdlib/LinearAlgebra/src/transpose.jl | 9 +++++++ stdlib/LinearAlgebra/test/cholesky.jl | 6 ++--- stdlib/LinearAlgebra/test/hessenberg.jl | 7 ++++++ stdlib/LinearAlgebra/test/lq.jl | 10 +++++--- stdlib/LinearAlgebra/test/svd.jl | 12 ++++++++++ 12 files changed, 64 insertions(+), 56 deletions(-) diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl index 917c32625adb5..8e5c85ac88948 100644 --- a/stdlib/LinearAlgebra/src/cholesky.jl +++ b/stdlib/LinearAlgebra/src/cholesky.jl @@ -178,10 +178,8 @@ Base.iterate(C::CholeskyPivoted, ::Val{:done}) = nothing # make a copy that allow inplace Cholesky factorization -@inline choltype(A) = promote_type(typeof(sqrt(oneunit(eltype(A)))), Float32) -@inline cholcopy(A::StridedMatrix) = copymutable_oftype(A, choltype(A)) -@inline cholcopy(A::RealHermSymComplexHerm) = copymutable_oftype(A, choltype(A)) -@inline cholcopy(A::AbstractMatrix) = copy_similar(A, choltype(A)) +choltype(A) = promote_type(typeof(sqrt(oneunit(eltype(A)))), Float32) +cholcopy(A::AbstractMatrix) = eigencopy_oftype(A, choltype(A)) # _chol!. Internal methods for calling unpivoted Cholesky ## BLAS/LAPACK element types diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 32687404752ff..3713deb026040 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -252,38 +252,19 @@ end rmul!(A::AbstractMatrix, D::Diagonal) = @inline mul!(A, A, D) lmul!(D::Diagonal, B::AbstractVecOrMat) = @inline mul!(B, D, B) -#TODO: It seems better to call (D' * adjA')' directly? -function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal) - A = adjA.parent - Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) - adjoint!(Ac, A) +function *(A::AdjOrTransAbsMat, D::Diagonal) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))) rmul!(Ac, D) end -function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal) - A = transA.parent - At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) - transpose!(At, A) - rmul!(At, D) -end - *(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ) -function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix}) - A = adjA.parent - Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) - adjoint!(Ac, A) +function *(D::Diagonal, A::AdjOrTransAbsMat) + Ac = copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))) lmul!(D, Ac) end -function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix}) - A = transA.parent - At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1))) - transpose!(At, A) - lmul!(D, At) -end - @inline function __muldiag!(out, D::Diagonal, B, alpha, beta) require_one_based_indexing(B) require_one_based_indexing(out) @@ -853,8 +834,8 @@ end inv(C::Cholesky{<:Any,<:Diagonal}) = Diagonal(map(inv∘abs2, C.factors.diag)) -@inline cholcopy(A::Diagonal) = copymutable_oftype(A, choltype(A)) -@inline cholcopy(A::RealHermSymComplexHerm{<:Real,<:Diagonal}) = copymutable_oftype(A, choltype(A)) +cholcopy(A::Diagonal) = copymutable_oftype(A, choltype(A)) +cholcopy(A::RealHermSymComplexHerm{<:Any,<:Diagonal}) = Diagonal(copy_similar(diag(A), choltype(A))) function getproperty(C::Cholesky{<:Any,<:Diagonal}, d::Symbol) Cfactors = getfield(C, :factors) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index a95a73dfc8819..d0013aa553929 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -502,7 +502,7 @@ true ``` """ hessenberg(A::AbstractMatrix{T}) where T = - hessenberg!(copymutable_oftype(A, eigtype(T))) + hessenberg!(eigencopy_oftype(A, eigtype(T))) function show(io::IO, mime::MIME"text/plain", F::Hessenberg) summary(io, F) diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index 52d4f944f682f..81c34447402d7 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -120,7 +120,7 @@ julia> l == S.L && q == S.Q true ``` """ -lq(A::AbstractMatrix{T}) where {T} = lq!(copymutable_oftype(A, lq_eltype(T))) +lq(A::AbstractMatrix{T}) where {T} = lq!(copy_similar(A, lq_eltype(T))) lq(x::Number) = lq!(fill(convert(lq_eltype(typeof(x)), x), 1, 1)) lq_eltype(::Type{T}) where {T} = typeof(zero(T) / sqrt(abs2(one(T)))) @@ -195,9 +195,9 @@ function lmul!(A::LQ, B::StridedVecOrMat) lmul!(LowerTriangular(A.L), view(lmul!(A.Q, B), 1:size(A,1), axes(B,2))) return B end -function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB} +function *(A::LQ{TA}, B::AbstractVecOrMat{TB}) where {TA,TB} TAB = promote_type(TA, TB) - _cut_B(lmul!(convert(Factorization{TAB}, A), copymutable_oftype(B, TAB)), 1:size(A,1)) + _cut_B(lmul!(convert(Factorization{TAB}, A), copy_similar(B, TAB)), 1:size(A,1)) end ## Multiplication by Q diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 47e3fbfcb0232..3577f11d23399 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -461,18 +461,18 @@ end function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(A), eltype(F)) - return adjoint(ldiv!(F.parent, copymutable_oftype(adjoint(A), T))) + return adjoint(ldiv!(F.parent, copy_similar(adjoint(A), T))) end # To avoid ambiguities with definitions in adjtrans.jl and factorizations.jl (/)(adjA::Adjoint{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent) (/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) = adjoint(F.parent \ adjA.parent) function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(trA.parent, T)))) + return adjoint(ldiv!(F.parent, conj!(copy_similar(trA.parent, T)))) end function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:LU}) T = promote_type(eltype(trA), eltype(F)) - return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(trA.parent, T)))) + return adjoint(ldiv!(F.parent, conj!(copy_similar(trA.parent, T)))) end function det(F::LU{T}) where T diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 7fcace8e4ef71..8af8625a0e817 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -43,8 +43,8 @@ Bidiagonal(A::AbstractTriangular) = isbanded(A, -1, 0) ? Bidiagonal(diag(A, 0), diag(A, -1), :L) : # is lower bidiagonal throw(ArgumentError("matrix cannot be represented as Bidiagonal")) -_lucopy(A::Bidiagonal, T) = copymutable_oftype(Tridiagonal(A), T) -_lucopy(A::Diagonal, T) = copymutable_oftype(Tridiagonal(A), T) +_lucopy(A::Bidiagonal, T) = copymutable_oftype(Tridiagonal(A), T) +_lucopy(A::Diagonal, T) = copymutable_oftype(Tridiagonal(A), T) function _lucopy(A::SymTridiagonal, T) du = copy_similar(_evview(A), T) dl = copy.(transpose.(du)) diff --git a/stdlib/LinearAlgebra/src/svd.jl b/stdlib/LinearAlgebra/src/svd.jl index c58c83bcb5a98..a4d83edb50f13 100644 --- a/stdlib/LinearAlgebra/src/svd.jl +++ b/stdlib/LinearAlgebra/src/svd.jl @@ -175,11 +175,11 @@ julia> Uonly == U true ``` """ -function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T} - svd!(copymutable_oftype(A, eigtype(T)), full = full, alg = alg) +function svd(A::AbstractVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T} + svd!(eigencopy_oftype(A, eigtype(T)), full = full, alg = alg) end -function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T <: Union{Float16,Complex{Float16}}} - A = svd!(copymutable_oftype(A, eigtype(T)), full = full, alg = alg) +function svd(A::AbstractVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T <: Union{Float16,Complex{Float16}}} + A = svd!(eigencopy_oftype(A, eigtype(T)), full = full, alg = alg) return SVD{T}(A) end function svd(x::Number; full::Bool = false, alg::Algorithm = default_svd_alg(x)) @@ -240,10 +240,8 @@ julia> svdvals(A) 0.0 ``` """ -svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(copymutable_oftype(A, eigtype(T))) +svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(eigencopy_oftype(A, eigtype(T))) svdvals(A::AbstractVector{T}) where {T} = [convert(eigtype(T), norm(A))] -svdvals(A::AbstractMatrix{<:BlasFloat}) = svdvals!(copy(A)) -svdvals(A::AbstractVector{<:BlasFloat}) = [norm(A)] svdvals(x::Number) = abs(x) svdvals(S::SVD{<:Any,T}) where {T} = (S.S)::Vector{T} @@ -457,9 +455,9 @@ julia> U == Uonly true ``` """ -function svd(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} +function svd(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB} S = promote_type(eigtype(TA),TB) - return svd!(copymutable_oftype(A, S), copymutable_oftype(B, S)) + return svd!(copy_similar(A, S), copy_similar(B, S)) end # This method can be heavily optimized but it is probably not critical # and might introduce bugs or inconsistencies relative to the 1x1 matrix @@ -541,7 +539,6 @@ function svdvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasFloat end a[1:k + l] ./ b[1:k + l] end -svdvals(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:BlasFloat} = svdvals!(copy(A),copy(B)) """ svdvals(A, B) @@ -567,9 +564,9 @@ julia> svdvals(A, B) 1.0 ``` """ -function svdvals(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB} +function svdvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB} S = promote_type(eigtype(TA), TB) - return svdvals!(copymutable_oftype(A, S), copymutable_oftype(B, S)) + return svdvals!(copy_similar(A, S), copy_similar(B, S)) end svdvals(x::Number, y::Number) = abs(x/y) diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index c7ca6339aac6a..f15f54669d124 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -201,3 +201,12 @@ function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_de end return B end + +function copy_similar(A::AdjointAbsMat, ::Type{T}) where {T} + C = similar(A, T, size(A)) + adjoint!(C, parent(A)) +end +function copy_similar(A::TransposeAbsMat, ::Type{T}) where {T} + C = similar(A, T, size(A)) + transpose!(C, parent(A)) +end diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl index 448ab19b2fb58..a3008a236df7b 100644 --- a/stdlib/LinearAlgebra/test/cholesky.jl +++ b/stdlib/LinearAlgebra/test/cholesky.jl @@ -390,9 +390,9 @@ end # complex D = complex(D) - CD = cholesky(D) - CM = cholesky(Matrix(D)) - @test CD isa Cholesky{ComplexF64} + CD = cholesky(Hermitian(D)) + CM = cholesky(Matrix(Hermitian(D))) + @test CD isa Cholesky{ComplexF64,<:Diagonal} @test CD.U ≈ Diagonal(.√d) ≈ CM.U @test D ≈ CD.L * CD.U @test CD.info == 0 diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl index b2b23caac6865..4b14179e644e5 100644 --- a/stdlib/LinearAlgebra/test/hessenberg.jl +++ b/stdlib/LinearAlgebra/test/hessenberg.jl @@ -191,6 +191,13 @@ let n = 10 end end +@testset "hessenberg(::AbstractMatrix)" begin + n = 10 + A = Tridiagonal(rand(n-1), rand(n), rand(n-1)) + H = hessenberg(A) + @test convert(Array, H) ≈ A +end + # check logdet on a matrix that has a positive determinant let A = [0.5 0.1 0.9 0.4; 0.9 0.7 0.5 0.4; 0.3 0.4 0.9 0.0; 0.4 0.0 0.0 0.5] @test logdet(hessenberg(A)) ≈ logdet(A) ≈ -3.5065578973199822 diff --git a/stdlib/LinearAlgebra/test/lq.jl b/stdlib/LinearAlgebra/test/lq.jl index 96f31ded78d6d..c340317a7cc23 100644 --- a/stdlib/LinearAlgebra/test/lq.jl +++ b/stdlib/LinearAlgebra/test/lq.jl @@ -37,10 +37,10 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) @testset for isview in (false,true) let a = isview ? view(a, 1:m - 1, 1:n - 1) : a, b = isview ? view(b, 1:m - 1) : b, m = m - isview, n = n - isview - lqa = lq(a) + lqa = lq(a) x = lqa\b - l,q = lqa.L, lqa.Q - qra = qr(a, ColumnNorm()) + l, q = lqa.L, lqa.Q + qra = qr(a, ColumnNorm()) @testset "Basic ops" begin @test size(lqa,1) == size(a,1) @test size(lqa,3) == 1 @@ -62,6 +62,10 @@ rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q) @test Array{eltya}(q) ≈ Matrix(q) end @testset "Binary ops" begin + k = size(a, 2) + T = Tridiagonal(rand(eltya, k-1), rand(eltya, k), rand(eltya, k-1)) + @test lq(T) * T ≈ T * T rtol=3000ε + @test lqa * T ≈ a * T rtol=3000ε @test a*x ≈ b rtol=3000ε @test x ≈ qra \ b rtol=3000ε @test lqa*x ≈ a*x rtol=3000ε diff --git a/stdlib/LinearAlgebra/test/svd.jl b/stdlib/LinearAlgebra/test/svd.jl index 8bd3edadc911d..7f2aad904a88f 100644 --- a/stdlib/LinearAlgebra/test/svd.jl +++ b/stdlib/LinearAlgebra/test/svd.jl @@ -127,8 +127,20 @@ aimg = randn(n,n)/2 gsvd = svd(b,c) @test gsvd.U*gsvd.D1*gsvd.R*gsvd.Q' ≈ b @test gsvd.V*gsvd.D2*gsvd.R*gsvd.Q' ≈ c + # AbstractMatrix svd + T = Tridiagonal(a) + asvd = svd(T, a) + @test asvd.U*asvd.D1*asvd.R*asvd.Q' ≈ T + @test asvd.V*asvd.D2*asvd.R*asvd.Q' ≈ a + @test all(≈(1), svdvals(T, T)) end end + @testset "singular value decomposition of AbstractMatrix" begin + A = Tridiagonal(aa) + F = svd(A) + @test Matrix(F) ≈ A + @test svdvals(A) ≈ F.S + end @testset "singular value decomposition of Hermitian/real-Symmetric" begin for T in (eltya <: Real ? (Symmetric, Hermitian) : (Hermitian,)) usv = svd(T(asym))