From 7275637466791b2554502c0860d9d8006ac90ad7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Sun, 1 Dec 2024 18:34:25 -0300
Subject: [PATCH 1/9] define lu!(F::LU,A) for general matrices

---
 src/lu.jl | 36 ++++++++++++++++++++++++++++++++++--
 1 file changed, 34 insertions(+), 2 deletions(-)

diff --git a/src/lu.jl b/src/lu.jl
index 0837ac08..16f980c2 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -102,6 +102,39 @@ function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupi
     end
     lu!(A.data, pivot; check, allowsingular)
 end
+
+#reusing LU object
+#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
+#but special care needs to be done in the HermOrSym case
+function _lu_copy!(A,x)
+    copyto!(A, x)
+end
+
+function lu_copy!(A,x::HermOrSym)
+    copytri!(A.data, x.uplo, isa(x, Hermitian))
+    @inbounds if isa(x, Hermitian) # realify diagonal
+        for i in axes(x, 1)
+            A[i,i] = x[i,i]
+        end
+    end
+    return A
+end
+
+function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
+    lu_copy!(F.factors,A)
+    lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
+    check && _check_lu_success(lpt[3], allowsingular)
+    return F
+end
+
+function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
+    lu_copy!(F.factors,A)
+    generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
+    return F
+end
+
+
+
 # for backward compatibility
 # TODO: remove towards Julia v2
 @deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, RowMaximum(); check=check)
@@ -149,7 +182,7 @@ Stacktrace:
 """
 lu!(A::AbstractMatrix, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(eltype(A));
     check::Bool = true, allowsingular::Bool = false) = generic_lufact!(A, pivot; check, allowsingular)
-function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T);
+function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T), ipiv::AbstractVector{BlasInt} = Vector{BlasInt}(undef,min(size(A)...));
                          check::Bool = true, allowsingular::Bool = false) where {T}
     check && LAPACK.chkfinite(A)
     # Extract values
@@ -158,7 +191,6 @@ function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,R
 
     # Initialize variables
     info = 0
-    ipiv = Vector{BlasInt}(undef, minmn)
     @inbounds begin
         for k = 1:minmn
             # find index max

From 14d2196b3746e265dbee30447b6af6d9244fa0fb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Mon, 2 Dec 2024 18:01:51 -0300
Subject: [PATCH 2/9] remove lu_copy!

---
 src/lu.jl | 18 ++----------------
 1 file changed, 2 insertions(+), 16 deletions(-)

diff --git a/src/lu.jl b/src/lu.jl
index 16f980c2..d44f5866 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -105,30 +105,16 @@ end
 
 #reusing LU object
 #lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
-#but special care needs to be done in the HermOrSym case
-function _lu_copy!(A,x)
-    copyto!(A, x)
-end
-
-function lu_copy!(A,x::HermOrSym)
-    copytri!(A.data, x.uplo, isa(x, Hermitian))
-    @inbounds if isa(x, Hermitian) # realify diagonal
-        for i in axes(x, 1)
-            A[i,i] = x[i,i]
-        end
-    end
-    return A
-end
 
 function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
-    lu_copy!(F.factors,A)
+    copyto!(F.factors,A)
     lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
     check && _check_lu_success(lpt[3], allowsingular)
     return F
 end
 
 function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
-    lu_copy!(F.factors,A)
+    copyto!(F.factors,A)
     generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
     return F
 end

From b1247380c34ef4913d068614e328eedca46f45ae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Mon, 2 Dec 2024 18:09:38 -0300
Subject: [PATCH 3/9] Update src/lu.jl

Co-authored-by: Steven G. Johnson <stevenj@mit.edu>
---
 src/lu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lu.jl b/src/lu.jl
index d44f5866..8e1a226e 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -110,7 +110,7 @@ function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, a
     copyto!(F.factors,A)
     lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
     check && _check_lu_success(lpt[3], allowsingular)
-    return F
+    return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
 end
 
 function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)

From 1bd4b06b61ccc2dabc0b9c80c2cd22e17a6a281c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Mon, 2 Dec 2024 18:10:30 -0300
Subject: [PATCH 4/9] add whitespace

---
 src/lu.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lu.jl b/src/lu.jl
index 8e1a226e..fd36281c 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -107,14 +107,14 @@ end
 #lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
 
 function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
-    copyto!(F.factors,A)
+    copyto!(F.factors, A)
     lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
     check && _check_lu_success(lpt[3], allowsingular)
     return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
 end
 
 function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
-    copyto!(F.factors,A)
+    copyto!(F.factors, A)
     generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
     return F
 end

From b8d1ecc050487931e7f79c8362f77d6a03c78d61 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Mon, 2 Dec 2024 18:11:39 -0300
Subject: [PATCH 5/9] also return a new LU object in the generic case

---
 src/lu.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/lu.jl b/src/lu.jl
index fd36281c..52da7b7f 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -115,8 +115,7 @@ end
 
 function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
     copyto!(F.factors, A)
-    generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
-    return F
+    return generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
 end
 
 

From 721f0778a8ee72705aeec53962c30973038eb228 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Thu, 12 Dec 2024 00:48:07 -0300
Subject: [PATCH 6/9] docs for lu!(F::LU,A)

---
 src/lu.jl | 57 +++++++++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 49 insertions(+), 8 deletions(-)

diff --git a/src/lu.jl b/src/lu.jl
index 52da7b7f..885f9826 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -103,22 +103,63 @@ function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupi
     lu!(A.data, pivot; check, allowsingular)
 end
 
-#reusing LU object
 #lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
+"""
+    lu!(F::LU, pivot = RowMaximum(); check = true, allowsingular = false) -> LU
 
