-
-
Notifications
You must be signed in to change notification settings - Fork 5.4k
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.
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
Fix LU factorization in-place operations #22774
Changes from 1 commit
8775d66
d08467c
5f2350a
944f02a
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
- Loading branch information
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -243,32 +243,50 @@ end | |
|
||
A_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = | ||
@assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info | ||
A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) = | ||
A_ldiv_B!(UpperTriangular(A.factors), | ||
A_ldiv_B!(UnitLowerTriangular(A.factors), b[ipiv2perm(A.ipiv, length(b))])) | ||
A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) = | ||
A_ldiv_B!(UpperTriangular(A.factors), | ||
A_ldiv_B!(UnitLowerTriangular(A.factors), B[ipiv2perm(A.ipiv, size(B, 1)),:])) | ||
|
||
function A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) | ||
b_permuted = b[ipiv2perm(A.ipiv, length(b))] | ||
A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), b_permuted)) | ||
copy!(b, b_permuted) | ||
end | ||
|
||
function A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) | ||
B_permuted = B[ipiv2perm(A.ipiv, size(B, 1)), :] | ||
A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), B_permuted)) | ||
copy!(B, B_permuted) | ||
end | ||
|
||
At_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = | ||
@assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, B) A.info | ||
At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) = | ||
At_ldiv_B!(UnitLowerTriangular(A.factors), | ||
At_ldiv_B!(UpperTriangular(A.factors), b))[invperm(ipiv2perm(A.ipiv, length(b)))] | ||
At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) = | ||
At_ldiv_B!(UnitLowerTriangular(A.factors), | ||
At_ldiv_B!(UpperTriangular(A.factors), B))[invperm(ipiv2perm(A.ipiv, size(B,1))),:] | ||
|
||
function At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) | ||
At_ldiv_B!(UnitLowerTriangular(A.factors), At_ldiv_B!(UpperTriangular(A.factors), b)) | ||
b_permuted = b[invperm(ipiv2perm(A.ipiv, length(b)))] | ||
copy!(b, b_permuted) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. couldn't this permute in place? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I was thinking the same, but |
||
end | ||
|
||
function At_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) | ||
At_ldiv_B!(UnitLowerTriangular(A.factors), At_ldiv_B!(UpperTriangular(A.factors), B)) | ||
B_permuted = B[invperm(ipiv2perm(A.ipiv, size(B,1))), :] | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Perhaps space-separate the arguments to |
||
copy!(B, B_permuted) | ||
end | ||
|
||
Ac_ldiv_B!(F::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:Real} = | ||
At_ldiv_B!(F, B) | ||
Ac_ldiv_B!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = | ||
@assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, B) A.info | ||
Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) = | ||
Ac_ldiv_B!(UnitLowerTriangular(A.factors), | ||
Ac_ldiv_B!(UpperTriangular(A.factors), b))[invperm(ipiv2perm(A.ipiv, length(b)))] | ||
Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) = | ||
Ac_ldiv_B!(UnitLowerTriangular(A.factors), | ||
Ac_ldiv_B!(UpperTriangular(A.factors), B))[invperm(ipiv2perm(A.ipiv, size(B,1))),:] | ||
|
||
function Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, b::StridedVector) | ||
Ac_ldiv_B!(UnitLowerTriangular(A.factors), Ac_ldiv_B!(UpperTriangular(A.factors), b)) | ||
b_permuted = b[invperm(ipiv2perm(A.ipiv, length(b)))] | ||
copy!(b, b_permuted) | ||
end | ||
|
||
function Ac_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedMatrix) | ||
Ac_ldiv_B!(UnitLowerTriangular(A.factors), Ac_ldiv_B!(UpperTriangular(A.factors), B)) | ||
B_permuted = B[invperm(ipiv2perm(A.ipiv, size(B,1))), :] | ||
copy!(B, B_permuted) | ||
end | ||
|
||
At_ldiv_Bt(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = | ||
@assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it is sufficient to do
A_ldiv_B!(UpperTriangular(A[:U]), A_ldiv_B!(LinAlg.UnitLowerTriangular(A[:L]), view(b, A[:p])))
here.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have tried exactly this, but the problem is that the user works with
b
itself after callingA_ldiv_B!(LU, b)
, andb
has to be permuted. So whether you do the permutation during the backward & forward substitutions implicitly with a view doesn't really matter, the permutation has to be applied in the end anyway.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I didn't know about
A[:p]
btw, thanks