RFC: Rework copy_oftype a bit (#44756)

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -356,44 +356,47 @@ control over the factorization of `B`.
 """
 rdiv!(A, B)
 
+"""
+    copy_oftype(A, T)
 
+Creates a copy of `A` with eltype `T`. No assertions about mutability of the result are
+made. When `eltype(A) == T`, then this calls `copy(A)` which may be overloaded for custom
+array types. Otherwise, this calls `convert(AbstractArray{T}, A)`.
+"""
+copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
+copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)
 
 """
-    copy_oftype(A, T)
+    copymutable_oftype(A, T)
 
 Copy `A` to a mutable array with eltype `T` based on `similar(A, T)`.
 
 The resulting matrix typically has similar algebraic structure as `A`. For
 example, supplying a tridiagonal matrix results in another tridiagonal matrix.
 In general, the type of the output corresponds to that of `similar(A, T)`.
 
-There are three often used methods in LinearAlgebra to create a mutable copy
-of an array with a given eltype. These copies can be passed to in-place
-algorithms (such as `ldiv!`, `rdiv!`, `lu!` and so on). Which one to use in practice
-depends on what is known (or assumed) about the structure of the array in that
-algorithm.
+In LinearAlgebra, mutable copies (of some desired eltype) are created to be passed
+to in-place algorithms (such as `ldiv!`, `rdiv!`, `lu!` and so on). If the specific
+algorithm is known to preserve the algebraic structure, use `copymutable_oftype`.
+If the algorithm is known to return a dense matrix (or some wrapper backed by a dense
+matrix), then use `copy_similar`.
 
-See also: `copy_similar`.
+See also: `Base.copymutable`, `copy_similar`.
 """
-copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T), A)
+copymutable_oftype(A::AbstractArray, ::Type{S}) where {S} = copyto!(similar(A, S), A)
 
 """
     copy_similar(A, T)
 
 Copy `A` to a mutable array with eltype `T` based on `similar(A, T, size(A))`.
 
-Compared to `copy_oftype`, the result can be more flexible. In general, the type
+Compared to `copymutable_oftype`, the result can be more flexible. In general, the type
 of the output corresponds to that of the three-argument method `similar(A, T, size(A))`.
 
-See also: `copy_oftype`.
+See also: `copymutable_oftype`.
 """
 copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A)
 
-# The three copy functions above return mutable arrays with eltype T.
-# To only ensure a certain eltype, and if a mutable copy is not needed, it is
-# more efficient to use:
-# convert(AbstractArray{T}, A)
-
 
 include("adjtrans.jl")
 include("transpose.jl")

stdlib/LinearAlgebra/src/bunchkaufman.jl

@@ -197,7 +197,7 @@ julia> S.L*S.D*S.L' - A[S.p, S.p]
 ```
 """
 bunchkaufman(A::AbstractMatrix{T}, rook::Bool=false; check::Bool = true) where {T} =
-    bunchkaufman!(copy_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
+    bunchkaufman!(copymutable_oftype(A, typeof(sqrt(oneunit(T)))), rook; check = check)
 
 BunchKaufman{T}(B::BunchKaufman) where {T} =
     BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)

stdlib/LinearAlgebra/src/cholesky.jl

@@ -179,8 +179,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) = copy_oftype(A, choltype(A))
-@inline cholcopy(A::RealHermSymComplexHerm) = copy_oftype(A, choltype(A))
+@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))
 
 # _chol!. Internal methods for calling unpivoted Cholesky

stdlib/LinearAlgebra/src/dense.jl

@@ -500,7 +500,7 @@ function (^)(A::AbstractMatrix{T}, p::Real) where T
     # Quicker return if A is diagonal
     if isdiag(A)
         TT = promote_op(^, T, typeof(p))
-        retmat = copy_oftype(A, TT)
+        retmat = copymutable_oftype(A, TT)
         for i in 1:n
             retmat[i, i] = retmat[i, i] ^ p
         end

stdlib/LinearAlgebra/src/diagonal.jl

@@ -792,8 +792,8 @@ end
 @deprecate cholesky!(A::Diagonal, ::Val{false}; check::Bool = true) cholesky!(A::Diagonal, NoPivot(); check) false
 @deprecate cholesky(A::Diagonal, ::Val{false}; check::Bool = true) cholesky(A::Diagonal, NoPivot(); check) false
 
