audit copymutable_oftype usage (#47063)

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

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)

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)

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

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

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

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)
 

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

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

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

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ε

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