From 0fb5d6f195a12d267dbca3c39484c6165562131c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Fauchereau?=
 <clement.fauchereau@ensta-bretagne.org>
Date: Thu, 29 Apr 2021 17:03:25 +0200
Subject: [PATCH] Optimisation of getproperty for UmfpackLU object (avoid
 copying irrelevant information)

---
 stdlib/SuiteSparse/src/umfpack.jl  | 299 ++++++++++++++++++++++-------
 stdlib/SuiteSparse/test/umfpack.jl |   5 +
 2 files changed, 231 insertions(+), 73 deletions(-)

diff --git a/stdlib/SuiteSparse/src/umfpack.jl b/stdlib/SuiteSparse/src/umfpack.jl
index 4d6593a34aad9..e1e5df876c9a5 100644
--- a/stdlib/SuiteSparse/src/umfpack.jl
+++ b/stdlib/SuiteSparse/src/umfpack.jl
@@ -428,61 +428,237 @@ for itype in UmfpackIndexTypes
                            lnz, unz, n_row, n_col, nz_diag, lu.numeric)
             (lnz[], unz[], n_row[], n_col[], nz_diag[])
         end
-        function umf_extract(lu::UmfpackLU{Float64,$itype})
-            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
-            (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
-            Lp = Vector{$itype}(undef, n_row + 1)
-            Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
-            Lx = Vector{Float64}(undef, lnz)
-            Up = Vector{$itype}(undef, n_col + 1)
-            Ui = Vector{$itype}(undef, unz)
-            Ux = Vector{Float64}(undef, unz)
-            P  = Vector{$itype}(undef, n_row)
-            Q  = Vector{$itype}(undef, n_col)
-            Rs = Vector{Float64}(undef, n_row)
-            @isok ccall(($get_num_r,:libumfpack), $itype,
-                        (Ptr{$itype},Ptr{$itype},Ptr{Float64},
-                         Ptr{$itype},Ptr{$itype},Ptr{Float64},
-                         Ptr{$itype},Ptr{$itype},Ptr{Cvoid},
-                         Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
-                        Lp,Lj,Lx,
-                        Up,Ui,Ux,
-                        P, Q, C_NULL,
-                        0, Rs, lu.numeric)
-            (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), Lx))),
-             SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), Ux),
-             increment!(P), increment!(Q), Rs)
+        function getproperty(lu::UmfpackLU{Float64, $itype}, d::Symbol)
+            if d === :L
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Lp = Vector{$itype}(undef, n_row + 1)
+                # L is returned in CSR (compressed sparse row) format
+                Lj = Vector{$itype}(undef, lnz)
+                Lx = Vector{Float64}(undef, lnz)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            Lp, Lj, Lx,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
+                                                      increment!(Lp), increment!(Lj), Lx)))
+            elseif d === :U
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Up = Vector{$itype}(undef, n_col + 1)
+                Ui = Vector{$itype}(undef, unz)
+                Ux = Vector{Float64}(undef, unz)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL,
+                            Up, Ui, Ux,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return  SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
+                                        increment!(Ui), Ux)
+            elseif d === :p
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                P  = Vector{$itype}(undef, n_row)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL,
+                            P, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return increment!(P)
+            elseif d === :q
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Q  = Vector{$itype}(undef, n_col)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{$itype}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, Q, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return increment!(Q)
+            elseif d === :Rs
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Rs = Vector{Float64}(undef, n_row)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL,
+                            C_NULL, Rs, lu.numeric)
+                return Rs
+            elseif d === :(:)
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Lp = Vector{$itype}(undef, n_row + 1)
+                # L is returned in CSR (compressed sparse row) format
+                Lj = Vector{$itype}(undef, lnz)
+                Lx = Vector{Float64}(undef, lnz)
+                Up = Vector{$itype}(undef, n_col + 1)
+                Ui = Vector{$itype}(undef, unz)
+                Ux = Vector{Float64}(undef, unz)
+                P  = Vector{$itype}(undef, n_row)
+                Q  = Vector{$itype}(undef, n_col)
+                Rs = Vector{Float64}(undef, n_row)
+                @isok ccall(($get_num_r, :libumfpack), $itype,
+                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Float64},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
+                            Lp, Lj, Lx,
+                            Up, Ui, Ux,
+                            P, Q, C_NULL,
+                            C_NULL, Rs, lu.numeric)
+                return (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
+                                                       increment!(Lp), increment!(Lj),
+                                                       Lx))),
+                        SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
+                                        increment!(Ui), Ux),
+                        increment!(P), increment!(Q), Rs)
+            else
+                return getfield(lu, d)
+            end
         end
