diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 8dfed4b5c89b9..ebdb78a0c63c2 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -354,8 +354,58 @@ 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) + + +""" + copy_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. + +See also: `copy_similar`, `copy_to_array`. +""" +copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A,T), 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. For example, +supplying a tridiagonal matrix results in a sparse array. In general, the type +of the output corresponds to that of the three-argument method `similar(A, T, size(s))`. + +See also: `copy_oftype`, `copy_to_array`. +""" +copy_similar(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T, size(A)), A) + +""" + copy_to_array(A, T) + +Copy `A` to a regular dense `Array` with element type `T`. + +The resulting array is mutable. It can be used, for example, to pass the data of +`A` to an efficient in-place method for a matrix factorization such as `lu!`, in +cases where a more specific implementation of `lu!` (or `lu`) is not available. + +See also: `copy_oftype`, `copy_similar` +""" +copy_to_array(A::AbstractArray, ::Type{T}) where {T} = copyto!(Array{T}(undef, 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") diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index dac216f0f072a..29f91e6edaed4 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -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) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 66f9ab6a96159..461d38c1e8759 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -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) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 8b1ca1694607e..07b96c00d54ca 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -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) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index b651e85512f6d..626a1ae7b1a74 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -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)) diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl index b8a34712dc9ae..1a71b0604b5a2 100644 --- a/stdlib/LinearAlgebra/src/givens.jl +++ b/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), 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) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index e3b5d4983b30d..9a864da2b0a37 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -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)) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index e55508db6b397..22ed1fc5a49dc 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -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) @@ -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) @@ -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) @@ -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 diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index d0ec430347193..d467300a4e660 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -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 @@ -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)) diff --git a/stdlib/LinearAlgebra/src/schur.jl b/stdlib/LinearAlgebra/src/schur.jl index 268758383b7f1..610067fe51452 100644 --- a/stdlib/LinearAlgebra/src/schur.jl +++ b/stdlib/LinearAlgebra/src/schur.jl @@ -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) diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index c141f6469e359..31cf04a2e5fae 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -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) @@ -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 @@ -1572,57 +1566,49 @@ 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 @@ -1630,8 +1616,7 @@ _inner_type_promotion(A,B) = promote_type(eltype(A), eltype(B), typeof(zero(elty ## 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 @@ -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 @@ -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 diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index 99de21d4bd76d..aa2df6bb41516 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -159,6 +159,9 @@ similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal(similar(S.dv, T # The method below is moved to SparseArrays for now # similar(S::SymTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) +copyto!(dest::SymTridiagonal, src::SymTridiagonal) = + (copyto!(dest.dv, src.dv); copyto!(dest.ev, src.ev); dest) + #Elementary operations for func in (:conj, :copy, :real, :imag) @eval ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev)) diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index f0f0461c23862..8226ddf004a72 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -4,6 +4,8 @@ module TestAdjointTranspose using Test, LinearAlgebra, SparseArrays +const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") + @testset "Adjoint and Transpose inner constructor basics" begin intvec, intmat = [1, 2], [1 2; 3 4] # Adjoint/Transpose eltype must match the type of the Adjoint/Transpose of the input eltype @@ -239,6 +241,25 @@ end @test convert(Transpose{Float64,Matrix{Float64}}, Transpose(intmat))::Transpose{Float64,Matrix{Float64}} == Transpose(intmat) end +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Adjoint and Transpose convert methods to AbstractArray" begin + # tests corresponding to #34995 + intvec, intmat = [1, 2], [1 2 3; 4 5 6] + statvec = ImmutableArray(intvec) + statmat = ImmutableArray(intmat) + + @test convert(AbstractArray{Float64}, Adjoint(statvec))::Adjoint{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Adjoint(statvec) + @test convert(AbstractArray{Float64}, Adjoint(statmat))::Array{Float64,2} == Adjoint(statmat) + @test convert(AbstractArray{Float64}, Transpose(statvec))::Transpose{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Transpose(statvec) + @test convert(AbstractArray{Float64}, Transpose(statmat))::Array{Float64,2} == Transpose(statmat) + @test convert(AbstractMatrix{Float64}, Adjoint(statvec))::Adjoint{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Adjoint(statvec) + @test convert(AbstractMatrix{Float64}, Adjoint(statmat))::Array{Float64,2} == Adjoint(statmat) + @test convert(AbstractMatrix{Float64}, Transpose(statvec))::Transpose{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Transpose(statvec) + @test convert(AbstractMatrix{Float64}, Transpose(statmat))::Array{Float64,2} == Transpose(statmat) +end + @testset "Adjoint and Transpose similar methods" begin intvec, intmat = [1, 2], [1 2 3; 4 5 6] # similar with no additional specifications, vector (rewrapping) semantics @@ -527,7 +548,6 @@ end @test pointer(Transpose(D)) === pointer(D) end -const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl")) using .Main.OffsetArrays diff --git a/stdlib/LinearAlgebra/test/bidiag.jl b/stdlib/LinearAlgebra/test/bidiag.jl index e4dcd14053778..161540a9883b4 100644 --- a/stdlib/LinearAlgebra/test/bidiag.jl +++ b/stdlib/LinearAlgebra/test/bidiag.jl @@ -649,4 +649,20 @@ end @test c \ A ≈ c \ Matrix(A) end +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Conversion to AbstractArray" begin + # tests corresponding to #34995 + dv = ImmutableArray([1, 2, 3, 4]) + ev = ImmutableArray([7, 8, 9]) + Bu = Bidiagonal(dv, ev, :U) + Bl = Bidiagonal(dv, ev, :L) + + @test convert(AbstractArray{Float64}, Bu)::Bidiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Bu + @test convert(AbstractMatrix{Float64}, Bu)::Bidiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Bu + @test convert(AbstractArray{Float64}, Bl)::Bidiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Bl + @test convert(AbstractMatrix{Float64}, Bl)::Bidiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Bl +end + end # module TestBidiagonal diff --git a/stdlib/LinearAlgebra/test/diagonal.jl b/stdlib/LinearAlgebra/test/diagonal.jl index 14d555f7de92b..7c67c8d36dc03 100644 --- a/stdlib/LinearAlgebra/test/diagonal.jl +++ b/stdlib/LinearAlgebra/test/diagonal.jl @@ -807,4 +807,17 @@ end end end +const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Conversion to AbstractArray" begin + # tests corresponding to #34995 + d = ImmutableArray([1, 2, 3, 4]) + D = Diagonal(d) + + @test convert(AbstractArray{Float64}, D)::Diagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == D + @test convert(AbstractMatrix{Float64}, D)::Diagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == D +end + end # module TestDiagonal diff --git a/stdlib/LinearAlgebra/test/hessenberg.jl b/stdlib/LinearAlgebra/test/hessenberg.jl index 0f9246c722349..65dc029060596 100644 --- a/stdlib/LinearAlgebra/test/hessenberg.jl +++ b/stdlib/LinearAlgebra/test/hessenberg.jl @@ -196,4 +196,16 @@ end end end +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Conversion to AbstractArray" begin + # tests corresponding to #34995 + A = ImmutableArray([1 2 3; 4 5 6; 7 8 9]) + H = UpperHessenberg(A) + + @test convert(AbstractArray{Float64}, H)::UpperHessenberg{Float64,ImmutableArray{Float64,2,Array{Float64,2}}} == H + @test convert(AbstractMatrix{Float64}, H)::UpperHessenberg{Float64,ImmutableArray{Float64,2,Array{Float64,2}}} == H +end + end # module TestHessenberg diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index d6ec50fe4348d..0dffe7fa1738f 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -398,4 +398,21 @@ end @test a == c end end + +@testset "lu(A) has a fallback for abstract matrices (#40831)" begin + # check that lu works for some structured arrays + A0 = rand(5, 5) + @test lu(Diagonal(A0)) isa LU + @test Matrix(lu(Diagonal(A0))) ≈ Diagonal(A0) + @test lu(Bidiagonal(A0, :U)) isa LU + @test Matrix(lu(Bidiagonal(A0, :U))) ≈ Bidiagonal(A0, :U) + + # lu(A) copies A and then invokes lu!, make sure that the most efficient + # implementation of lu! continues to be used + A1 = Tridiagonal(rand(2), rand(3), rand(2)) + @test lu(A1) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} + @test lu(A1, RowMaximum()) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} + @test lu(A1, RowMaximum(); check = false) isa LU{Float64, Tridiagonal{Float64, Vector{Float64}}} +end + end # module TestLU diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index bdd003b57913e..f20b6fe2acc97 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -551,6 +551,26 @@ end end end +const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test") +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Conversion to AbstractArray" begin + # tests corresponding to #34995 + immutablemat = ImmutableArray([1 2 3; 4 5 6; 7 8 9]) + for SymType in (Symmetric, Hermitian) + S = Float64 + symmat = SymType(immutablemat) + @test convert(AbstractArray{S}, symmat).data isa ImmutableArray{S} + @test convert(AbstractMatrix{S}, symmat).data isa ImmutableArray{S} + @test AbstractArray{S}(symmat).data isa ImmutableArray{S} + @test AbstractMatrix{S}(symmat).data isa ImmutableArray{S} + @test convert(AbstractArray{S}, symmat) == symmat + @test convert(AbstractMatrix{S}, symmat) == symmat + end +end + + @testset "#24572: eltype(A::HermOrSym) === eltype(parent(A))" begin A = rand(Float32, 3, 3) @test_throws TypeError Symmetric{Float64,Matrix{Float32}}(A, 'U') diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl index 030e4a27625f6..abc166a86e150 100644 --- a/stdlib/LinearAlgebra/test/triangular.jl +++ b/stdlib/LinearAlgebra/test/triangular.jl @@ -662,6 +662,25 @@ end end end +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "AbstractArray constructor should preserve underlying storage type" begin + # tests corresponding to #34995 + local m = 4 + local T, S = Float32, Float64 + immutablemat = ImmutableArray(randn(T,m,m)) + for TriType in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular) + trimat = TriType(immutablemat) + @test convert(AbstractArray{S}, trimat).data isa ImmutableArray{S} + @test convert(AbstractMatrix{S}, trimat).data isa ImmutableArray{S} + @test AbstractArray{S}(trimat).data isa ImmutableArray{S} + @test AbstractMatrix{S}(trimat).data isa ImmutableArray{S} + @test convert(AbstractArray{S}, trimat) == trimat + @test convert(AbstractMatrix{S}, trimat) == trimat + end +end + @testset "inplace mul of appropriate types should preserve triagular structure" begin for elty1 in (Float64, ComplexF32), elty2 in (Float64, ComplexF32) T = promote_type(elty1, elty2) diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index 079c7aa0d5a91..1f8ca042db0b7 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -198,10 +198,8 @@ end @testset "similar, size, and copyto!" begin B = similar(A) @test size(B) == size(A) - if mat_type == Tridiagonal # doesn't work for SymTridiagonal yet - copyto!(B, A) - @test B == A - end + copyto!(B, A) + @test B == A @test isa(similar(A), mat_type{elty}) @test isa(similar(A, Int), mat_type{Int}) @test isa(similar(A, (3, 2)), SparseMatrixCSC) @@ -651,4 +649,21 @@ end @test ishermitian(S) end +isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl")) +using .Main.ImmutableArrays + +@testset "Conversion to AbstractArray" begin + # tests corresponding to #34995 + v1 = ImmutableArray([1, 2]) + v2 = ImmutableArray([3, 4, 5]) + v3 = ImmutableArray([6, 7]) + T = Tridiagonal(v1, v2, v3) + Tsym = SymTridiagonal(v2, v1) + + @test convert(AbstractArray{Float64}, T)::Tridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == T + @test convert(AbstractMatrix{Float64}, T)::Tridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == T + @test convert(AbstractArray{Float64}, Tsym)::SymTridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Tsym + @test convert(AbstractMatrix{Float64}, Tsym)::SymTridiagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == Tsym +end + end # module TestTridiagonal diff --git a/test/testhelpers/ImmutableArrays.jl b/test/testhelpers/ImmutableArrays.jl new file mode 100644 index 0000000000000..df2a78387e07b --- /dev/null +++ b/test/testhelpers/ImmutableArrays.jl @@ -0,0 +1,28 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +# ImmutableArrays (arrays that implement getindex but not setindex!) + +# This test file defines an array wrapper that is immutable. It can be used to +# test the action of methods on immutable arrays. + +module ImmutableArrays + +export ImmutableArray + +"An immutable wrapper type for arrays." +struct ImmutableArray{T,N,A<:AbstractArray} <: AbstractArray{T,N} + data::A +end + +ImmutableArray(data::AbstractArray{T,N}) where {T,N} = ImmutableArray{T,N,typeof(data)}(data) + +# Minimal AbstractArray interface +Base.size(A::ImmutableArray) = size(A.data) +Base.size(A::ImmutableArray, d) = size(A.data, d) +Base.getindex(A::ImmutableArray, i...) = getindex(A.data, i...) + +# The immutable array remains immutable after conversion to AbstractArray +AbstractArray{T}(A::ImmutableArray) where {T} = ImmutableArray(AbstractArray{T}(A.data)) +AbstractArray{T,N}(A::ImmutableArray{S,N}) where {S,T,N} = ImmutableArray(AbstractArray{T,N}(A.data)) + +end