-@inline cholcopy(A::Diagonal) = copy_oftype(A, choltype(A))
-@inline cholcopy(A::RealHermSymComplexHerm{<:Real,<:Diagonal}) = copy_oftype(A, choltype(A))
+@inline cholcopy(A::Diagonal) = copymutable_oftype(A, choltype(A))
+@inline cholcopy(A::RealHermSymComplexHerm{<:Real,<:Diagonal}) = copymutable_oftype(A, choltype(A))
 
 function getproperty(C::Cholesky{<:Any,<:Diagonal}, d::Symbol)
     Cfactors = getfield(C, :factors)

stdlib/LinearAlgebra/src/eigen.jl

@@ -233,12 +233,12 @@ true
 ```
 """
 function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
-    AA = copy_oftype(A, eigtype(T))
+    AA = copymutable_oftype(A, eigtype(T))
     isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
     return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
 end
 function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where {T <: Union{Float16,Complex{Float16}}}
-    AA = copy_oftype(A, eigtype(T))
+    AA = copymutable_oftype(A, eigtype(T))
     isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
     A = eigen!(AA; permute, scale, sortby)
     values = convert(AbstractVector{isreal(A.values) ? Float16 : Complex{Float16}}, A.values)
@@ -333,7 +333,7 @@ julia> eigvals(diag_matrix)
 ```
 """
 eigvals(A::AbstractMatrix{T}; kws...) where T =
-    eigvals!(copy_oftype(A, eigtype(T)); kws...)
+    eigvals!(copymutable_oftype(A, eigtype(T)); kws...)
 
 """
 For a scalar input, `eigvals` will return a scalar.
@@ -508,7 +508,7 @@ true
 """
 function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
-    eigen!(copy_oftype(A, S), copy_oftype(B, S); kws...)
+    eigen!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
 end
 
 eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
@@ -587,7 +587,7 @@ julia> eigvals(A,B)
 """
 function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
-    return eigvals!(copy_oftype(A, S), copy_oftype(B, S); kws...)
+    return eigvals!(copymutable_oftype(A, S), copymutable_oftype(B, S); kws...)
 end
 
 """

stdlib/LinearAlgebra/src/generic.jl

@@ -1779,7 +1779,7 @@ julia> normalize(a)
 function normalize(a::AbstractArray, p::Real = 2)
     nrm = norm(a, p)
     if !isempty(a)
-        aa = copy_oftype(a, typeof(first(a)/nrm))
+        aa = copymutable_oftype(a, typeof(first(a)/nrm))
         return __normalize!(aa, nrm)
     else
         T = typeof(zero(eltype(a))/nrm)

stdlib/LinearAlgebra/src/givens.jl

@@ -8,7 +8,7 @@ transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R
 
 function (*)(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S}
     TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
-    lmul!(convert(AbstractRotation{TS}, R), copy_oftype(A, TS))
+    lmul!(convert(AbstractRotation{TS}, R), copy_similar(A, TS))
 end
 (*)(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
 (*)(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)

stdlib/LinearAlgebra/src/hessenberg.jl

@@ -502,7 +502,7 @@ true
 ```
 """
 hessenberg(A::AbstractMatrix{T}) where T =
-    hessenberg!(copy_oftype(A, eigtype(T)))
+    hessenberg!(copymutable_oftype(A, eigtype(T)))
 
 function show(io::IO, mime::MIME"text/plain", F::Hessenberg)
     summary(io, F)

stdlib/LinearAlgebra/src/ldlt.jl

@@ -162,7 +162,7 @@ julia> S \\ b
 """
 function ldlt(M::SymTridiagonal{T}; shift::Number=false) where T
     S = typeof((zero(T)+shift)/one(T))
-    Mₛ = SymTridiagonal{S}(copy_oftype(M.dv, S), copy_oftype(M.ev, S))
+    Mₛ = SymTridiagonal{S}(copymutable_oftype(M.dv, S), copymutable_oftype(M.ev, S))
     if !iszero(shift)
         Mₛ.dv .+= shift
     end

stdlib/LinearAlgebra/src/lq.jl

@@ -120,7 +120,7 @@ julia> l == S.L &&  q == S.Q
 true
 ```
 """
