missing/bad methods for ldiv! with dense triangular?

[chriscoey]

I am on Julia nightly. The docs say ldiv!(Y, A, B) -> Y computes A \ B in-place and stores the result in Y. I need to use this 3-arg in-place function to solve L*x = y for x, where L is a LowerTriangular matrix. I run:

L = LowerTriangular(rand(5,5))
y = rand(5)
x = L\y    # works correctly, but not in-place
ldiv!(x, L, y)    # MethodError: no method matching ldiv!(::Array{Float64,1}, 
    ::LowerTriangular{Float64,Array{Float64,2}}, ::Array{Float64,1})

The method error says that an available method is:

ldiv!(::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::Factorization, 
    ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) 

so I tried ldiv!(x, factorize(L) , y). But this gives the same method error, because factorize(L) returns L and LowerTriangular is not a subtype of Factorization
(Aside: it feels a little counterintuitive to me that for most matrices A, factorize(A) returns a Factorization, but for some (special) matrices A (like a LowerTriangular or a Diagonal) it doesn't.)

I saw the definition:

julia/stdlib/LinearAlgebra/src/triangular.jl

Line 1274 in d55b044

function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector)

but this corresponds to ldiv!(A, B, Y) -> Y, which is inconsistent with the standard documented ldiv!(Y, A, B) -> Y (the latter is consistent with mul!(Y, A, B) -> Y). I run:

ldiv!(copy(L')', y, x)    # stores L\y in x
ldiv!(L, y, x)    # method error

which convinced me this is an issue. That inconsistent ldiv! method shouldn't be allowed anywhere, right?
(Aside: the names of lmul! and rmul! (which take 2 arguments and destroy one) and ldiv! and rdiv! (which take 3 arguments and update one) are a little confusing to me because they imply an analogy that doesn't exist.)

[chriscoey]

I think the following methods (modified from @Sacha0's code in JuliaLang/LinearAlgebra.jl#366) correctly implement ldiv!(Y,A,B) -> Y for A::LowerTriangular. ie they do not destroy B, only change Y, consistent with the generic docstring.

function ldiv!(Y::AbstractVector, A::LowerTriangular, B::AbstractVector)
    m = length(Y)
    @assert m == size(A, 1) == length(B)
    Y .= B
    for i in 1:m
        @inbounds iszero(A.data[i,i]) && throw(SingularException(i))
        @inbounds Yi = Y[i] /= A.data[i,i]
        for k in i+1:m
            @inbounds Y[k] -= A.data[k,i] * Yi
        end
    end
    return Y
end

function ldiv!(Y::AbstractMatrix, A::LowerTriangular, B::AbstractMatrix)
    (m, n) = size(Y)
    @assert m == size(A, 1) == size(B, 1)
    @assert n == size(B, 2)
    Y .= B
    for i in 1:m
        @inbounds iszero(A.data[i,i]) && throw(SingularException(i))
        Aii = A.data[i,i]
        for j in 1:n
            @inbounds Yij = Y[i,j] /= Aii
            for k in i+1:m
                @inbounds Y[k,j] -= A.data[k,i] * Yij
            end
        end
    end
    return Y
end

[andreasnoack]

This is no longer an issue

julia> ldiv!(x, L, y)
5-element Array{Float64,1}:
  0.6643442988466239
  0.9702264861304267
  0.4002865909862719
 -0.08221672211173842
 -1.8072733080260066