Diagonal for lower triangular of LU decomposition set incorrectly (#1813)

Author: tgymnich

lib/cusolver/linalg.jl

@@ -455,7 +455,7 @@ function Base.getproperty(F::LU{T,<:StridedCuMatrix}, d::Symbol) where T
     m, n = size(F)
     if d === :L
         L = tril!(getfield(F, :factors)[1:m, 1:min(m,n)])
-        L[1:min(m,n)+1:end] .= one(T)   # set the diagonal (linear indexing trick)
+        L[1:m+1:end] .= one(T)   # set the diagonal (linear indexing trick)
         return L
     elseif VERSION >= v"1.9.0-DEV.1775"
         invoke(getproperty, Tuple{LU{T}, Symbol}, F, d)

test/cusolver/dense.jl

@@ -602,6 +602,18 @@ k = 1
             lu_gpu = lu(A_d)
             @test ldiv!(lu_cpu, B) ≈ collect(ldiv!(lu_gpu, B_d))
         end
+        
+        A = CuMatrix(rand(1024, 1024))
+        lua = lu(A)
+        @test Matrix(lua.L) * Matrix(lua.U) ≈ Matrix(lua.P) * Matrix(A)
+
+        A = CuMatrix(rand(1024, 512))
+        lua = lu(A)
+        @test Matrix(lua.L) * Matrix(lua.U) ≈ Matrix(lua.P) * Matrix(A)
+    
+        A = CuMatrix(rand(512, 1024))
+        lua = lu(A)
+        @test Matrix(lua.L) * Matrix(lua.U) ≈ Matrix(lua.P) * Matrix(A)
     end
 end