-        function umf_extract(lu::UmfpackLU{ComplexF64,$itype})
-            umfpack_numeric!(lu)        # ensure the numeric decomposition exists
-            (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
-            Lp = Vector{$itype}(undef, n_row + 1)
-            Lj = Vector{$itype}(undef, lnz) # L is returned in CSR (compressed sparse row) format
-            Lx = Vector{Float64}(undef, lnz)
-            Lz = Vector{Float64}(undef, lnz)
-            Up = Vector{$itype}(undef, n_col + 1)
-            Ui = Vector{$itype}(undef, unz)
-            Ux = Vector{Float64}(undef, unz)
-            Uz = Vector{Float64}(undef, unz)
-            P  = Vector{$itype}(undef, n_row)
-            Q  = Vector{$itype}(undef, n_col)
-            Rs = Vector{Float64}(undef, n_row)
-            @isok ccall(($get_num_z,:libumfpack), $itype,
-                        (Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
-                         Ptr{$itype},Ptr{$itype},Ptr{Float64},Ptr{Float64},
-                         Ptr{$itype},Ptr{$itype},Ptr{Cvoid}, Ptr{Cvoid},
-                         Ref{$itype},Ptr{Float64},Ptr{Cvoid}),
-                        Lp,Lj,Lx,Lz,
-                        Up,Ui,Ux,Uz,
-                        P, Q, C_NULL, C_NULL,
-                        0, Rs, lu.numeric)
-            (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row, increment!(Lp), increment!(Lj), complex.(Lx, Lz)))),
-             SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up), increment!(Ui), complex.(Ux, Uz)),
-             increment!(P), increment!(Q), Rs)
+        function getproperty(lu::UmfpackLU{ComplexF64, $itype}, d::Symbol)
+            if d === :L
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Lp = Vector{$itype}(undef, n_row + 1)
+                # L is returned in CSR (compressed sparse row) format
+                Lj = Vector{$itype}(undef, lnz)
+                Lx = Vector{Float64}(undef, lnz)
+                Lz = Vector{Float64}(undef, lnz)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            Lp, Lj, Lx, Lz,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
+                                                      increment!(Lp), increment!(Lj),
+                                                      complex.(Lx, Lz))))
+            elseif d === :U
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Up = Vector{$itype}(undef, n_col + 1)
+                Ui = Vector{$itype}(undef, unz)
+                Ux = Vector{Float64}(undef, unz)
+                Uz = Vector{Float64}(undef, unz)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            Up, Ui, Ux, Uz,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
+                                       increment!(Ui), complex.(Ux, Uz))
+            elseif d === :p
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                P  = Vector{$itype}(undef, n_row)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            P, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return increment!(P)
+            elseif d === :q
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Q  = Vector{$itype}(undef, n_col)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, Q, C_NULL, C_NULL,
+                            C_NULL, C_NULL, lu.numeric)
+                return increment!(Q)
+            elseif d === :Rs
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Rs = Vector{Float64}(undef, n_row)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, C_NULL, C_NULL, C_NULL,
+                            C_NULL, Rs, lu.numeric)
+                return Rs
+            elseif d === :(:)
+                umfpack_numeric!(lu)        # ensure the numeric decomposition exists
+                (lnz, unz, n_row, n_col, nz_diag) = umf_lunz(lu)
+                Lp = Vector{$itype}(undef, n_row + 1)
+                # L is returned in CSR (compressed sparse row) format
+                Lj = Vector{$itype}(undef, lnz)
+                Lx = Vector{Float64}(undef, lnz)
+                Lz = Vector{Float64}(undef, lnz)
+                Up = Vector{$itype}(undef, n_col + 1)
+                Ui = Vector{$itype}(undef, unz)
+                Ux = Vector{Float64}(undef, unz)
+                Uz = Vector{Float64}(undef, unz)
+                P  = Vector{$itype}(undef, n_row)
+                Q  = Vector{$itype}(undef, n_col)
+                Rs = Vector{Float64}(undef, n_row)
+                @isok ccall(($get_num_z, :libumfpack), $itype,
+                            (Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Float64}, Ptr{Float64},
+                             Ptr{$itype}, Ptr{$itype}, Ptr{Cvoid}, Ptr{Cvoid},
+                             Ptr{Cvoid}, Ptr{Float64}, Ptr{Cvoid}),
+                            Lp, Lj, Lx, Lz,
+                            Up, Ui, Ux, Uz,
+                            P, Q, C_NULL, C_NULL,
+                            C_NULL, Rs, lu.numeric)
+                return (copy(transpose(SparseMatrixCSC(min(n_row, n_col), n_row,
+                                                       increment!(Lp), increment!(Lj),
+                                                       complex.(Lx, Lz)))),
+                        SparseMatrixCSC(min(n_row, n_col), n_col, increment!(Up),
+                                        increment!(Ui), complex.(Ux, Uz)),
+                        increment!(P), increment!(Q), Rs)
+            else
+                return getfield(lu, d)
+            end
         end
     end
 end
 
