Fix LU factorization in-place operations

[haampie]

Fixes #22683.

[tkelman]

couldn't this permute in place?

[haampie]

I was thinking the same, but permute! copies the permutation array, so a new allocation of roughly the same size is done anyway. Also the docs suggest using b[perm]. But if the array or number type is special, then permute! might be better.

[Sacha0]

Regarding three-argument [i]permute! methods, having run into the same lack last week I opened an issue: #22815. Best! (Wonderful meeting you last month Harmen! :) )

[Sacha0]

Perhaps space-separate the arguments to size as in the similar method above? :)

[Sacha0]

Were these eltype restrictions unnecessary? :)

[haampie]

' and .' are equivalent for real matrices but not for complex ones, and I think the idea has been to restrict the .' test to complex matrices. But then it should have read <: Complex. When fixing that I thought it wasn't really necessary to do that test conditionally; makes it easier to read what's happening.

[Sacha0]

Perhaps DRY out these tests? For example, something along the lines of

for (Aq_ldiv_B!, q) in (
        (A_ldiv_B!, identity),
        (At_ldiv_B!, transpose),
        (Ac_ldiv_B!, ctranspose) )
    Aq_ldiv_B!(b_dest, lua, b)
    Aq_ldiv_B!(c_dest, lua, c)
    @test norm(b_dest - q(lua) \ b, 1) < ε*κ*2n
    @test norm(c_dest - q(lua) \ c, 1) < ε*κ*n
end

might serve? (Edit: Missed the secondary op the first time around! :) )

[haampie]

I've tried it, but it doesn't seem to work :(

> transpose(lufact(rand(3, 3))) \ rand(3)
ERROR: transpose not implemented for Base.LinAlg.LU{Float64,Array{Float64,2}}

[fredrikekre]

It is because A.' \ B becomes At_ldiv_B(A, B) (not transpose(A) \ B). Perhaps this will do?

for (Aq_ldiv_B!, f) in (
        (A_ldiv_B!, \),
        (At_ldiv_B!, At_ldiv_B),
        (Ac_ldiv_B!, Ac_ldiv_B) )
    Aq_ldiv_B!(b_dest, lua, b)
    Aq_ldiv_B!(c_dest, lua, c)
    @test norm(b_dest - f(lua, b), 1) < ε*κ*2n
    @test norm(c_dest - f(lua, c), 1) < ε*κ*n
end

[haampie]

Thx, so is .' \ parsed as a single operator or something?

[fredrikekre]

Yes,

julia> expand(Main, :(A.'\B))
:(At_ldiv_B(A, B))

and similarly for * as well.

[StefanKarpinski]

This is definitely a hack and we want to get rid of it, but we're not quite there yet.

[Sacha0]

Not necessary in this pull request, but DRYing out these and the similar method definitions below might be nice. For example, as a first step perhaps something along the lines of

function A_ldiv_B!(A::LU{<:Any,<:StridedMatrix}, B::StridedVecOrMat)
    B_permuted = _ipivpermute(B, A.ipiv)
    A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), B_permuted))
    copy!(B, B_permuted)
end
_ipivpermute(b::StridedVector, ipiv) = b[ipiv2perm(ipiv, length(b))]
_ipivpermute(B::StridedMatrix, ipiv) = B[ipiv2perm(ipiv, size(B, 1)), :]

would do the trick? :)

[haampie]

Thanks @Sacha0! Nice meeting you as well. I'll condense it a bit :)

[andreasnoack]

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.

[haampie]

I have tried exactly this, but the problem is that the user works with b itself after calling A_ldiv_B!(LU, b), and b 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.

[haampie]

I didn't know about A[:p] btw, thanks

[haampie]

OK... We should have used the ipiv array (I didn't even consider it), since that provides a way to sequentially apply the permutation in-place:

function my_A_ldiv_B!(A::LU, b::StridedVector)
    for i = eachindex(A.ipiv)
        b[i], b[A.ipiv[i]] = b[A.ipiv[i]], b[i]
    end

    A_ldiv_B!(UpperTriangular(A.factors), A_ldiv_B!(UnitLowerTriangular(A.factors), b))
end