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

make AbstractArray{T}(...) constructor in LinearAlgebra more consistent #40831

Merged
merged 4 commits into from
Jun 19, 2021
Merged

make AbstractArray{T}(...) constructor in LinearAlgebra more consistent #40831

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
87891d7
improve AbstractArray{T} constructor in LinearAlgebra
daanhb Jun 17, 2021
0ee05bc
document the copy functions of LinearAlgebra (#40831)
daanhb Jun 17, 2021
72ac4bf
test lu factorization fallback (#40831)
daanhb Jun 17, 2021
ada4020
tests for lu, docstrings (#40831)
daanhb Jun 18, 2021
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
26 changes: 24 additions & 2 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -354,8 +354,30 @@ control over the factorization of `B`.
"""
rdiv!(A, B)

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)

# 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.

# copy_oftype: make a mutable copy 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.
copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A,T), A)
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be renamed copymutable_oftype as it behaves very different than copy (which is not in general mutable)


# copy_similar: make a mutable copy based on `similar(A, T, size(A))`.
# Compared to copy_oftype, the result can be more flexible. For example,
# supplying a tridiagonal matrix results in a sparse array.
copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A)

# copy_to_array: copy the contents of A to a standard dense array with element type T
copy_to_array(A::AbstractArray, ::Type{T}) where {T} = copyto!(Array{T}(undef, size(A)...), A)

# 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")
Expand Down
3 changes: 3 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 @@ -210,6 +210,9 @@ 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)

# AbstractMatrix{T} constructor for adjtrans vector: preserve wrapped type
AbstractMatrix{T}(A::AdjOrTransAbsVec) where {T} = wrapperop(A)(AbstractVector{T}(A.parent))

# sundry basic definitions
parent(A::AdjOrTrans) = A.parent
vec(v::TransposeAbsVec{<:Number}) = parent(v)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/dense.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -451,7 +451,7 @@ function (^)(A::AbstractMatrix{T}, p::Integer) where T<:Integer
end
function integerpow(A::AbstractMatrix{T}, p) where T
TT = promote_op(^, T, typeof(p))
return (TT == T ? A : copyto!(similar(A, TT), A))^Integer(p)
return (TT == T ? A : convert(AbstractMatrix{TT}, A))^Integer(p)
end
function schurpow(A::AbstractMatrix, p)
if istriu(A)
Expand Down
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,14 +215,14 @@ function (*)(D::Diagonal, V::AbstractVector)
end

(*)(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
8 changes: 2 additions & 6 deletions stdlib/LinearAlgebra/src/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,17 +102,13 @@ end
function \(F::Union{Factorization, Adjoint{<:Any,<: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)
ldiv!(F, BB)
ldiv!(F, copy_similar(B, TFB))
end

function /(B::AbstractMatrix, F::Union{Factorization, Adjoint{<:Any,<:Factorization}})
require_one_based_indexing(B)
TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F)))
BB = similar(B, TFB, size(B))
copyto!(BB, B)
rdiv!(BB, F)
rdiv!(copy_similar(B, TFB), F)
end
/(adjB::AdjointAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjB.parent)
/(B::TransposeAbsVec, adjF::Adjoint{<:Any,<:Factorization}) = adjoint(adjF.parent \ adjoint(B))
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
12 changes: 11 additions & 1 deletion stdlib/LinearAlgebra/src/lu.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,9 @@ adjoint(F::LU) = Adjoint(F)
transpose(F::LU) = Transpose(F)

# StridedMatrix
lu(A::StridedMatrix, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) =
lu!(copy_oftype(A, lutype(eltype(A))), pivot; check=check)

lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMaximum(); check=check)
function lu!(A::StridedMatrix{T}, ::RowMaximum; check::Bool = true) where {T<:BlasFloat}
lpt = LAPACK.getrf!(A)
Expand All @@ -85,6 +88,10 @@ end
function lu!(A::StridedMatrix{<:BlasFloat}, pivot::NoPivot; check::Bool = true)
return generic_lufact!(A, pivot; check = check)
end

lu(A::HermOrSym, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) =
lu!(copy_oftype(A, lutype(eltype(A))), pivot; check=check)

function lu!(A::HermOrSym, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true)
copytri!(A.data, A.uplo, isa(A, Hermitian))
lu!(A.data, pivot; check = check)
Expand Down Expand Up @@ -276,7 +283,7 @@ true
"""
function lu(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T}
S = lutype(T)
lu!(copy_oftype(A, S), pivot; check = check)
lu!(copy_to_array(A, S), pivot; check = check)
end
# TODO: remove for Julia v2.0
@deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, RowMaximum(); check=check)
Expand Down Expand Up @@ -490,6 +497,9 @@ inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A))

# Tridiagonal