-lq(A::AbstractMatrix{T}) where {T}  = lq!(copy_oftype(A, lq_eltype(T)))
+lq(A::AbstractMatrix{T}) where {T}  = lq!(copymutable_oftype(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))))
@@ -197,15 +197,15 @@ function lmul!(A::LQ, B::StridedVecOrMat)
 end
 function *(A::LQ{TA}, B::StridedVecOrMat{TB}) where {TA,TB}
     TAB = promote_type(TA, TB)
-    _cut_B(lmul!(convert(Factorization{TAB}, A), copy_oftype(B, TAB)), 1:size(A,1))
+    _cut_B(lmul!(convert(Factorization{TAB}, A), copymutable_oftype(B, TAB)), 1:size(A,1))
 end
 
 ## Multiplication by Q
 ### QB
 lmul!(A::LQPackedQ{T}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = LAPACK.ormlq!('L','N',A.factors,A.τ,B)
 function (*)(A::LQPackedQ, B::StridedVecOrMat)
     TAB = promote_type(eltype(A), eltype(B))
-    lmul!(AbstractMatrix{TAB}(A), copy_oftype(B, TAB))
+    lmul!(AbstractMatrix{TAB}(A), copymutable_oftype(B, TAB))
 end
 
 ### QcB
@@ -218,7 +218,7 @@ function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
     A = adjA.parent
     TAB = promote_type(eltype(A), eltype(B))
     if size(B,1) == size(A.factors,2)
-        lmul!(adjoint(AbstractMatrix{TAB}(A)), copy_oftype(B, TAB))
+        lmul!(adjoint(AbstractMatrix{TAB}(A)), copymutable_oftype(B, TAB))
     elseif size(B,1) == size(A.factors,1)
         lmul!(adjoint(AbstractMatrix{TAB}(A)), [B; zeros(TAB, size(A.factors, 2) - size(A.factors, 1), size(B, 2))])
     else
@@ -269,7 +269,7 @@ rmul!(A::StridedMatrix{T}, adjB::Adjoint{<:Any,<:LQPackedQ{T}}) where {T<:BlasCo
 function *(A::StridedVecOrMat, adjQ::Adjoint{<:Any,<:LQPackedQ})
     Q = adjQ.parent
     TR = promote_type(eltype(A), eltype(Q))
-    return rmul!(copy_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
+    return rmul!(copymutable_oftype(A, TR), adjoint(AbstractMatrix{TR}(Q)))
 end
 function *(adjA::Adjoint{<:Any,<:StridedMatrix}, adjQ::Adjoint{<:Any,<:LQPackedQ})
     A, Q = adjA.parent, adjQ.parent
@@ -293,7 +293,7 @@ end
 function *(A::StridedVecOrMat, Q::LQPackedQ)
     TR = promote_type(eltype(A), eltype(Q))
     if size(A, 2) == size(Q.factors, 2)
-        C = copy_oftype(A, TR)
+        C = copymutable_oftype(A, TR)
     elseif size(A, 2) == size(Q.factors, 1)
         C = zeros(TR, size(A, 1), size(Q.factors, 2))
         copyto!(C, 1, A, 1, length(A))

stdlib/LinearAlgebra/src/lu.jl

@@ -283,8 +283,8 @@ end
 @deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, NoPivot(); check=check)
 
 _lucopy(A::AbstractMatrix, T) = copy_similar(A, T)
-_lucopy(A::HermOrSym, T)      = copy_oftype(A, T)
-_lucopy(A::Tridiagonal, T)    = copy_oftype(A, T)
+_lucopy(A::HermOrSym, T)      = copymutable_oftype(A, T)
+_lucopy(A::Tridiagonal, T)    = copymutable_oftype(A, T)
 
 lu(S::LU) = S
 function lu(x::Number; check::Bool=true)
@@ -438,18 +438,18 @@ end
 
 function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:LU})
     T = promote_type(eltype(A), eltype(F))
