Skip to content

Commit

Permalink
Avoid unwrapping-wrapping combination in LinearAlgebra (#41252)
Browse files Browse the repository at this point in the history
  • Loading branch information
dkarrasch authored Jun 19, 2021
1 parent 0e5fc4c commit ae1b469
Show file tree
Hide file tree
Showing 10 changed files with 142 additions and 235 deletions.
45 changes: 19 additions & 26 deletions stdlib/LinearAlgebra/src/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -731,47 +731,40 @@ ldiv!(A::Adjoint{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(copy(A), b)
function ldiv!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
require_one_based_indexing(A, B)
nA,mA = size(A)
tmp = similar(B,size(B,1))
n = size(B, 1)
if nA != n
throw(DimensionMismatch("size of A is ($nA,$mA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
copyto!(tmp, 1, B, (i - 1)*n + 1, n)
ldiv!(A, tmp)
copyto!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
for b in eachcol(B)
ldiv!(A, b)
end
B
end
function ldiv!(adjA::Adjoint{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix)
function ldiv!(adjA::Adjoint{<:Any,<:Bidiagonal}, B::AbstractMatrix)
require_one_based_indexing(adjA, B)
A = adjA.parent
nA,mA = size(A)
tmp = similar(B,size(B,1))
mA, nA = size(adjA)
n = size(B, 1)
tmp = similar(B, n)
Ac = copy(adjA)
if mA != n
throw(DimensionMismatch("size of adjoint of A is ($mA,$nA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
copyto!(tmp, 1, B, (i - 1)*n + 1, n)
ldiv!(adjoint(A), tmp)
copyto!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
for b in eachcol(B)
ldiv!(Ac, b)
end
B
end
function ldiv!(transA::Transpose{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix)
require_one_based_indexing(transA, B)
A = transA.parent
nA,mA = size(A)
tmp = similar(B,size(B,1))
function ldiv!(tA::Transpose{<:Any,<:Bidiagonal}, B::AbstractMatrix)
require_one_based_indexing(tA, B)
mA, nA = size(tA)
n = size(B, 1)
tmp = similar(B, n)
At = copy(tA)
if mA != n
throw(DimensionMismatch("size of transpose of A is ($mA,$nA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
copyto!(tmp, 1, B, (i - 1)*n + 1, n)
ldiv!(transpose(A), tmp)
copyto!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
for b in eachcol(B)
ldiv!(At, b)
end
B
end
Expand Down Expand Up @@ -824,20 +817,20 @@ function \(A::Bidiagonal{<:Number}, B::AbstractVecOrMat{<:Number})
ldiv!(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
end
\(A::Bidiagonal, B::AbstractVecOrMat) = ldiv!(A, copy(B))
function \(transA::Transpose{<:Number,<:Bidiagonal{<:Number}}, B::AbstractVecOrMat{<:Number})
A = transA.parent
function \(tA::Transpose{<:Number,<:Bidiagonal{<:Number}}, B::AbstractVecOrMat{<:Number})
A = tA.parent
TA, TB = eltype(A), eltype(B)
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
ldiv!(transpose(convert(AbstractArray{TAB}, A)), copy_oftype(B, TAB))
end
\(transA::Transpose{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(transpose(transA.parent), copy(B))
\(tA::Transpose{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(tA, copy(B))
function \(adjA::Adjoint{<:Number,<:Bidiagonal{<:Number}}, B::AbstractVecOrMat{<:Number})
A = adjA.parent
TA, TB = eltype(A), eltype(B)
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
ldiv!(adjoint(convert(AbstractArray{TAB}, A)), copy_oftype(B, TAB))
end
\(adjA::Adjoint{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(adjoint(adjA.parent), copy(B))
\(adjA::Adjoint{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(adjA, copy(B))

factorize(A::Bidiagonal) = A

Expand Down
3 changes: 2 additions & 1 deletion stdlib/LinearAlgebra/src/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -289,7 +289,8 @@ function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
rmul!(At, D)
end

*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ)
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ)

function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
A = adjA.parent
Expand Down
10 changes: 5 additions & 5 deletions stdlib/LinearAlgebra/src/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -207,9 +207,9 @@ end

### QcB
lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasReal} =
(A = adjA.parent; LAPACK.ormlq!('L','T',A.factors,A.τ,B))
(A = adjA.parent; LAPACK.ormlq!('L', 'T', A.factors, A.τ, B))
lmul!(adjA::Adjoint{<:Any,<:LQPackedQ{T}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} =
(A = adjA.parent; LAPACK.ormlq!('L','C',A.factors,A.τ,B))
(A = adjA.parent; LAPACK.ormlq!('L', 'C', A.factors, A.τ, B))

function *(adjA::Adjoint{<:Any,<:LQPackedQ}, B::StridedVecOrMat)
A = adjA.parent
Expand All @@ -232,11 +232,11 @@ function *(A::LQPackedQ, adjB::Adjoint{<:Any,<:StridedVecOrMat})
return lmul!(A, BB)
end
function *(adjA::Adjoint{<:Any,<:LQPackedQ}, adjB::Adjoint{<:Any,<:StridedVecOrMat})
A, B = adjA.parent, adjB.parent
TAB = promote_type(eltype(A), eltype(B))
B = adjB.parent
TAB = promote_type(eltype(adjA.parent), eltype(B))
BB = similar(B, TAB, (size(B, 2), size(B, 1)))
adjoint!(BB, B)
return lmul!(adjoint(A), BB)
return lmul!(adjA, BB)
end

# in-place right-application of LQPackedQs
Expand Down
Loading

2 comments on commit ae1b469

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

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

Executing the daily package evaluation, I will reply here when finished:

@nanosoldier runtests(ALL, isdaily = true)

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

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

Your package evaluation job has completed - possible new issues were detected. A full report can be found here. cc @maleadt

Please sign in to comment.