Fix generic triangular solves with empty matrices (JuliaLang#54201)

Author: dkarrasch

stdlib/LinearAlgebra/src/triangular.jl

@@ -1296,6 +1296,7 @@ function generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function,
     if size(C) != size(B)
         throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
     end
+    iszero(mA) && return C
     oA = oneunit(eltype(A))
     @inbounds if uploc == 'U'
         if isunitc == 'N'
@@ -1432,6 +1433,7 @@ function generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, xA:
     if size(C) != size(B)
         throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))"))
     end
+    iszero(mA) && return C
     oA = oneunit(eltype(A))
     @inbounds if uploc == 'U'
         if isunitc == 'N'

stdlib/LinearAlgebra/test/lu.jl

@@ -486,4 +486,17 @@ end
     LinearAlgebra.generic_lufact!(fill(Inf, 2, 2), check=false)
 end
 
+@testset "lu for empty matrices" begin
+    for T in (Float64, BigFloat)
+        A = fill(T(0.0), 0, 0)
+        v = fill(T(1.0), 0, 10)
+        @test A \ v ≈ lu(A) \ v
+        vt = permutedims(v)
+        @test vt / A ≈ vt / lu(A)
+        B = UpperTriangular(transpose(fill(complex(T(0.0)), 0, 0)'))
+        @test B \ v ≈ v
+        @test vt / B ≈ vt
+    end
+end
+
 end # module TestLU