LU decomposition without pivoting in JULIA

Question

Trying to rewrite the lu_nopivot from this answer https://stackoverflow.com/a/41151228 into JULIA and use only one loop. The julia code I wrote

using LinearAlgebra
function lu_nopivot(A)
    n = size(A, 1)
    L = Matrix{eltype(A)}(I, n, n)
    U = copy(A)
    for k = 1:n
        L[k+1:n,k] = U[k+1:n,k] / U[k,k]
        U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k]*U[k,:]  
    end
    return L, U
end

But calling the function L, U = lu_nopivot(A) gives an error MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64}) on L[k+1:n,k]*U[k,:]. Any reason why this is?

Answer 1

The problem is that when you do U[k, :], even though you're extracting out a row, what you get back is a column vector. So L[k+1:n,k]*U[k,:] becomes an attempt to multiply a column vector by a column vector.

One way to get a row vector aka a 1 x N Matrix in Julia (though I don't know if this is the idiomatic way) is to do U[k:k, :] instead:

U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k] * U[k:k,:]  

Note however that Julia's implementation of lu already has a no-pivot option:

Pivoting can be turned off by passing pivot = NoPivot().

julia> M = rand(1.0:100.0, 3, 3)
3×3 Matrix{Float64}:
 71.0  50.0  23.0
 82.0  63.0  37.0
 96.0  28.0   5.0

julia> L, U = lu_nopivot(M); # after the above k:k change has been made

julia> L
3×3 Matrix{Float64}:
 1.0       0.0      0.0
 1.15493   1.0      0.0
 1.35211  -7.53887  1.0

julia> U
3×3 Matrix{Float64}:
 71.0  50.0      23.0
  0.0   5.25352  10.4366
  0.0   0.0      52.5818

julia> lu(M, NoPivot())
LU{Float64, Matrix{Float64}}
L factor:
3×3 Matrix{Float64}:
 1.0       0.0      0.0
 1.15493   1.0      0.0
 1.35211  -7.53887  1.0
U factor:
3×3 Matrix{Float64}:
 71.0  50.0      23.0
  0.0   5.25352  10.4366
  0.0   0.0      52.5818

Using this is likely more performant, and more robust too (eg. your current implementation cannot handle matrices of eltype Int, since L and U are given the same type as the input but will need to contain Floats).