-function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
-    copyto!(F.factors, A)
-    lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
-    check && _check_lu_success(lpt[3], allowsingular)
-    return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
-end
+`lu!` is the same as [`lu`](@ref), but saves space by overwriting the
+input `F`, instead of creating a copy.
+
+!!! compat "Julia 1.12"
+    reusing `LU` factorizations in `lu!` require Julia 1.12 or later.
+
+# Examples
+```jldoctest
+julia> A = [4. 3.; 6. 3.]
+2×2 Matrix{Float64}:
+ 4.0  3.0
+ 6.0  3.0
+
+julia> F = lu(A)
+LU{Float64, Matrix{Float64}, Vector{Int64}}
+L factor:
+2×2 Matrix{Float64}:
+ 1.0       0.0
+ 0.666667  1.0
+U factor:
+2×2 Matrix{Float64}:
+ 6.0  3.0
+ 0.0  1.0
+
+julia> B = [8 3; 12 3]
+2×2 Matrix{Int64}:
+ 8  3
+ 12  3
 
+julia> F2 = lu!(F,A)
+LU{Float64, Matrix{Float64}, Vector{Int64}}
+L factor:
+2×2 Matrix{Float64}:
+ 1.0       0.0
+ 0.666667  1.0
+U factor:
+2×2 Matrix{Float64}:
+ 12.0  3.0
+  0.0  1.0
+```
+"""
 function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
     copyto!(F.factors, A)
     return generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
 end
 
-
+#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
+function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
+    copyto!(F.factors, A)
+    lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
+    check && _check_lu_success(lpt[3], allowsingular)
+    return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
+end
 
 # for backward compatibility
 # TODO: remove towards Julia v2

From 8a5576125b88e77538f2ba5dc7fba5ee484ccb17 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Thu, 12 Dec 2024 00:52:50 -0300
Subject: [PATCH 7/9] add a square lu!(F::LU,A) test

---
 test/lu.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/lu.jl b/test/lu.jl
index 56a402d7..20358f30 100644
--- a/test/lu.jl
+++ b/test/lu.jl
@@ -67,6 +67,7 @@ dimg  = randn(n)/2
         @test (l*u)[invperm(p),:] ≈ a
         @test a * inv(lua) ≈ Matrix(I, n, n)
         @test copy(lua) == lua
+        @test lu!(copy(lua),a) == lua
         if eltya <: BlasFloat
             # test conversion of LU factorization's numerical type
             bft = eltya <: Real ? LinearAlgebra.LU{BigFloat} : LinearAlgebra.LU{Complex{BigFloat}}

From fe08274a1c392e4292973f0da9823dc7f57976b7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Fri, 13 Dec 2024 17:48:52 -0300
Subject: [PATCH 8/9] typo

---
 src/lu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lu.jl b/src/lu.jl
index 885f9826..601e6288 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -158,7 +158,7 @@ function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, a
     copyto!(F.factors, A)
     lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
     check && _check_lu_success(lpt[3], allowsingular)
-    return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
+    return LU{eltype(lpt[1]),typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
 end
 
 # for backward compatibility

From 4839167e30a153547af7134c0e8f82d12861f97a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9s=20Riedemann?=
 <38795484+longemen3000@users.noreply.github.com>
Date: Sat, 14 Dec 2024 15:42:10 -0300
Subject: [PATCH 9/9] Update src/lu.jl

Co-authored-by: Steven G. Johnson <stevenj@mit.edu>
---
 src/lu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/lu.jl b/src/lu.jl
index 601e6288..b2d8380f 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -111,7 +111,7 @@ end
 input `F`, instead of creating a copy.
 
 !!! compat "Julia 1.12"
-    reusing `LU` factorizations in `lu!` require Julia 1.12 or later.
+    Reusing dense `LU` factorizations in `lu!` require Julia 1.12 or later.
 
 # Examples
 ```jldoctest