From 78a9d47de96531db4ca19ca154fd3953131852b2 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 23 Sep 2022 19:03:21 +0200 Subject: [PATCH 01/12] Introduce `AdjointFactorization` --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 6 +- stdlib/LinearAlgebra/src/adjtrans.jl | 12 +-- stdlib/LinearAlgebra/src/factorization.jl | 90 ++++++++++++++-------- stdlib/LinearAlgebra/src/hessenberg.jl | 8 +- stdlib/LinearAlgebra/src/lq.jl | 5 +- stdlib/LinearAlgebra/src/lu.jl | 49 ++++-------- stdlib/LinearAlgebra/src/qr.jl | 4 +- stdlib/LinearAlgebra/test/adjtrans.jl | 8 +- stdlib/LinearAlgebra/test/factorization.jl | 21 ++++- stdlib/LinearAlgebra/test/qr.jl | 6 +- stdlib/LinearAlgebra/test/testgroups | 2 +- 11 files changed, 115 insertions(+), 96 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index b107777e7ae85..4c6bbfffbf9cd 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -493,7 +493,7 @@ const LAPACKFactorizations{T,S} = Union{ QRCompactWY{T,S}, QRPivoted{T,S}, SVD{T,<:Real,S}} -function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorizations}}, B::AbstractVecOrMat) +function (\)(F::Union{<:LAPACKFactorizations,AdjointFactorization{<:Any,<:LAPACKFactorizations}}, B::AbstractVecOrMat) require_one_based_indexing(B) m, n = size(F) if m != size(B, 1) @@ -523,7 +523,9 @@ function (\)(F::Union{<:LAPACKFactorizations,Adjoint{<:Any,<:LAPACKFactorization end # disambiguate (\)(F::LAPACKFactorizations{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = - invoke(\, Tuple{Factorization{T}, VecOrMat{Complex{T}}}, F, B) + @invoke \(F::Factorization{T}, B::VecOrMat{Complex{T}}) +(\)(F::AdjointFactorization{T,<:LAPACKFactorizations}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = + @invoke \(F::AdjointFactorization{<:Any,<:LAPACKFactorizations}, B::AbstractVecOrMat) """ LinearAlgebra.peakflops(n::Integer=2000; parallel::Bool=false) diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index 1a67b7f69e24a..46133c6fd9fd4 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -32,7 +32,7 @@ julia> Adjoint(A) 9-2im 0+0im ``` """ -struct Adjoint{T,S} <: AbstractMatrix{T} +struct Adjoint{T,S<:AbstractVecOrMat} <: AbstractMatrix{T} parent::S end """ @@ -59,13 +59,13 @@ julia> Transpose(A) 3 0 ``` """ -struct Transpose{T,S} <: AbstractMatrix{T} +struct Transpose{T,S<:AbstractVecOrMat} <: AbstractMatrix{T} parent::S end # basic outer constructors -Adjoint(A) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A) -Transpose(A) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A) +Adjoint(A::AbstractVecOrMat) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A) +Transpose(A::AbstractVecOrMat) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A) Base.dataids(A::Union{Adjoint, Transpose}) = Base.dataids(A.parent) Base.unaliascopy(A::Union{Adjoint,Transpose}) = typeof(A)(Base.unaliascopy(A.parent)) @@ -460,8 +460,8 @@ pinv(v::TransposeAbsVec, tol::Real = 0) = pinv(conj(v.parent)).parent ## right-division / /(u::AdjointAbsVec, A::AbstractMatrix) = adjoint(adjoint(A) \ u.parent) /(u::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A) \ u.parent) -/(u::AdjointAbsVec, A::Transpose{<:Any,<:AbstractMatrix}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) -/(u::TransposeAbsVec, A::Adjoint{<:Any,<:AbstractMatrix}) = transpose(conj(A.parent) \ u.parent) # technically should be transpose(copy(transpose(copy(A))) \ u.parent) +/(u::AdjointAbsVec, A::TransposeAbsMat) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent) +/(u::TransposeAbsVec, A::AdjointAbsMat) = transpose(conj(A.parent) \ u.parent) # technically should be transpose(copy(transpose(copy(A))) \ u.parent) ## complex conjugate conj(A::Transpose) = adjoint(A.parent) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index a5bcdf66ecc7c..58832f9b2a06b 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -11,9 +11,32 @@ matrix factorizations. """ abstract type Factorization{T} end +struct AdjointFactorization{T,S<:Factorization} <: Factorization{T} + parent::S +end +AdjointFactorization(F::Factorization) = + AdjointFactorization{Base.promote_op(adjoint,eltype(F)),typeof(F)}(F) + +struct TransposeFactorization{T,S<:Factorization} <: Factorization{T} + parent::S +end +TransposeFactorization(F::Factorization) = + TransposeFactorization{Base.promote_op(adjoint,eltype(F)),typeof(F)}(F) + eltype(::Type{<:Factorization{T}}) where {T} = T -size(F::Adjoint{<:Any,<:Factorization}) = reverse(size(parent(F))) -size(F::Transpose{<:Any,<:Factorization}) = reverse(size(parent(F))) +size(F::AdjointFactorization) = reverse(size(parent(F))) +size(F::TransposeFactorization) = reverse(size(parent(F))) +size(F::Union{AdjointFactorization,TransposeFactorization}, d::Integer) = d in (1, 2) ? size(F)[d] : 1 +parent(F::Union{AdjointFactorization,TransposeFactorization}) = F.parent + +adjoint(F::Factorization) = AdjointFactorization(F) +transpose(F::Factorization) = TransposeFactorization(F) +transpose(F::Factorization{<:Real}) = AdjointFactorization(F) +adjoint(F::AdjointFactorization) = F.parent +transpose(F::TransposeFactorization) = F.parent +transpose(F::AdjointFactorization{<:Real}) = F.parent +conj(A::TransposeFactorization) = adjoint(A.parent) +conj(A::AdjointFactorization) = transpose(A.parent) checkpositivedefinite(info) = info == 0 || throw(PosDefException(info)) checknonsingular(info, ::RowMaximum) = info == 0 || throw(SingularException(info)) @@ -60,8 +83,8 @@ convert(::Type{T}, f::Factorization) where {T<:AbstractArray} = T(f)::T ### General promotion rules Factorization{T}(F::Factorization{T}) where {T} = F -# This is a bit odd since the return is not a Factorization but it works well in generic code -Factorization{T}(A::Adjoint{<:Any,<:Factorization}) where {T} = +# This no longer looks odd since the return _is_ a Factorization! +Factorization{T}(A::AdjointFactorization) where {T} = adjoint(Factorization{T}(parent(A))) inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n))) @@ -69,55 +92,67 @@ Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash Base.:(==)( F::T, G::T) where {T<:Factorization} = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F)) Base.isequal(F::T, G::T) where {T<:Factorization} = all(f -> isequal(getfield(F, f), getfield(G, f)), 1:nfields(F))::Bool -function Base.show(io::IO, x::Adjoint{<:Any,<:Factorization}) - print(io, "Adjoint of ") +function Base.show(io::IO, x::AdjointFactorization) + print(io, "adjoint of ") show(io, parent(x)) end -function Base.show(io::IO, x::Transpose{<:Any,<:Factorization}) - print(io, "Transpose of ") +function Base.show(io::IO, x::TransposeFactorization) + print(io, "transpose of ") show(io, parent(x)) end -function Base.show(io::IO, ::MIME"text/plain", x::Adjoint{<:Any,<:Factorization}) - print(io, "Adjoint of ") +function Base.show(io::IO, ::MIME"text/plain", x::AdjointFactorization) + print(io, "adjoint of ") show(io, MIME"text/plain"(), parent(x)) end -function Base.show(io::IO, ::MIME"text/plain", x::Transpose{<:Any,<:Factorization}) - print(io, "Transpose of ") +function Base.show(io::IO, ::MIME"text/plain", x::TransposeFactorization) + print(io, "transpose of ") show(io, MIME"text/plain"(), parent(x)) end # With a real lhs and complex rhs with the same precision, we can reinterpret # the complex rhs as a real rhs with twice the number of columns or rows -function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal +function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} require_one_based_indexing(B) c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2)) x = ldiv!(F, c2r) return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))), _ret_size(F, B)) end -function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where T<:BlasReal +# don't do the reinterpretation for [Adjoint/Transpose]Factorization +(\)(F::TransposeFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = + conj!(adjoint(parent(F)) \ conj.(B)) +(\)(F::AdjointFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = + @invoke \(F::typeof(F), B::AbstractVecOrMat) # send to meta-function in LinearAlgebra.jl + +function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where {T<:BlasReal} require_one_based_indexing(B) x = rdiv!(copy(reinterpret(T, B)), F) return copy(reinterpret(Complex{T}, x)) end +# don't do the reinterpretation for [Adjoint/Transpose]Factorization +(/)(B::VecOrMat{Complex{T}}, F::TransposeFactorization{T}) where {T<:BlasReal} = + conj!(adjoint(parent(F)) \ conj.(B)) +(/)(B::VecOrMat{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = + invoke(/, Tuple{VecOrMat{Complex{T}}, Factorization{T}}, B, F) -function \(F::Union{Factorization, Adjoint{<:Any,<:Factorization}}, B::AbstractVecOrMat) +function (\)(F::Factorization, B::AbstractVecOrMat) require_one_based_indexing(B) - TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F))) + TFB = typeof(oneunit(eltype(F)) / oneunit(eltype(B))) ldiv!(F, copy_similar(B, TFB)) end +(\)(F::TransposeFactorization{<:Real}, B::AbstractVecOrMat) = adjoint(F.parent) \ B +(\)(F::TransposeFactorization, B::AbstractVecOrMat) = conj!(adjoint(F.parent) \ conj.(B)) -function /(B::AbstractMatrix, F::Union{Factorization, Adjoint{<:Any,<:Factorization}}) +function (/)(B::AbstractMatrix, F::Factorization) require_one_based_indexing(B) TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(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)) - +(/)(A::AbstractMatrix, F::AdjointFactorization) = adjoint(adjoint(F) \ adjoint(A)) +(/)(A::AbstractMatrix, F::TransposeFactorization) = transpose(transpose(F) \ transpose(A)) function ldiv!(Y::AbstractVector, A::Factorization, B::AbstractVector) require_one_based_indexing(Y, B) - m, n = size(A, 1), size(A, 2) + m, n = size(A) if m > n Bc = copy(B) ldiv!(A, Bc) @@ -128,7 +163,7 @@ function ldiv!(Y::AbstractVector, A::Factorization, B::AbstractVector) end function ldiv!(Y::AbstractMatrix, A::Factorization, B::AbstractMatrix) require_one_based_indexing(Y, B) - m, n = size(A, 1), size(A, 2) + m, n = size(A) if m > n Bc = copy(B) ldiv!(A, Bc) @@ -138,14 +173,3 @@ function ldiv!(Y::AbstractMatrix, A::Factorization, B::AbstractMatrix) return ldiv!(A, Y) end end - -# fallback methods for transposed solves -\(F::Transpose{<:Any,<:Factorization{<:Real}}, B::AbstractVecOrMat) = adjoint(F.parent) \ B -\(F::Transpose{<:Any,<:Factorization}, B::AbstractVecOrMat) = conj.(adjoint(F.parent) \ conj.(B)) - -/(B::AbstractMatrix, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent) -/(B::AbstractMatrix, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent)) -/(B::AdjointAbsVec, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent) -/(B::TransposeAbsVec, F::Transpose{<:Any,<:Factorization{<:Real}}) = B / adjoint(F.parent) -/(B::AdjointAbsVec, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent)) -/(B::TransposeAbsVec, F::Transpose{<:Any,<:Factorization}) = conj.(conj.(B) / adjoint(F.parent)) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index a95a73dfc8819..81b0c001141e8 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -421,11 +421,9 @@ Hessenberg(F::Hessenberg, μ::Number) = Hessenberg(F.factors, F.τ, F.H, F.uplo; copy(F::Hessenberg{<:Any,<:UpperHessenberg}) = Hessenberg(copy(F.factors), copy(F.τ); μ=F.μ) copy(F::Hessenberg{<:Any,<:SymTridiagonal}) = Hessenberg(copy(F.factors), copy(F.τ), copy(F.H), F.uplo; μ=F.μ) -size(F::Hessenberg, d) = size(F.H, d) +size(F::Hessenberg, d::Integer) = size(F.H, d) size(F::Hessenberg) = size(F.H) -adjoint(F::Hessenberg) = Adjoint(F) - # iteration for destructuring into components Base.iterate(S::Hessenberg) = (S.Q, Val(:H)) Base.iterate(S::Hessenberg, ::Val{:H}) = (S.H, Val(:μ)) @@ -686,8 +684,8 @@ function rdiv!(B::AbstractVecOrMat{<:Complex}, F::Hessenberg{<:Complex,<:Any,<:A return B .= Complex.(Br,Bi) end -ldiv!(F::Adjoint{<:Any,<:Hessenberg}, B::AbstractVecOrMat) = rdiv!(B', F')' -rdiv!(B::AbstractMatrix, F::Adjoint{<:Any,<:Hessenberg}) = ldiv!(F', B')' +ldiv!(F::AdjointFactorization{<:Any,<:Hessenberg}, B::AbstractVecOrMat) = rdiv!(B', F')' +rdiv!(B::AbstractMatrix, F::AdjointFactorization{<:Any,<:Hessenberg}) = ldiv!(F', B')' det(F::Hessenberg) = det(F.H; shift=F.μ) logabsdet(F::Hessenberg) = logabsdet(F.H; shift=F.μ) diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index 52d4f944f682f..86b8a9674c69f 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -135,8 +135,7 @@ AbstractArray(A::LQ) = AbstractMatrix(A) Matrix(A::LQ) = Array(AbstractArray(A)) Array(A::LQ) = Matrix(A) -adjoint(A::LQ) = Adjoint(A) -Base.copy(F::Adjoint{T,<:LQ{T}}) where {T} = +Base.copy(F::AdjointFactorization{T,<:LQ{T}}) where {T} = QR{T,typeof(F.parent.factors),typeof(F.parent.τ)}(copy(adjoint(F.parent.factors)), copy(F.parent.τ)) function getproperty(F::LQ, d::Symbol) @@ -343,7 +342,7 @@ function ldiv!(A::LQ, B::StridedVecOrMat) return lmul!(adjoint(A.Q), B) end -function ldiv!(Fadj::Adjoint{<:Any,<:LQ}, B::StridedVecOrMat) +function ldiv!(Fadj::AdjointFactorization{<:Any,<:LQ}, B::StridedVecOrMat) require_one_based_indexing(B) m, n = size(Fadj) m >= n || throw(DimensionMismatch("solver does not support underdetermined systems (more columns than rows)")) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 47e3fbfcb0232..da3281946bbf4 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -72,8 +72,11 @@ Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p)) Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done)) Base.iterate(S::LU, ::Val{:done}) = nothing -adjoint(F::LU) = Adjoint(F) -transpose(F::LU) = Transpose(F) +# LU prefers transpose over adjoint in the real case +adjoint(F::LU) = AdjointFactorization(F) +transpose(F::LU) = TransposeFactorization(F) +# override the generic fallback +transpose(F::LU{<:Real}) = TransposeFactorization(F) # StridedMatrix lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMaximum(); check=check) @@ -327,7 +330,7 @@ Factorization{T}(F::LU) where {T} = LU{T}(F) copy(A::LU{T,S,P}) where {T,S,P} = LU{T,S,P}(copy(A.factors), copy(A.ipiv), A.info) size(A::LU) = size(getfield(A, :factors)) -size(A::LU, i) = size(getfield(A, :factors), i) +size(A::LU, i::Integer) = size(getfield(A, :factors), i) function ipiv2perm(v::AbstractVector{T}, maxi::Integer) where T require_one_based_indexing(v) @@ -432,49 +435,29 @@ function ldiv!(A::LU{<:Any,<:StridedMatrix}, B::StridedVecOrMat) ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B)) end -ldiv!(transA::Transpose{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = +ldiv!(transA::TransposeFactorization{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = (A = transA.parent; LAPACK.getrs!('T', A.factors, A.ipiv, B)) -function ldiv!(transA::Transpose{<:Any,<:LU{<:Any,<:StridedMatrix}}, B::StridedVecOrMat) +function ldiv!(transA::TransposeFactorization{<:Any,<:LU{<:Any,<:StridedMatrix}}, B::StridedVecOrMat) A = transA.parent ldiv!(transpose(UnitLowerTriangular(A.factors)), ldiv!(transpose(UpperTriangular(A.factors)), B)) _apply_inverse_ipiv_rows!(A, B) end -ldiv!(adjF::Adjoint{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:Real} = - (F = adjF.parent; ldiv!(transpose(F), B)) -ldiv!(adjA::Adjoint{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = +ldiv!(adjA::AdjointFactorization{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = (A = adjA.parent; LAPACK.getrs!('C', A.factors, A.ipiv, B)) -function ldiv!(adjA::Adjoint{<:Any,<:LU{<:Any,<:StridedMatrix}}, B::StridedVecOrMat) +function ldiv!(adjA::AdjointFactorization{<:Any,<:LU{<:Any,<:StridedMatrix}}, B::StridedVecOrMat) A = adjA.parent ldiv!(adjoint(UnitLowerTriangular(A.factors)), ldiv!(adjoint(UpperTriangular(A.factors)), B)) _apply_inverse_ipiv_rows!(A, B) end -(\)(A::Adjoint{<:Any,<:LU}, B::Adjoint{<:Any,<:StridedVecOrMat}) = A \ copy(B) -(\)(A::Transpose{<:Any,<:LU}, B::Transpose{<:Any,<:StridedVecOrMat}) = A \ copy(B) -(\)(A::Adjoint{T,<:LU{T,<:StridedMatrix}}, B::Adjoint{T,<:StridedVecOrMat{T}}) where {T<:BlasComplex} = +(\)(A::AdjointFactorization{T,<:LU{T,<:StridedMatrix}}, B::Adjoint{T,<:StridedVecOrMat{T}}) where {T<:BlasComplex} = LAPACK.getrs!('C', A.parent.factors, A.parent.ipiv, copy(B)) -(\)(A::Transpose{T,<:LU{T,<:StridedMatrix}}, B::Transpose{T,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = +(\)(A::TransposeFactorization{T,<:LU{T,<:StridedMatrix}}, B::Transpose{T,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = LAPACK.getrs!('T', A.parent.factors, A.parent.ipiv, copy(B)) -function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:LU}) - T = promote_type(eltype(A), eltype(F)) - 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!(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!(copymutable_oftype(trA.parent, T)))) -end - function det(F::LU{T}) where T n = checksquare(F) issuccess(F) || return zero(T) @@ -650,7 +633,7 @@ function ldiv!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V} return B end -function ldiv!(transA::Transpose{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} +function ldiv!(transA::TransposeFactorization{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} require_one_based_indexing(B) A = transA.parent n = size(A,1) @@ -687,7 +670,7 @@ function ldiv!(transA::Transpose{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVec end # Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B) -function ldiv!(adjA::Adjoint{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} +function ldiv!(adjA::AdjointFactorization{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} require_one_based_indexing(B) A = adjA.parent n = size(A,1) @@ -724,8 +707,8 @@ function ldiv!(adjA::Adjoint{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMa end rdiv!(B::AbstractMatrix, A::LU) = transpose(ldiv!(transpose(A), transpose(B))) -rdiv!(B::AbstractMatrix, A::Transpose{<:Any,<:LU}) = transpose(ldiv!(A.parent, transpose(B))) -rdiv!(B::AbstractMatrix, A::Adjoint{<:Any,<:LU}) = adjoint(ldiv!(A.parent, adjoint(B))) +rdiv!(B::AbstractMatrix, A::TransposeFactorization{<:Any,<:LU}) = transpose(ldiv!(A.parent, transpose(B))) +rdiv!(B::AbstractMatrix, A::AdjointFactorization{<:Any,<:LU}) = adjoint(ldiv!(A.parent, adjoint(B))) # Conversions AbstractMatrix(F::LU) = (F.L * F.U)[invperm(F.p),:] diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 741b53bcd56a9..f7a290464a205 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -514,7 +514,7 @@ end Base.propertynames(F::QRPivoted, private::Bool=false) = (:R, :Q, :p, :P, (private ? fieldnames(typeof(F)) : ())...) -adjoint(F::Union{QR,QRPivoted,QRCompactWY}) = Adjoint(F) +transpose(F::Union{QR,QRPivoted,QRCompactWY}) = throw(ArgumentError("transpose of QR decomposition is not supported")) abstract type AbstractQ{T} <: AbstractMatrix{T} end @@ -1001,7 +1001,7 @@ function _apply_permutation!(F::QRPivoted, B::AbstractVecOrMat) end _apply_permutation!(F::Factorization, B::AbstractVecOrMat) = B -function ldiv!(Fadj::Adjoint{<:Any,<:Union{QR,QRCompactWY,QRPivoted}}, B::AbstractVecOrMat) +function ldiv!(Fadj::AdjointFactorization{<:Any,<:Union{QR,QRCompactWY,QRPivoted}}, B::AbstractVecOrMat) require_one_based_indexing(B) m, n = size(Fadj) diff --git a/stdlib/LinearAlgebra/test/adjtrans.jl b/stdlib/LinearAlgebra/test/adjtrans.jl index e96ea28531d37..cb713bd86340e 100644 --- a/stdlib/LinearAlgebra/test/adjtrans.jl +++ b/stdlib/LinearAlgebra/test/adjtrans.jl @@ -489,13 +489,13 @@ end @test B == A .* A' end -@testset "test show methods for $t of Factorizations" for t in (Adjoint, Transpose) - A = randn(4, 4) +@testset "test show methods for $t of Factorizations" for t in (adjoint, transpose) + A = randn(ComplexF64, 4, 4) F = lu(A) Fop = t(F) - @test "LinearAlgebra."*sprint(show, Fop) == + @test sprint(show, Fop) == "$t of "*sprint(show, parent(Fop)) - @test "LinearAlgebra."*sprint((io, t) -> show(io, MIME"text/plain"(), t), Fop) == + @test sprint((io, t) -> show(io, MIME"text/plain"(), t), Fop) == "$t of "*sprint((io, t) -> show(io, MIME"text/plain"(), t), parent(Fop)) end diff --git a/stdlib/LinearAlgebra/test/factorization.jl b/stdlib/LinearAlgebra/test/factorization.jl index d200eff2f17bf..72233293ff515 100644 --- a/stdlib/LinearAlgebra/test/factorization.jl +++ b/stdlib/LinearAlgebra/test/factorization.jl @@ -56,11 +56,24 @@ end A = randn(3, 3) A = A * A' # ensure A is pos. def. and symmetric F = f(A) - tF = Transpose(F) - aF = Adjoint(F) @test size(F) == size(A) - @test size(tF) == size(Transpose(A)) - @test size(aF) == size(Adjoint(A)) + @test size(F') == size(A') +end + +@testset "size for transpose factorizations - $f" for f in Any[ + bunchkaufman, + cholesky, + x -> cholesky(x, RowMaximum()), + hessenberg, + lq, + lu, + svd, +] + A = randn(3, 3) + A = A * A' # ensure A is pos. def. and symmetric + F = f(A) + @test size(F) == size(A) + @test size(transpose(F)) == size(transpose(A)) end @testset "equality of QRCompactWY" begin diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index c8db95b8c34b6..b78e56ac0ff1b 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -211,9 +211,9 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) end @testset "transpose errors" begin - @test_throws MethodError transpose(qr(randn(3,3))) - @test_throws MethodError transpose(qr(randn(3,3), NoPivot())) - @test_throws MethodError transpose(qr(big.(randn(3,3)))) + @test_throws ArgumentError transpose(qr(randn(3,3))) + @test_throws ArgumentError transpose(qr(randn(3,3), NoPivot())) + @test_throws ArgumentError transpose(qr(big.(randn(3,3)))) end @testset "Issue 7304" begin diff --git a/stdlib/LinearAlgebra/test/testgroups b/stdlib/LinearAlgebra/test/testgroups index de082d8e7dce0..82aff7615da28 100644 --- a/stdlib/LinearAlgebra/test/testgroups +++ b/stdlib/LinearAlgebra/test/testgroups @@ -1,3 +1,4 @@ +addmul triangular qr dense @@ -23,6 +24,5 @@ adjtrans pinv givens structuredbroadcast -addmul ldlt factorization From ffd8c797f9c2bc4a6cac84d6d17dc6d63a6eec48 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 24 Sep 2022 13:43:33 +0200 Subject: [PATCH 02/12] for the time being, don't restrict to matrix --- stdlib/LinearAlgebra/src/adjtrans.jl | 8 ++++---- stdlib/LinearAlgebra/src/lu.jl | 3 +-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index 46133c6fd9fd4..ef815b3ad708b 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -32,7 +32,7 @@ julia> Adjoint(A) 9-2im 0+0im ``` """ -struct Adjoint{T,S<:AbstractVecOrMat} <: AbstractMatrix{T} +struct Adjoint{T,S} <: AbstractMatrix{T} parent::S end """ @@ -59,13 +59,13 @@ julia> Transpose(A) 3 0 ``` """ -struct Transpose{T,S<:AbstractVecOrMat} <: AbstractMatrix{T} +struct Transpose{T,S} <: AbstractMatrix{T} parent::S end # basic outer constructors -Adjoint(A::AbstractVecOrMat) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A) -Transpose(A::AbstractVecOrMat) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A) +Adjoint(A) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A) +Transpose(A) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A) Base.dataids(A::Union{Adjoint, Transpose}) = Base.dataids(A.parent) Base.unaliascopy(A::Union{Adjoint,Transpose}) = typeof(A)(Base.unaliascopy(A.parent)) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index da3281946bbf4..0ad9ee2e385ea 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -72,10 +72,9 @@ Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p)) Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done)) Base.iterate(S::LU, ::Val{:done}) = nothing -# LU prefers transpose over adjoint in the real case adjoint(F::LU) = AdjointFactorization(F) transpose(F::LU) = TransposeFactorization(F) -# override the generic fallback +# LU prefers transpose over adjoint in the real case, override the generic fallback transpose(F::LU{<:Real}) = TransposeFactorization(F) # StridedMatrix From 297bcf2464fabac8d21ec4f4880b807b0b07955f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 25 Sep 2022 22:15:15 +0200 Subject: [PATCH 03/12] provide generic fallbacks --- stdlib/LinearAlgebra/src/adjtrans.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index ef815b3ad708b..331b875aa6422 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -291,6 +291,10 @@ wrapperop(_) = identity wrapperop(::Adjoint) = adjoint wrapperop(::Transpose) = transpose +# the following fallbacks can be removed if Adjoint/Transpose are restricted to AbstractVecOrMat +length(A::AdjOrTrans) = length(A.parent) +size(A::AdjOrTrans) = reverse(size(A.parent)) +axes(A::AdjOrTrans) = reverse(axes(A.parent)) # AbstractArray interface, basic definitions length(A::AdjOrTrans) = length(A.parent) size(v::AdjOrTransAbsVec) = (1, length(v.parent)) From e679db8f013faddf7788ae2f5aafc45672233070 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 9 Jan 2023 22:06:18 +0100 Subject: [PATCH 04/12] clarity --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 13 ++++++++++--- stdlib/LinearAlgebra/src/adjtrans.jl | 1 - stdlib/LinearAlgebra/src/factorization.jl | 3 ++- stdlib/LinearAlgebra/src/lu.jl | 3 +-- 4 files changed, 13 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 5f240f2cb9dbb..61ffd47279dc9 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -501,7 +501,7 @@ _initarray(op, ::Type{TA}, ::Type{TB}, C) where {TA,TB} = # While this definition is pretty general, it does e.g. promote to common element type of lhs and rhs # which is required by LAPACK but not SuiteSparse which allows real-complex solves in some cases. Hence, # we restrict this method to only the LAPACK factorizations in LinearAlgebra. -# The definition is put here since it explicitly references all the Factorizion structs so it has +# The definition is put here since it explicitly references all the Factorization structs so it has # to be located after all the files that define the structs. const LAPACKFactorizations{T,S} = Union{ BunchKaufman{T,S}, @@ -512,7 +512,12 @@ const LAPACKFactorizations{T,S} = Union{ QRCompactWY{T,S}, QRPivoted{T,S}, SVD{T,<:Real,S}} -function (\)(F::Union{<:LAPACKFactorizations,AdjointFactorization{<:Any,<:LAPACKFactorizations}}, B::AbstractVecOrMat) + +(\)(F::LAPACKFactorizations, B::AbstractVecOrMat) = ldiv(F, B) +(\)(F::AdjointFactorization{<:Any,<:LAPACKFactorizations}, B::AbstractVecOrMat) = ldiv(F, B) +(\)(F::TransposeFactorization{<:Any,<:LU}, B::AbstractVecOrMat) = ldiv(F, B) + +function ldiv(F::Factorization, B::AbstractVecOrMat) require_one_based_indexing(B) m, n = size(F) if m != size(B, 1) @@ -544,7 +549,9 @@ end (\)(F::LAPACKFactorizations{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = @invoke \(F::Factorization{T}, B::VecOrMat{Complex{T}}) (\)(F::AdjointFactorization{T,<:LAPACKFactorizations}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = - @invoke \(F::AdjointFactorization{<:Any,<:LAPACKFactorizations}, B::AbstractVecOrMat) + ldiv(F, B) +(\)(F::TransposeFactorization{T,<:LU}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = + ldiv(F, B) """ LinearAlgebra.peakflops(n::Integer=2000; parallel::Bool=false) diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index 0ebf528e35cda..9a872497fdbae 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -292,7 +292,6 @@ wrapperop(::Adjoint) = adjoint wrapperop(::Transpose) = transpose # the following fallbacks can be removed if Adjoint/Transpose are restricted to AbstractVecOrMat -length(A::AdjOrTrans) = length(A.parent) size(A::AdjOrTrans) = reverse(size(A.parent)) axes(A::AdjOrTrans) = reverse(axes(A.parent)) # AbstractArray interface, basic definitions diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 58832f9b2a06b..4f7e316444a25 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -86,6 +86,8 @@ Factorization{T}(F::Factorization{T}) where {T} = F # This no longer looks odd since the return _is_ a Factorization! Factorization{T}(A::AdjointFactorization) where {T} = adjoint(Factorization{T}(parent(A))) +Factorization{T}(A::TransposeFactorization) where {T} = + transpose(Factorization{T}(parent(A))) inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n))) Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash, 1:nfields(F); init=h) @@ -139,7 +141,6 @@ function (\)(F::Factorization, B::AbstractVecOrMat) TFB = typeof(oneunit(eltype(F)) / oneunit(eltype(B))) ldiv!(F, copy_similar(B, TFB)) end -(\)(F::TransposeFactorization{<:Real}, B::AbstractVecOrMat) = adjoint(F.parent) \ B (\)(F::TransposeFactorization, B::AbstractVecOrMat) = conj!(adjoint(F.parent) \ conj.(B)) function (/)(B::AbstractMatrix, F::Factorization) diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 73bfe608824fe..a93803ca2ea45 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -72,9 +72,8 @@ Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p)) Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done)) Base.iterate(S::LU, ::Val{:done}) = nothing -adjoint(F::LU) = AdjointFactorization(F) -transpose(F::LU) = TransposeFactorization(F) # LU prefers transpose over adjoint in the real case, override the generic fallback +adjoint(F::LU{<:Real}) = TransposeFactorization(F) transpose(F::LU{<:Real}) = TransposeFactorization(F) # the following method is meant to catch calls to lu!(A::LAPACKArray) without a pivoting stategy From f1a50959fab5293652af335bc55db6f4dfb3a4f2 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 14 Jan 2023 13:50:19 +0100 Subject: [PATCH 05/12] fix typo --- stdlib/LinearAlgebra/src/factorization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 4f7e316444a25..7f32110ce2983 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -138,7 +138,7 @@ end function (\)(F::Factorization, B::AbstractVecOrMat) require_one_based_indexing(B) - TFB = typeof(oneunit(eltype(F)) / oneunit(eltype(B))) + TFB = typeof(oneunit(eltype(F)) \ oneunit(eltype(B))) ldiv!(F, copy_similar(B, TFB)) end (\)(F::TransposeFactorization, B::AbstractVecOrMat) = conj!(adjoint(F.parent) \ conj.(B)) From 40b9c1850979023be8a1983b7633e3ac19d71a41 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 9 Feb 2023 13:13:10 +0100 Subject: [PATCH 06/12] cody style --- stdlib/LinearAlgebra/src/factorization.jl | 4 ++-- stdlib/LinearAlgebra/src/qr.jl | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index 7f32110ce2983..bd7a2cf397594 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -123,7 +123,7 @@ end (\)(F::TransposeFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = conj!(adjoint(parent(F)) \ conj.(B)) (\)(F::AdjointFactorization{T}, B::VecOrMat{Complex{T}}) where {T<:BlasReal} = - @invoke \(F::typeof(F), B::AbstractVecOrMat) # send to meta-function in LinearAlgebra.jl + @invoke \(F::typeof(F), B::VecOrMat) function (/)(B::VecOrMat{Complex{T}}, F::Factorization{T}) where {T<:BlasReal} require_one_based_indexing(B) @@ -134,7 +134,7 @@ end (/)(B::VecOrMat{Complex{T}}, F::TransposeFactorization{T}) where {T<:BlasReal} = conj!(adjoint(parent(F)) \ conj.(B)) (/)(B::VecOrMat{Complex{T}}, F::AdjointFactorization{T}) where {T<:BlasReal} = - invoke(/, Tuple{VecOrMat{Complex{T}}, Factorization{T}}, B, F) + @invoke /(B::VecOrMat{Complex{T}}, F::Factorization{T}) function (\)(F::Factorization, B::AbstractVecOrMat) require_one_based_indexing(B) diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 107a5c6d4259d..0da25232af1ba 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -514,7 +514,8 @@ end Base.propertynames(F::QRPivoted, private::Bool=false) = (:R, :Q, :p, :P, (private ? fieldnames(typeof(F)) : ())...) -transpose(F::Union{QR,QRPivoted,QRCompactWY}) = throw(ArgumentError("transpose of QR decomposition is not supported")) +transpose(::Union{QR,QRPivoted,QRCompactWY}) = + throw(ArgumentError("transpose of QR decomposition is not supported, consider using adjoint")) abstract type AbstractQ{T} <: AbstractMatrix{T} end From 10a1201e00d237b83e965b7f2d343a714b7302c5 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 10 Feb 2023 18:39:37 +0100 Subject: [PATCH 07/12] transpose of real qr and lq --- stdlib/LinearAlgebra/src/lq.jl | 4 ++++ stdlib/LinearAlgebra/src/qr.jl | 1 + 2 files changed, 5 insertions(+) diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index fe72aad9a1b10..a162333be339a 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -135,6 +135,10 @@ AbstractArray(A::LQ) = AbstractMatrix(A) Matrix(A::LQ) = Array(AbstractArray(A)) Array(A::LQ) = Matrix(A) +transpose(F::LQ{<:Real}) = F' +transpose(::LQ) = + throw(ArgumentError("transpose of LQ decomposition is not supported, consider using adjoint")) + Base.copy(F::AdjointFactorization{T,<:LQ{T}}) where {T} = QR{T,typeof(F.parent.factors),typeof(F.parent.τ)}(copy(adjoint(F.parent.factors)), copy(F.parent.τ)) diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 0da25232af1ba..30deb785e6019 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -514,6 +514,7 @@ end Base.propertynames(F::QRPivoted, private::Bool=false) = (:R, :Q, :p, :P, (private ? fieldnames(typeof(F)) : ())...) +transpose(F::Union{QR{<:Real},QRPivoted{<:Real},QRCompactWY{<:Real}}) = F' transpose(::Union{QR,QRPivoted,QRCompactWY}) = throw(ArgumentError("transpose of QR decomposition is not supported, consider using adjoint")) From d83ebcb34c27cb6f2686f06b0ac7ecb686c689c6 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 13 Feb 2023 12:24:03 +0100 Subject: [PATCH 08/12] add NEWS and docs --- NEWS.md | 5 +++++ stdlib/LinearAlgebra/docs/src/index.md | 9 +++++++-- stdlib/LinearAlgebra/src/factorization.jl | 12 ++++++++++++ 3 files changed, 24 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 1e9f86ebb14a4..141f697b2cb55 100644 --- a/NEWS.md +++ b/NEWS.md @@ -45,8 +45,13 @@ Standard library changes packages being loaded in the Julia session. This has similar applications as the Requires.jl package but also supports precompilation and setting compatibility. + #### LinearAlgebra +* Adjoints and transposes of `Factorization` objects are no longer wrapped in `Adjoint` + and `Transpose` wrappers, respectively. Instead, they are wrapped in + `AdjointFactorization` and `TranposeFactorization` types, which themselves subtype + `Factorization` ([#46874]). #### Printf * Format specifiers now support dynamic width and precision, e.g. `%*s` and `%*.*g` ([#40105]). diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 9f12af174a4ff..62a7601fa88e0 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -276,12 +276,11 @@ to first compute the Hessenberg factorization `F` of `A` via the [`hessenberg`]( Given `F`, Julia employs an efficient algorithm for `(F+μ*I) \ b` (equivalent to `(A+μ*I)x \ b`) and related operations like determinants. - ## [Matrix factorizations](@id man-linalg-factorizations) [Matrix factorizations (a.k.a. matrix decompositions)](https://en.wikipedia.org/wiki/Matrix_decomposition) compute the factorization of a matrix into a product of matrices, and are one of the central concepts -in linear algebra. +in (numerical) linear algebra. The following table summarizes the types of matrix factorizations that have been implemented in Julia. Details of their associated methods can be found in the [Standard functions](@ref) section @@ -306,6 +305,10 @@ of the Linear Algebra documentation. | `Schur` | [Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition) | | `GeneralizedSchur` | [Generalized Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition#Generalized_Schur_decomposition) | +Adjoints and transposes of [`Factorization`](@ref) objects are lazily wrapped in +`AdjointFactorization` and `TransposeFactorization` objects, respectively. Generically, +transpose of real `Factorization`s are wrapped as `AdjointFactorization`. + ## Standard functions Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](http://www.netlib.org/lapack/). @@ -460,9 +463,11 @@ LinearAlgebra.ishermitian Base.transpose LinearAlgebra.transpose! LinearAlgebra.Transpose +LinearAlgebra.TransposeFactorization Base.adjoint LinearAlgebra.adjoint! LinearAlgebra.Adjoint +LinearAlgebra.AdjointFactorization Base.copy(::Union{Transpose,Adjoint}) LinearAlgebra.stride1 LinearAlgebra.checksquare diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index bd7a2cf397594..df4deb53164e0 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -29,7 +29,19 @@ size(F::TransposeFactorization) = reverse(size(parent(F))) size(F::Union{AdjointFactorization,TransposeFactorization}, d::Integer) = d in (1, 2) ? size(F)[d] : 1 parent(F::Union{AdjointFactorization,TransposeFactorization}) = F.parent +``` + adjoint(F::Factorization) -> AdjointFactorization + +Lazy adjoint (conjugate transposition) of the factorization `F`. +``` adjoint(F::Factorization) = AdjointFactorization(F) +``` + transpose(F::Factorization) -> TransposeFactorization + transpose(F::Factorization{<:Real}) -> AdjointFactorization + +Lazy transpose of the factorization `F`. Generically, for `Factorization`s with real +`eltype`, return an [`AdjointFactorization`](@ref). +``` transpose(F::Factorization) = TransposeFactorization(F) transpose(F::Factorization{<:Real}) = AdjointFactorization(F) adjoint(F::AdjointFactorization) = F.parent From 5927a9d725ff6e37cbee22d397c2044dd5684210 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 14 Feb 2023 11:24:35 +0100 Subject: [PATCH 09/12] fix docstrings --- stdlib/LinearAlgebra/src/factorization.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index df4deb53164e0..ae66e68bf8c5e 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -30,14 +30,14 @@ size(F::Union{AdjointFactorization,TransposeFactorization}, d::Integer) = d in ( parent(F::Union{AdjointFactorization,TransposeFactorization}) = F.parent ``` - adjoint(F::Factorization) -> AdjointFactorization + adjoint(F::Factorization)::AdjointFactorization Lazy adjoint (conjugate transposition) of the factorization `F`. ``` adjoint(F::Factorization) = AdjointFactorization(F) ``` - transpose(F::Factorization) -> TransposeFactorization - transpose(F::Factorization{<:Real}) -> AdjointFactorization + transpose(F::Factorization)::TransposeFactorization + transpose(F::Factorization{<:Real})::AdjointFactorization Lazy transpose of the factorization `F`. Generically, for `Factorization`s with real `eltype`, return an [`AdjointFactorization`](@ref). From b31311ae4ceb44738fe60137d1efcc42ce266b66 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 Feb 2023 13:05:00 +0100 Subject: [PATCH 10/12] learn how to write docstrings --- stdlib/LinearAlgebra/src/factorization.jl | 34 ++++++++++++++++------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index ae66e68bf8c5e..8c35a23e6b6d5 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -11,12 +11,26 @@ matrix factorizations. """ abstract type Factorization{T} end +""" + AdjointFactorization + +Lazy wrapper type for the adjoint of the underlying `Factorization` object. Usually, the +`AdjointFactorization` constructor should not be called directly, use +[`adjoint(:: Factorization)`](@ref) instead. +""" struct AdjointFactorization{T,S<:Factorization} <: Factorization{T} parent::S end AdjointFactorization(F::Factorization) = AdjointFactorization{Base.promote_op(adjoint,eltype(F)),typeof(F)}(F) +""" + TransposeFactorization + +Lazy wrapper type for the transpose of the underlying `Factorization` object. Usually, the +`TransposeFactorization` constructor should not be called directly, use +[`transpose(:: Factorization)`](@ref) instead. +""" struct TransposeFactorization{T,S<:Factorization} <: Factorization{T} parent::S end @@ -29,19 +43,19 @@ size(F::TransposeFactorization) = reverse(size(parent(F))) size(F::Union{AdjointFactorization,TransposeFactorization}, d::Integer) = d in (1, 2) ? size(F)[d] : 1 parent(F::Union{AdjointFactorization,TransposeFactorization}) = F.parent -``` - adjoint(F::Factorization)::AdjointFactorization +""" + adjoint(F::Factorization) -Lazy adjoint (conjugate transposition) of the factorization `F`. -``` +Lazy adjoint of the factorization `F`. By default, returns an +[`AdjointFactorization`](@ref) wrapper. +""" adjoint(F::Factorization) = AdjointFactorization(F) -``` - transpose(F::Factorization)::TransposeFactorization - transpose(F::Factorization{<:Real})::AdjointFactorization +""" + transpose(F::Factorization) -Lazy transpose of the factorization `F`. Generically, for `Factorization`s with real -`eltype`, return an [`AdjointFactorization`](@ref). -``` +Lazy transpose of the factorization `F`. By default, returns a [`TransposeFactorization`](@ref), +except for `Factorization`s with real `eltype`, in which case returns an [`AdjointFactorization`](@ref). +""" transpose(F::Factorization) = TransposeFactorization(F) transpose(F::Factorization{<:Real}) = AdjointFactorization(F) adjoint(F::AdjointFactorization) = F.parent From ba3eabf555c2915e2ab8f839cc24e1f266c0f4c0 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 Feb 2023 14:00:53 +0100 Subject: [PATCH 11/12] fix transpose qr tests --- stdlib/LinearAlgebra/test/qr.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/test/qr.jl b/stdlib/LinearAlgebra/test/qr.jl index b78e56ac0ff1b..9ccc42dfb3259 100644 --- a/stdlib/LinearAlgebra/test/qr.jl +++ b/stdlib/LinearAlgebra/test/qr.jl @@ -211,9 +211,9 @@ rectangularQ(Q::LinearAlgebra.AbstractQ) = convert(Array, Q) end @testset "transpose errors" begin - @test_throws ArgumentError transpose(qr(randn(3,3))) - @test_throws ArgumentError transpose(qr(randn(3,3), NoPivot())) - @test_throws ArgumentError transpose(qr(big.(randn(3,3)))) + @test_throws ArgumentError transpose(qr(randn(ComplexF64,3,3))) + @test_throws ArgumentError transpose(qr(randn(ComplexF64,3,3), NoPivot())) + @test_throws ArgumentError transpose(qr(big.(randn(ComplexF64,3,3)))) end @testset "Issue 7304" begin From 4b25a15ece00800dc39624b47237c53442cb85e7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 15 Feb 2023 16:31:15 +0100 Subject: [PATCH 12/12] allow transpose of real Hessenberg --- stdlib/LinearAlgebra/src/hessenberg.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/stdlib/LinearAlgebra/src/hessenberg.jl b/stdlib/LinearAlgebra/src/hessenberg.jl index 72d475d55fa44..44cd0873c1ee5 100644 --- a/stdlib/LinearAlgebra/src/hessenberg.jl +++ b/stdlib/LinearAlgebra/src/hessenberg.jl @@ -425,6 +425,10 @@ copy(F::Hessenberg{<:Any,<:SymTridiagonal}) = Hessenberg(copy(F.factors), copy(F size(F::Hessenberg, d::Integer) = size(F.H, d) size(F::Hessenberg) = size(F.H) +transpose(F::Hessenberg{<:Real}) = F' +transpose(::Hessenberg) = + throw(ArgumentError("transpose of Hessenberg decomposition is not supported, consider using adjoint")) + # iteration for destructuring into components Base.iterate(S::Hessenberg) = (S.Q, Val(:H)) Base.iterate(S::Hessenberg, ::Val{:H}) = (S.H, Val(:μ))