-    return adjoint(ldiv!(F.parent, copy_oftype(adjoint(A), T)))
+    return adjoint(ldiv!(F.parent, copymutable_oftype(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!(copy_oftype(trA.parent, T))))
+    return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(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!(copy_oftype(trA.parent, T))))
+    return adjoint(ldiv!(F.parent, conj!(copymutable_oftype(trA.parent, T))))
 end
 
 function det(F::LU{T}) where T

stdlib/LinearAlgebra/src/matmul.jl

@@ -171,7 +171,7 @@ end
 function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMat{<:BlasReal})
     TS = promote_type(eltype(A), eltype(B))
     mul!(similar(B, TS, (size(A, 1), size(B, 2))),
-         copy_oftype(A, TS), # remove AdjOrTrans to use reinterpret trick below
+         copymutable_oftype(A, TS), # remove AdjOrTrans to use reinterpret trick below
          wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B))))
 end
 # the following case doesn't seem to benefit from the translation A*B = (B' * A')'

stdlib/LinearAlgebra/src/qr.jl

@@ -657,7 +657,7 @@ function (*)(A::AbstractQ, b::StridedVector)
     TAb = promote_type(eltype(A), eltype(b))
     Anew = convert(AbstractMatrix{TAb}, A)
     if size(A.factors, 1) == length(b)
-        bnew = copy_oftype(b, TAb)
+        bnew = copymutable_oftype(b, TAb)
     elseif size(A.factors, 2) == length(b)
         bnew = [b; zeros(TAb, size(A.factors, 1) - length(b))]
     else
@@ -669,7 +669,7 @@ function (*)(A::AbstractQ, B::StridedMatrix)
     TAB = promote_type(eltype(A), eltype(B))
     Anew = convert(AbstractMatrix{TAB}, A)
     if size(A.factors, 1) == size(B, 1)
-        Bnew = copy_oftype(B, TAB)
+        Bnew = copymutable_oftype(B, TAB)
     elseif size(A.factors, 2) == size(B, 1)
         Bnew = [B; zeros(TAB, size(A.factors, 1) - size(B,1), size(B, 2))]
     else
@@ -723,7 +723,7 @@ end
 function *(adjQ::Adjoint{<:Any,<:AbstractQ}, B::StridedVecOrMat)
     Q = adjQ.parent
     TQB = promote_type(eltype(Q), eltype(B))
-    return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copy_oftype(B, TQB))
+    return lmul!(adjoint(convert(AbstractMatrix{TQB}, Q)), copymutable_oftype(B, TQB))
 end
 
 ### QBc/QcBc
@@ -775,7 +775,7 @@ end
 function (*)(A::StridedMatrix, Q::AbstractQ)
     TAQ = promote_type(eltype(A), eltype(Q))
 
-    return rmul!(copy_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
+    return rmul!(copymutable_oftype(A, TAQ), convert(AbstractMatrix{TAQ}, Q))
 end
 
 function (*)(a::Number, B::AbstractQ)

stdlib/LinearAlgebra/src/special.jl

@@ -50,8 +50,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)     = copy_oftype(Tridiagonal(A), T)
-_lucopy(A::Diagonal, T)       = copy_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

@@ -176,10 +176,10 @@ true
 ```
 """
 function svd(A::StridedVecOrMat{T}; full::Bool = false, alg::Algorithm = default_svd_alg(A)) where {T}
-    svd!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
+    svd!(copymutable_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!(copy_oftype(A, eigtype(T)), full = full, alg = alg)
+    A = svd!(copymutable_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,7 +240,7 @@ julia> svdvals(A)
  0.0
 ```
 """
-svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(copy_oftype(A, eigtype(T)))
+svdvals(A::AbstractMatrix{T}) where {T} = svdvals!(copymutable_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)]
@@ -459,7 +459,7 @@ true
 """
 function svd(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
-    return svd!(copy_oftype(A, S), copy_oftype(B, S))
+    return svd!(copymutable_oftype(A, S), copymutable_oftype(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
@@ -569,7 +569,7 @@ julia> svdvals(A, B)
 """
 function svdvals(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
     S = promote_type(eigtype(TA), TB)
-    return svdvals!(copy_oftype(A, S), copy_oftype(B, S))
+    return svdvals!(copymutable_oftype(A, S), copymutable_oftype(B, S))
 end
 svdvals(x::Number, y::Number) = abs(x/y)