+# backward compatibility
+umfpack_extract(lu::UmfpackLU) = getproperty(lu, :(:))
+
 function nnz(lu::UmfpackLU)
     lnz, unz, = umf_lunz(lu)
     return Int(lnz + unz)
@@ -555,29 +731,6 @@ function _AqldivB_kernel!(X::StridedMatrix{Tb}, lu::UmfpackLU{Float64},
     end
 end
 
-
-@inline function getproperty(lu::UmfpackLU, d::Symbol)
-    if d === :L || d === :U || d === :p || d === :q || d === :Rs || d === :(:)
-        # Guard the call to umf_extract behaind a branch to avoid infinite recursion
-        L, U, p, q, Rs = umf_extract(lu)
-        if d === :L
-            return L
-        elseif d === :U
-            return U
-        elseif d === :p
-            return p
-        elseif d === :q
-            return q
-        elseif d === :Rs
-            return Rs
-        elseif d === :(:)
-            return (L, U, p, q, Rs)
-        end
-    else
-        getfield(lu, d)
-    end
-end
-
 for Tv in (:Float64, :ComplexF64), Ti in UmfpackIndexTypes
     f = Symbol(umf_nm("free_symbolic", Tv, Ti))
     @eval begin
diff --git a/stdlib/SuiteSparse/test/umfpack.jl b/stdlib/SuiteSparse/test/umfpack.jl
index a4f749f1ce58b..4dcf7007aa1c7 100644
--- a/stdlib/SuiteSparse/test/umfpack.jl
+++ b/stdlib/SuiteSparse/test/umfpack.jl
@@ -27,6 +27,11 @@ using SparseArrays: nnz, sparse, sprand, sprandn, SparseMatrixCSC
             @test nnz(lua) == 18
             @test_throws ErrorException lua.Z
             L,U,p,q,Rs = lua.:(:)
+            @test L == lua.L
+            @test U == lua.U
+            @test p == lua.p
+            @test q == lua.q
+            @test Rs == lua.Rs
             @test (Diagonal(Rs) * A)[p,q] ≈ L * U
 
             det(lua) ≈ det(Array(A))