Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Differentiate between copy_oftype(A, T) and convert(AbstractArray{T}, A) for structured arrays #35165

Closed
wants to merge 4 commits into from
Closed

Differentiate between copy_oftype(A, T) and convert(AbstractArray{T}, A) for structured arrays #35165

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eb30779
Redefine copy_oftype, add AbstractArray constructores
daanhb Mar 18, 2020
e383d7e
implement and use copy_similar
daanhb Mar 18, 2020
2a61d67
add tests for abstractarray conversion
daanhb Mar 25, 2020
8e81917
move ImmutableArray to testhelpers
daanhb Apr 25, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 16 additions & 1 deletion stdlib/LinearAlgebra/src/LinearAlgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,8 +346,22 @@ control over the factorization of `B`.
"""
rdiv!(A, B)

"""
Copy array `A` and make sure the result has element type `T`.

This method often preserves the structure of `A`. It does so if both `copy(A)`
and the two-argument function `similar(A, T)` do.
"""
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::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = copyto!(similar(A, S), A)

"""
Copy array `A` and make sure the result has element type `T`.

This method typically discards the structure of `A`. It does so if the
three-argument function `similar(A, T, size(A))` does.
"""
copy_similar(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = copyto!(similar(A, S, size(A)), A)

include("adjtrans.jl")
include("transpose.jl")
Expand Down Expand Up @@ -450,4 +464,5 @@ function __init__()
Base.at_disable_library_threading(() -> BLAS.set_num_threads(1))
end


end # module LinearAlgebra
4 changes: 4 additions & 0 deletions stdlib/LinearAlgebra/src/adjtrans.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,10 @@ similar(A::AdjOrTrans) = similar(A.parent, eltype(A), axes(A))
similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, axes(A))
similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims)

# mimick the style of similar for the AbstractMatrix{T} constructor:
# for vectors we maintain the parent and wrapper type here, matrices use the default
AbstractMatrix{T}(A::AdjOrTransAbsVec) where {T} = wrapperop(A)(AbstractVector{T}(A.parent))

# sundry basic definitions
parent(A::AdjOrTrans) = A.parent
vec(v::TransposeAbsVec) = parent(v)
Expand Down
16 changes: 8 additions & 8 deletions stdlib/LinearAlgebra/src/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,14 +178,14 @@ end
(*)(D::Diagonal, V::AbstractVector) = D.diag .* V

(*)(A::AbstractTriangular, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A), D)
rmul!(copy_oftype(A, promote_op(*, eltype(A), eltype(D.diag))), D)
(*)(D::Diagonal, B::AbstractTriangular) =
lmul!(D, copyto!(similar(B, promote_op(*, eltype(B), eltype(D.diag))), B))
lmul!(D, copy_oftype(B, promote_op(*, eltype(B), eltype(D.diag))))

(*)(A::AbstractMatrix, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A), D)
rmul!(copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))), D)
(*)(D::Diagonal, A::AbstractMatrix) =
lmul!(D, copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A))
lmul!(D, copy_similar(A, promote_op(*, eltype(A), eltype(D.diag))))

function rmul!(A::AbstractMatrix, D::Diagonal)
require_one_based_indexing(A)
Expand Down Expand Up @@ -232,7 +232,7 @@ end

*(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag)
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A), D)
rmul!(copy_oftype(A, promote_op(*, eltype(A), eltype(D.diag))), D)
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)))
Expand All @@ -242,7 +242,7 @@ end

*(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag)
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) =
rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A), D)
rmul!(copy_oftype(A, promote_op(*, eltype(A), eltype(D.diag))), D)
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)))
Expand All @@ -252,7 +252,7 @@ end

*(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag))
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) =
lmul!(D, copyto!(similar(B, promote_op(*, eltype(B), eltype(D.diag))), B))
lmul!(D, copy_oftype(B, promote_op(*, eltype(B), eltype(D.diag))))
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; rmul!(Array(D), adjoint(Q)))
function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
A = adjA.parent
Expand All @@ -263,7 +263,7 @@ end

*(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag))
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) =
lmul!(D, copyto!(similar(B, promote_op(*, eltype(B), eltype(D.diag))), B))
lmul!(D, copy_oftype(B, promote_op(*, eltype(B), eltype(D.diag))))
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)))
Expand Down
12 changes: 4 additions & 8 deletions stdlib/LinearAlgebra/src/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,32 +96,28 @@ end
function \(F::Factorization, B::AbstractVecOrMat)
require_one_based_indexing(B)
TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
BB = similar(B, TFB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TFB)
ldiv!(F, BB)
end
function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat)
require_one_based_indexing(B)
F = adjF.parent
TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
BB = similar(B, TFB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TFB)
ldiv!(adjoint(F), BB)
end

function /(B::AbstractMatrix, F::Factorization)
require_one_based_indexing(B)
TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
BB = similar(B, TFB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TFB)
rdiv!(BB, F)
end
function /(B::AbstractMatrix, adjF::Adjoint{<:Any,<:Factorization})
require_one_based_indexing(B)
F = adjF.parent
TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
BB = similar(B, TFB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TFB)
rdiv!(BB, adjoint(F))
end
/(adjB::AdjointAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjB.parent)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/givens.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
lmul!(convert(AbstractRotation{TS}, R), copy_oftype(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)
Expand Down
2 changes: 2 additions & 0 deletions stdlib/LinearAlgebra/src/hessenberg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,8 @@ parent(H::UpperHessenberg) = H.data
similar(H::UpperHessenberg, ::Type{T}) where {T} = UpperHessenberg(similar(H.data, T))
similar(H::UpperHessenberg, ::Type{T}, dims::Dims{N}) where {T,N} = similar(H.data, T, dims)

AbstractMatrix{T}(H::UpperHessenberg) where {T} = UpperHessenberg(AbstractMatrix{T}(H.data))

copy(H::UpperHessenberg) = UpperHessenberg(copy(H.data))
real(H::UpperHessenberg{<:Real}) = H
real(H::UpperHessenberg{<:Complex}) = UpperHessenberg(triu!(real(H.data),-1))
Expand Down
5 changes: 3 additions & 2 deletions stdlib/LinearAlgebra/src/lu.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,11 +432,12 @@ end
(/)(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, convert(AbstractVector{T}, conj(trA.parent))))
# Make a copy here to avoid overwriting trA.parent if it is real, #35071
return adjoint(ldiv!(F.parent, copy_oftype(conj(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, convert(AbstractMatrix{T}, conj(trA.parent))))
return adjoint(ldiv!(F.parent, copy_oftype(conj(trA.parent), T)))
end

function det(F::LU{T}) where T
Expand Down
6 changes: 2 additions & 4 deletions stdlib/LinearAlgebra/src/qr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -378,8 +378,7 @@ true
"""
function qr(A::AbstractMatrix{T}, arg...; kwargs...) where T
require_one_based_indexing(A)
AA = similar(A, _qreltype(T), size(A))
copyto!(AA, A)
AA = copy_similar(A, _qreltype(T))
return qr!(AA, arg...; kwargs...)
end
qr(x::Number) = qr(fill(x,1,1))
Expand Down Expand Up @@ -720,8 +719,7 @@ function *(A::StridedMatrix, adjB::Adjoint{<:Any,<:AbstractQ})
TAB = promote_type(eltype(A),eltype(B))
BB = convert(AbstractMatrix{TAB}, B)
if size(A,2) == size(B.factors, 1)
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return rmul!(AA, adjoint(BB))
elseif size(A,2) == size(B.factors,2)
return rmul!([A zeros(TAB, size(A, 1), size(B.factors, 1) - size(B.factors, 2))], adjoint(BB))
Expand Down
12 changes: 6 additions & 6 deletions stdlib/LinearAlgebra/src/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -674,7 +674,7 @@ eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; sortby::Union{Func
function eigen(A::RealHermSymComplexHerm; sortby::Union{Function,Nothing}=nothing)
T = eltype(A)
S = eigtype(T)
eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby=sortby)
eigen!(copy_oftype(A, S), sortby=sortby)
end

eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}, irange::UnitRange) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.data, 0.0, 0.0, irange.start, irange.stop, -1.0)...)
Expand All @@ -699,7 +699,7 @@ The [`UnitRange`](@ref) `irange` specifies indices of the sorted eigenvalues to
function eigen(A::RealHermSymComplexHerm, irange::UnitRange)
T = eltype(A)
S = eigtype(T)
eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
eigen!(copy_oftype(A, S), irange)
end

eigen!(A::RealHermSymComplexHerm{T,<:StridedMatrix}, vl::Real, vh::Real) where {T<:BlasReal} =
Expand All @@ -725,7 +725,7 @@ The following functions are available for `Eigen` objects: [`inv`](@ref), [`det`
function eigen(A::RealHermSymComplexHerm, vl::Real, vh::Real)
T = eltype(A)
S = eigtype(T)
eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
eigen!(copy_oftype(A, S), vl, vh)
end

eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) =
Expand All @@ -734,7 +734,7 @@ eigvals!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) =
function eigvals(A::RealHermSymComplexHerm)
T = eltype(A)
S = eigtype(T)
eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A))
eigvals!(copy_oftype(A, S))
end

"""
Expand Down Expand Up @@ -775,7 +775,7 @@ julia> eigvals(A)
function eigvals(A::RealHermSymComplexHerm, irange::UnitRange)
T = eltype(A)
S = eigtype(T)
eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
eigvals!(copy_oftype(A, S), irange)
end

"""
Expand Down Expand Up @@ -815,7 +815,7 @@ julia> eigvals(A)
function eigvals(A::RealHermSymComplexHerm, vl::Real, vh::Real)
T = eltype(A)
S = eigtype(T)
eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
eigvals!(copy_oftype(A, S), vl, vh)
end

eigmax(A::RealHermSymComplexHerm{<:Real,<:StridedMatrix}) = eigvals(A, size(A, 1):size(A, 1))[1]
Expand Down
Loading