lu(A::Tridiagonal{T}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where T =
lu!(copy_oftype(A, lutype(T)), pivot; check = check)

# See dgttrf.f
function lu!(A::Tridiagonal{T,V}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T,V}
# Extract values
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 @@ -381,8 +381,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
# TODO: remove in Julia v2.0
Expand Down Expand Up @@ -770,8 +769,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
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/schur.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ true
schur(A::StridedMatrix{<:BlasFloat}) = schur!(copy(A))
schur(A::StridedMatrix{T}) where T = schur!(copy_oftype(A, eigtype(T)))

schur(A::AbstractMatrix{T}) where {T} = schur!(copyto!(Matrix{eigtype(T)}(undef, size(A)...), A))
schur(A::AbstractMatrix{T}) where {T} = schur!(copy_to_array(A, eigtype(T)))
function schur(A::RealHermSymComplexHerm)
F = eigen(A; sortby=nothing)
return Schur(typeof(F.vectors)(Diagonal(F.values)), F.vectors, F.values)
Expand Down
71 changes: 25 additions & 46 deletions stdlib/LinearAlgebra/src/triangular.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@ for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular,
end
Matrix(A::$t{T}) where {T} = Matrix{T}(A)

AbstractMatrix{T}(A::$t) where {T} = $t{T}(A)

size(A::$t, d) = size(A.data, d)
size(A::$t) = size(A.data)

Expand Down Expand Up @@ -1506,64 +1508,56 @@ for (f, f2!) in ((:*, :lmul!), (:\, :ldiv!))
function ($f)(A::LowerTriangular, B::LowerTriangular)
TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
($f)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function $(f)(A::UnitLowerTriangular, B::LowerTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function $(f)(A::LowerTriangular, B::UnitLowerTriangular)
TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
($f)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return LowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function $(f)(A::UnitLowerTriangular, B::UnitLowerTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return UnitLowerTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function ($f)(A::UpperTriangular, B::UpperTriangular)
TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
($f)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function ($f)(A::UnitUpperTriangular, B::UpperTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function ($f)(A::UpperTriangular, B::UnitUpperTriangular)
TAB = typeof(($f)(zero(eltype(A)), zero(eltype(B))) +
($f)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return UpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end

function ($f)(A::UnitUpperTriangular, B::UnitUpperTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
return UnitUpperTriangular($f2!(convert(AbstractMatrix{TAB}, A), BB))
end
end
Expand All @@ -1572,66 +1566,57 @@ end
function (/)(A::LowerTriangular, B::LowerTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return LowerTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UnitLowerTriangular, B::LowerTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return LowerTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::LowerTriangular, B::UnitLowerTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return LowerTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UnitLowerTriangular, B::UnitLowerTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return UnitLowerTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UpperTriangular, B::UpperTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return UpperTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UnitUpperTriangular, B::UpperTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return UpperTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UpperTriangular, B::UnitUpperTriangular)
TAB = typeof((/)(zero(eltype(A)), one(eltype(B))) +
(/)(zero(eltype(A)), one(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return UpperTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UnitUpperTriangular, B::UnitUpperTriangular)
TAB = typeof((*)(zero(eltype(A)), zero(eltype(B))) +
(*)(zero(eltype(A)), zero(eltype(B))))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
return UnitUpperTriangular(rdiv!(AA, convert(AbstractMatrix{TAB}, B)))
end

_inner_type_promotion(A,B) = promote_type(eltype(A), eltype(B), typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))))
## The general promotion methods
function *(A::AbstractTriangular, B::AbstractTriangular)
TAB = _inner_type_promotion(A,B)
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
lmul!(convert(AbstractArray{TAB}, A), BB)
end

Expand All @@ -1640,40 +1625,35 @@ for mat in (:AbstractVector, :AbstractMatrix)
@eval function *(A::AbstractTriangular, B::$mat)
require_one_based_indexing(B)
TAB = _inner_type_promotion(A,B)
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
lmul!(convert(AbstractArray{TAB}, A), BB)
end
### Left division with triangle to the left hence rhs cannot be transposed. No quotients.
@eval function \(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::$mat)
require_one_based_indexing(B)
TAB = _inner_type_promotion(A,B)
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
ldiv!(convert(AbstractArray{TAB}, A), BB)
end
### Left division with triangle to the left hence rhs cannot be transposed. Quotients.
@eval function \(A::Union{UpperTriangular,LowerTriangular}, B::$mat)
require_one_based_indexing(B)
TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A)))
BB = similar(B, TAB, size(B))
copyto!(BB, B)
BB = copy_similar(B, TAB)
ldiv!(convert(AbstractArray{TAB}, A), BB)
end
### Right division with triangle to the right hence lhs cannot be transposed. No quotients.
@eval function /(A::$mat, B::Union{UnitUpperTriangular, UnitLowerTriangular})
require_one_based_indexing(A)
TAB = _inner_type_promotion(A,B)
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
rdiv!(AA, convert(AbstractArray{TAB}, B))
end
### Right division with triangle to the right hence lhs cannot be transposed. Quotients.
@eval function /(A::$mat, B::Union{UpperTriangular,LowerTriangular})
require_one_based_indexing(A)
TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A)))
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
rdiv!(AA, convert(AbstractArray{TAB}, B))
end
end
Expand All @@ -1682,8 +1662,7 @@ end
function *(A::AbstractMatrix, B::AbstractTriangular)
require_one_based_indexing(A)
TAB = _inner_type_promotion(A,B)
AA = similar(A, TAB, size(A))
copyto!(AA, A)
AA = copy_similar(A, TAB)
rmul!(AA, convert(AbstractArray{TAB}, B))
end
# ambiguity resolution with definitions in linalg/rowvector.jl
Expand Down
Loading