From d508952853f87f8415ed8c6b3ff79f749f3977da Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Wed, 9 Jul 2025 02:13:13 +0530
Subject: [PATCH 01/22] Implement RowMaximum Pivoting Strategy for Distributed
 LU Factorization

- **Implemented RowMaximum Pivoting**: Added a new LU factorization strategy using the RowMaximum pivoting method for distributed matrices

- **Custom Pivot Search and Swapping**: Introduced helper functions for searching row maxima, updating pivot indices, and swapping rows in both panel and trailing submatrices

- **Blockwise Distributed Algorithm**: Ensured compatibility with block-partitioned distributed matrices, supporting only equal block sizes for now

- **Non-breaking Addition**: Existing NoPivot LU functionality remains unchanged; RowMaximum is an additional strategy selectable via the LinearAlgebra interface.
---
 src/Dagger.jl           |  2 +-
 src/array/lu.jl         | 96 ++++++++++++++++++++++++++++++++++++++++-
 test/array/linalg/lu.jl | 28 ++++++------
 3 files changed, 111 insertions(+), 15 deletions(-)

diff --git a/src/Dagger.jl b/src/Dagger.jl
index 0c3761c44..68fd0d89c 100644
--- a/src/Dagger.jl
+++ b/src/Dagger.jl
@@ -7,7 +7,7 @@ import SparseArrays: sprand, SparseMatrixCSC
 import MemPool
 import MemPool: DRef, FileRef, poolget, poolset
 
-import Base: collect, reduce
+import Base: collect, reduce, view
 
 import LinearAlgebra
 import LinearAlgebra: Adjoint, BLAS, Diagonal, Bidiagonal, Tridiagonal, LAPACK, LowerTriangular, PosDefException, Transpose, UpperTriangular, UnitLowerTriangular, UnitUpperTriangular, diagind, ishermitian, issymmetric
diff --git a/src/array/lu.jl b/src/array/lu.jl
index 6a0b14680..979eb4648 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -7,8 +7,6 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
     mzone = -one(T)
     Ac = A.chunks
     mt, nt = size(Ac)
-    iscomplex = T <: Complex
-    trans = iscomplex ? 'C' : 'T'
 
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
@@ -29,3 +27,97 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
 
     return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)
 end
+
+function searchmax_pivot!(piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}, A::AbstractArray{T}, offset::Int=0) where T
+    max_idx = argmax(abs.(A[:]))
+    piv_idx[1] = offset+max_idx
+    piv_val[1] = A[max_idx]
+    println("searchmax_pivot: ", piv_idx[1], "\n", abs(piv_val[1]))
+end
+
+function update_ipiv!(ipivl, piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}, k::Int, nb::Int) where T
+    max_piv_idx = argmax(abs.(piv_val))
+    ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
+    println("update_ipiv: ", ipivl[1])
+end
+
+function swaprows_panel!(A::AbstractArray{T}, M::AbstractArray{T}, ipivl::AbstractVector{Int}, m::Int, p::Int, nb::Int) where T
+    q = div(ipivl[1]-1,nb) + 1
+    r = (ipivl[1]-1)%nb+1
+    if m == q
+        A[p,:], M[r,:] = M[r,:], A[p,:]
+        println("swaprows_panel: ", imag.(A[p,:]), "\n", imag.(M[r,:]))
+    end
+end
+
+function update_panel!(M::AbstractArray{T}, A::AbstractArray{T}, p::Int) where T
+    Acinv = one(T) / A[p,p]  
+    LinearAlgebra.BLAS.scal!(Acinv, view(M, :, p))
+    LinearAlgebra.BLAS.ger!(-one(T), view(M, :, p), conj.(view(A, p, p+1:size(A,2))), view(M, :, p+1:size(M,2)))
+end
+
+function swaprows_trail!(A::AbstractArray{T}, M::AbstractArray{T}, ipiv::AbstractVector{Int}, m::Int, nb::Int) where T
+    for p in eachindex(ipiv)
+        q = div(ipiv[p]-1,nb) + 1
+        r = (ipiv[p]-1)%nb+1
+        if m == q
+            A[p,:], M[r,:] = M[r,:], A[p,:]
+        println("swaprows_trail: ", imag.(A[p,:]), "\n", imag.(M[r,:]))
+        end
+    end
+end
+
+function LinearAlgebra.lu(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool=true) where T
+    A_copy = LinearAlgebra._lucopy(A, LinearAlgebra.lutype(T))
+    return LinearAlgebra.lu!(A_copy, LinearAlgebra.RowMaximum(); check=check)
+end
+function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool=true) where T
+    zone = one(T)
+    mzone = -one(T)
+
+    Ac = A.chunks
+    mt, nt = size(Ac)
+    m,  n  = size(A)
+    mb, nb = A.partitioning.blocksize
+    
+    mb != nb && error("Unequal block sizes are not supported: mb = $mb, nb = $nb")
+
+    ipiv = DVector(collect(1:min(m, n)), Blocks(mb))
+    ipivc = ipiv.chunks
+
+    max_piv_idx = zeros(Int,mt)
+    max_piv_val = zeros(T, mt)
+ 
+    Dagger.spawn_datadeps() do
+        for k in 1:min(mt, nt)
+            for p in 1:min(nb, m-(k-1)*nb, n-(k-1)*nb)
+                Dagger.@spawn searchmax_pivot!(Out(view(max_piv_idx, k:k)), Out(view(max_piv_val, k:k)), In(view(Ac[k,k],p:min(nb,m-(k-1)*nb),p:p)), p-1)
+                for i in k+1:mt
+                    Dagger.@spawn searchmax_pivot!(Out(view(max_piv_idx, i:i)), Out(view(max_piv_val, i:i)), In(view(Ac[i,k],:,p:p)))
+                end
+                Dagger.@spawn update_ipiv!(InOut(view(ipivc[k],p:p)), In(view(max_piv_idx, k:mt)), In(view(max_piv_val, k:mt)), k, nb)
+                for i in k:mt
+                    Dagger.@spawn swaprows_panel!(InOut(Ac[k, k]), InOut(Ac[i, k]), InOut(view(ipivc[k],p:p)), i, p, nb) 
+                end
+                Dagger.@spawn update_panel!(InOut(view(Ac[k,k],p+1:min(nb,m-(k-1)*nb),:)), In(Ac[k,k]), p)
+                for i in k+1:mt
+                    Dagger.@spawn update_panel!(InOut(Ac[i, k]), In(Ac[k,k]), p)
+                end
+
+            end
+            for j in Iterators.flatten((1:k-1, k+1:nt))
+                for i in k:mt
+                    Dagger.@spawn swaprows_trail!(InOut(Ac[k, j]), InOut(Ac[i, j]), In(ipivc[k]), i, mb) 
+                end
+            end
+            for j in k+1:nt
+                Dagger.@spawn BLAS.trsm!('L', 'L', 'N', 'U', zone, In(Ac[k, k]), InOut(Ac[k, j]))
+                for i in k+1:mt
+                    Dagger.@spawn BLAS.gemm!('N', 'N', mzone, In(Ac[i, k]), In(Ac[k, j]), zone, InOut(Ac[i, j]))
+                end
+            end
+        end
+    end
+
+    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)    
+end 
\ No newline at end of file
diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 7f6993a02..3127b5cf7 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -1,33 +1,37 @@
-@testset "$T" for T in (Float32, Float64, ComplexF32, ComplexF64)
+@testset "$T with $pivot" for T in (Float32, Float64, ComplexF32, ComplexF64), pivot in (NoPivot(), RowMaximum())
     A = rand(T, 128, 128)
     B = copy(A)
     DA = view(A, Blocks(64, 64))
 
     # Out-of-place
-    lu_A = lu(A, NoPivot())
-    lu_DA = lu(DA, NoPivot())
+    lu_A = lu(A, pivot)
+    lu_DA = lu(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32, ComplexF64)) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
-    @test lu_A.P ≈ lu_DA.P
-    @test lu_A.p ≈ lu_DA.p
+    if !(T in (ComplexF32, ComplexF64))
+        @test lu_A.P ≈ lu_DA.P
+        @test lu_A.p ≈ lu_DA.p
+    end
     # Check that lu did not modify A or DA
     @test A ≈ DA ≈ B
 
     # In-place
     A_copy = copy(A)
-    lu_A = lu!(A_copy, NoPivot())
-    lu_DA = lu!(DA, NoPivot())
+    lu_A = lu!(A_copy, pivot)
+    lu_DA = lu!(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32, ComplexF64)) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
-    @test lu_A.P ≈ lu_DA.P
-    @test lu_A.p ≈ lu_DA.p
+    if !(T in (ComplexF32, ComplexF64))
+        @test lu_A.P ≈ lu_DA.P
+        @test lu_A.p ≈ lu_DA.p
+    end
     # Check that changes propagated to A
     @test DA ≈ A
     @test !(B ≈ A)
-end
+end
\ No newline at end of file

From 1242052a9c3912688820700f5aa4035780f34f57 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Wed, 9 Jul 2025 02:33:54 +0530
Subject: [PATCH 02/22] Remove debug print statements from pivoting functions
 in LU factorization

---
 src/array/lu.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 979eb4648..4fd192d00 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -32,13 +32,11 @@ function searchmax_pivot!(piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}
     max_idx = argmax(abs.(A[:]))
     piv_idx[1] = offset+max_idx
     piv_val[1] = A[max_idx]
-    println("searchmax_pivot: ", piv_idx[1], "\n", abs(piv_val[1]))
 end
 
 function update_ipiv!(ipivl, piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}, k::Int, nb::Int) where T
     max_piv_idx = argmax(abs.(piv_val))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
-    println("update_ipiv: ", ipivl[1])
 end
 
 function swaprows_panel!(A::AbstractArray{T}, M::AbstractArray{T}, ipivl::AbstractVector{Int}, m::Int, p::Int, nb::Int) where T
@@ -46,7 +44,6 @@ function swaprows_panel!(A::AbstractArray{T}, M::AbstractArray{T}, ipivl::Abstra
     r = (ipivl[1]-1)%nb+1
     if m == q
         A[p,:], M[r,:] = M[r,:], A[p,:]
-        println("swaprows_panel: ", imag.(A[p,:]), "\n", imag.(M[r,:]))
     end
 end
 
@@ -62,7 +59,6 @@ function swaprows_trail!(A::AbstractArray{T}, M::AbstractArray{T}, ipiv::Abstrac
         r = (ipiv[p]-1)%nb+1
         if m == q
             A[p,:], M[r,:] = M[r,:], A[p,:]
-        println("swaprows_trail: ", imag.(A[p,:]), "\n", imag.(M[r,:]))
         end
     end
 end

From bb1620d93cae7a2d80bb6d19ab62c325d43abda5 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Wed, 9 Jul 2025 15:07:12 +0530
Subject: [PATCH 03/22] Refactor function signatures for improved clarity and
 consistency

---
 src/array/lu.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 4fd192d00..7939ef474 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -28,18 +28,18 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
     return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)
 end
 
-function searchmax_pivot!(piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}, A::AbstractArray{T}, offset::Int=0) where T
+function searchmax_pivot!(piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, A::AbstractMatrix{T}, offset::Int=0) where T
     max_idx = argmax(abs.(A[:]))
     piv_idx[1] = offset+max_idx
     piv_val[1] = A[max_idx]
 end
 
-function update_ipiv!(ipivl, piv_idx::AbstractArray{Int}, piv_val::AbstractArray{T}, k::Int, nb::Int) where T
+function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
     max_piv_idx = argmax(abs.(piv_val))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
 end
 
-function swaprows_panel!(A::AbstractArray{T}, M::AbstractArray{T}, ipivl::AbstractVector{Int}, m::Int, p::Int, nb::Int) where T
+function swaprows_panel!(A::AbstractMatrix{T}, M::AbstractMatrix{T}, ipivl::AbstractVector{Int}, m::Int, p::Int, nb::Int) where T
     q = div(ipivl[1]-1,nb) + 1
     r = (ipivl[1]-1)%nb+1
     if m == q
@@ -47,13 +47,13 @@ function swaprows_panel!(A::AbstractArray{T}, M::AbstractArray{T}, ipivl::Abstra
     end
 end
 
-function update_panel!(M::AbstractArray{T}, A::AbstractArray{T}, p::Int) where T
+function update_panel!(M::AbstractMatrix{T}, A::AbstractMatrix{T}, p::Int) where T
     Acinv = one(T) / A[p,p]  
     LinearAlgebra.BLAS.scal!(Acinv, view(M, :, p))
     LinearAlgebra.BLAS.ger!(-one(T), view(M, :, p), conj.(view(A, p, p+1:size(A,2))), view(M, :, p+1:size(M,2)))
 end
 
-function swaprows_trail!(A::AbstractArray{T}, M::AbstractArray{T}, ipiv::AbstractVector{Int}, m::Int, nb::Int) where T
+function swaprows_trail!(A::AbstractMatrix{T}, M::AbstractMatrix{T}, ipiv::AbstractVector{Int}, m::Int, nb::Int) where T
     for p in eachindex(ipiv)
         q = div(ipiv[p]-1,nb) + 1
         r = (ipiv[p]-1)%nb+1

From 11ccee700f42c41867fba6888f6c22404435a6c6 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Wed, 9 Jul 2025 22:32:42 +0530
Subject: [PATCH 04/22] Add singularity check for robust LU factorization

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

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 7939ef474..bc335d717 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -36,6 +36,8 @@ end
 
 function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
     max_piv_idx = argmax(abs.(piv_val))
+    max_piv_val = abs(piv_val[max_piv_idx])
+    isapprox(max_piv_val, zero(T); atol=eps(T)) && throw(SingularException(k))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
 end
 

From 11042f9ae725d3e6009f8cbed479a31258e4fcd9 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <41839847+AkhilAkkapelli@users.noreply.github.com>
Date: Thu, 10 Jul 2025 06:18:12 +0530
Subject: [PATCH 05/22] Fix lu.jl

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

diff --git a/src/array/lu.jl b/src/array/lu.jl
index bc335d717..3fe6bb7c1 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -37,7 +37,7 @@ end
 function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
     max_piv_idx = argmax(abs.(piv_val))
     max_piv_val = abs(piv_val[max_piv_idx])
-    isapprox(max_piv_val, zero(T); atol=eps(T)) && throw(SingularException(k))
+    isapprox(max_piv_val, zero(T); atol=eps(T)) && throw(LinearAlgebra.SingularException(k))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
 end
 
@@ -118,4 +118,4 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     end
 
     return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)    
-end 
\ No newline at end of file
+end 

From 794e55488ee6a6fc852841def0aea5a483908433 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Thu, 10 Jul 2025 19:14:43 +0530
Subject: [PATCH 06/22] Fix precision issue in singularity check

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

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 3fe6bb7c1..623349163 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -37,7 +37,7 @@ end
 function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
     max_piv_idx = argmax(abs.(piv_val))
     max_piv_val = abs(piv_val[max_piv_idx])
-    isapprox(max_piv_val, zero(T); atol=eps(T)) && throw(LinearAlgebra.SingularException(k))
+    isapprox(max_piv_val, zero(T); atol=eps(real(T))) && throw(LinearAlgebra.SingularException(k))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
 end
 

From afe96f74b13cd4f16028b3d79431398553e71a5f Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Thu, 10 Jul 2025 20:23:56 +0530
Subject: [PATCH 07/22] Fix pivot search by replacing argmax with BLAS.iamax

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

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 623349163..8ea1c7d6c 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -29,15 +29,16 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
 end
 
 function searchmax_pivot!(piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, A::AbstractMatrix{T}, offset::Int=0) where T
-    max_idx = argmax(abs.(A[:]))
+    max_idx = LinearAlgebra.BLAS.iamax(A[:])
     piv_idx[1] = offset+max_idx
     piv_val[1] = A[max_idx]
 end
 
 function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
-    max_piv_idx = argmax(abs.(piv_val))
-    max_piv_val = abs(piv_val[max_piv_idx])
-    isapprox(max_piv_val, zero(T); atol=eps(real(T))) && throw(LinearAlgebra.SingularException(k))
+    max_piv_idx = LinearAlgebra.BLAS.iamax(piv_val)
+    max_piv_val = piv_val[max_piv_idx]
+    abs_max_piv_val = max_piv_val isa Real ? abs(max_piv_val) : abs(real(max_piv_val)) + abs(imag(max_piv_val))
+    isapprox(abs_max_piv_val, zero(T); atol=eps(real(T))) && throw(LinearAlgebra.SingularException(k))
     ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
 end
 

From e636bbb9713afafd73c49a4394507a524824442b Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Thu, 10 Jul 2025 20:36:29 +0530
Subject: [PATCH 08/22] Refactor LU tests for consistency and clarity in pivot
 checks

---
 test/array/linalg/lu.jl | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 3127b5cf7..3ffe27d35 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -6,15 +6,13 @@
     # Out-of-place
     lu_A = lu(A, pivot)
     lu_DA = lu(DA, pivot)
-    @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32, ComplexF64)) # FIXME: NoPivot is unstable for FP32
+    @test lu_DA isa LU{T,DMatrix{T},DVector{Int}} 
+    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
-    if !(T in (ComplexF32, ComplexF64))
-        @test lu_A.P ≈ lu_DA.P
-        @test lu_A.p ≈ lu_DA.p
-    end
+    @test lu_A.P ≈ lu_DA.P
+    @test lu_A.p ≈ lu_DA.p
     # Check that lu did not modify A or DA
     @test A ≈ DA ≈ B
 
@@ -23,14 +21,12 @@
     lu_A = lu!(A_copy, pivot)
     lu_DA = lu!(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32, ComplexF64)) # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
-    if !(T in (ComplexF32, ComplexF64))
-        @test lu_A.P ≈ lu_DA.P
-        @test lu_A.p ≈ lu_DA.p
-    end
+    @test lu_A.P ≈ lu_DA.P
+    @test lu_A.p ≈ lu_DA.p
     # Check that changes propagated to A
     @test DA ≈ A
     @test !(B ≈ A)

From 337e59d733379f28b87958cafe9a8359fca843f5 Mon Sep 17 00:00:00 2001
From: Julian P Samaroo <jpsamaroo@gmail.com>
Date: Thu, 10 Jul 2025 14:00:45 -0700
Subject: [PATCH 09/22] Minor LU cleanup

---
 src/array/lu.jl         | 20 +++++++++++---------
 test/array/linalg/lu.jl |  2 +-
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 8ea1c7d6c..2abc9544e 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -7,6 +7,9 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
     mzone = -one(T)
     Ac = A.chunks
     mt, nt = size(Ac)
+    mb, nb = A.partitioning.blocksize
+
+    mb != nb && error("Unequal block sizes are not supported: mb = $mb, nb = $nb")
 
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
@@ -39,7 +42,7 @@ function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int},
     max_piv_val = piv_val[max_piv_idx]
     abs_max_piv_val = max_piv_val isa Real ? abs(max_piv_val) : abs(real(max_piv_val)) + abs(imag(max_piv_val))
     isapprox(abs_max_piv_val, zero(T); atol=eps(real(T))) && throw(LinearAlgebra.SingularException(k))
-    ipivl[1] = (max_piv_idx+k-2)*nb +  piv_idx[max_piv_idx]
+    ipivl[1] = (max_piv_idx+k-2)*nb + piv_idx[max_piv_idx]
 end
 
 function swaprows_panel!(A::AbstractMatrix{T}, M::AbstractMatrix{T}, ipivl::AbstractVector{Int}, m::Int, p::Int, nb::Int) where T
@@ -51,7 +54,7 @@ function swaprows_panel!(A::AbstractMatrix{T}, M::AbstractMatrix{T}, ipivl::Abst
 end
 
 function update_panel!(M::AbstractMatrix{T}, A::AbstractMatrix{T}, p::Int) where T
-    Acinv = one(T) / A[p,p]  
+    Acinv = one(T) / A[p,p]
     LinearAlgebra.BLAS.scal!(Acinv, view(M, :, p))
     LinearAlgebra.BLAS.ger!(-one(T), view(M, :, p), conj.(view(A, p, p+1:size(A,2))), view(M, :, p+1:size(M,2)))
 end
@@ -78,7 +81,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     mt, nt = size(Ac)
     m,  n  = size(A)
     mb, nb = A.partitioning.blocksize
-    
+
     mb != nb && error("Unequal block sizes are not supported: mb = $mb, nb = $nb")
 
     ipiv = DVector(collect(1:min(m, n)), Blocks(mb))
@@ -86,7 +89,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
 
     max_piv_idx = zeros(Int,mt)
     max_piv_val = zeros(T, mt)
- 
+
     Dagger.spawn_datadeps() do
         for k in 1:min(mt, nt)
             for p in 1:min(nb, m-(k-1)*nb, n-(k-1)*nb)
@@ -96,17 +99,16 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
                 end
                 Dagger.@spawn update_ipiv!(InOut(view(ipivc[k],p:p)), In(view(max_piv_idx, k:mt)), In(view(max_piv_val, k:mt)), k, nb)
                 for i in k:mt
-                    Dagger.@spawn swaprows_panel!(InOut(Ac[k, k]), InOut(Ac[i, k]), InOut(view(ipivc[k],p:p)), i, p, nb) 
+                    Dagger.@spawn swaprows_panel!(InOut(Ac[k, k]), InOut(Ac[i, k]), In(view(ipivc[k],p:p)), i, p, nb)
                 end
                 Dagger.@spawn update_panel!(InOut(view(Ac[k,k],p+1:min(nb,m-(k-1)*nb),:)), In(Ac[k,k]), p)
                 for i in k+1:mt
                     Dagger.@spawn update_panel!(InOut(Ac[i, k]), In(Ac[k,k]), p)
                 end
-
             end
             for j in Iterators.flatten((1:k-1, k+1:nt))
                 for i in k:mt
-                    Dagger.@spawn swaprows_trail!(InOut(Ac[k, j]), InOut(Ac[i, j]), In(ipivc[k]), i, mb) 
+                    Dagger.@spawn swaprows_trail!(InOut(Ac[k, j]), InOut(Ac[i, j]), In(ipivc[k]), i, mb)
                 end
             end
             for j in k+1:nt
@@ -118,5 +120,5 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
         end
     end
 
-    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)    
-end 
+    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)
+end
diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 3ffe27d35..4e558572b 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -6,7 +6,7 @@
     # Out-of-place
     lu_A = lu(A, pivot)
     lu_DA = lu(DA, pivot)
-    @test lu_DA isa LU{T,DMatrix{T},DVector{Int}} 
+    @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
     if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U

From ffdf186e40dbe589f646e136449ba4457d87044d Mon Sep 17 00:00:00 2001
From: Julian P Samaroo <jpsamaroo@gmail.com>
Date: Thu, 10 Jul 2025 14:04:33 -0700
Subject: [PATCH 10/22] test/array/linalg/lu: Test for invalid block sizes

---
 src/array/lu.jl         | 4 ++--
 test/array/linalg/lu.jl | 4 ++++
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 2abc9544e..4c58d9c2e 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -9,7 +9,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
     mt, nt = size(Ac)
     mb, nb = A.partitioning.blocksize
 
-    mb != nb && error("Unequal block sizes are not supported: mb = $mb, nb = $nb")
+    mb != nb && throw(ArgumentError("Unequal block sizes are not supported: mb = $mb, nb = $nb"))
 
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
@@ -82,7 +82,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     m,  n  = size(A)
     mb, nb = A.partitioning.blocksize
 
-    mb != nb && error("Unequal block sizes are not supported: mb = $mb, nb = $nb")
+    mb != nb && throw(ArgumentError("Unequal block sizes are not supported: mb = $mb, nb = $nb"))
 
     ipiv = DVector(collect(1:min(m, n)), Blocks(mb))
     ipivc = ipiv.chunks
diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 4e558572b..75bacb17b 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -30,4 +30,8 @@
     # Check that changes propagated to A
     @test DA ≈ A
     @test !(B ≈ A)
+
+    # Non-square block sizes are not supported
+    @test_throws ArgumentError lu(rand(Blocks(64, 32), T, 128, 128), pivot)
+    @test_throws ArgumentError lu!(rand(Blocks(64, 32), T, 128, 128), pivot)
 end
\ No newline at end of file

From 67fe600dee9173dace6f3d41d6e80f29f91b6803 Mon Sep 17 00:00:00 2001
From: Julian P Samaroo <jpsamaroo@gmail.com>
Date: Thu, 10 Jul 2025 14:25:22 -0700
Subject: [PATCH 11/22] docs: Add LU RowMaximum to supported ops

---
 docs/src/darray.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/darray.md b/docs/src/darray.md
index 956436b6d..8bca67488 100644
--- a/docs/src/darray.md
+++ b/docs/src/darray.md
@@ -694,7 +694,7 @@ From `LinearAlgebra`:
 - `*` (Out-of-place Matrix-(Matrix/Vector) multiply)
 - `mul!` (In-place Matrix-Matrix multiply)
 - `cholesky`/`cholesky!` (In-place/Out-of-place Cholesky factorization)
-- `lu`/`lu!` (In-place/Out-of-place LU factorization (`NoPivot` only))
+- `lu`/`lu!` (In-place/Out-of-place LU factorization (`NoPivot` and `RowMaximum` only))
 
 From `AbstractFFTs`:
 - `fft`/`fft!`

From adf11fcbc8771c1b0dbdef62421e740832973fe3 Mon Sep 17 00:00:00 2001
From: Julian P Samaroo <jpsamaroo@gmail.com>
Date: Thu, 10 Jul 2025 14:25:40 -0700
Subject: [PATCH 12/22] docs: Add SparseArrays sprand to supported ops

---
 docs/src/darray.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/docs/src/darray.md b/docs/src/darray.md
index 8bca67488..933994864 100644
--- a/docs/src/darray.md
+++ b/docs/src/darray.md
@@ -684,6 +684,9 @@ From `Base`:
 From `Random`:
 - `rand!`/`randn!`
 
+From `SparseArrays`:
+- `sprand`
+
 From `Statistics`:
 - `mean`
 - `var`

From 330a9e3011763437a458c77ca4bb68902822d2bb Mon Sep 17 00:00:00 2001
From: Julian P Samaroo <jpsamaroo@gmail.com>
Date: Thu, 10 Jul 2025 14:26:07 -0700
Subject: [PATCH 13/22] test/array/linalg/lu: Fix skipped tests

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 75bacb17b..cbf98dc87 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -7,7 +7,7 @@
     lu_A = lu(A, pivot)
     lu_DA = lu(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
@@ -21,7 +21,7 @@
     lu_A = lu!(A_copy, pivot)
     lu_DA = lu!(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end

From f06b594cda5520bdaaab1fefbb6f22bc42bc17b4 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 08:10:20 +0530
Subject: [PATCH 14/22] Fix logical condition for NoPivot stability check

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index cbf98dc87..75bacb17b 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -7,7 +7,7 @@
     lu_A = lu(A, pivot)
     lu_DA = lu(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
@@ -21,7 +21,7 @@
     lu_A = lu!(A_copy, pivot)
     lu_DA = lu!(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end

From ffb03bb334c419f4769543330641a12c5b4eac20 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 10:20:13 +0530
Subject: [PATCH 15/22] Enhance LU decomposition functions to support
 allowsingular parameter and improve error handling

---
 src/array/lu.jl         | 44 ++++++++++++++++++++++++++++++-----------
 test/array/linalg/lu.jl |  3 +++
 2 files changed, 35 insertions(+), 12 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 4c58d9c2e..dac2a2ff6 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -1,19 +1,28 @@
-function LinearAlgebra.lu(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=true) where T
+LinearAlgebra.lu(A::DMatrix{T}, pivot::Union{LinearAlgebra.RowMaximum,LinearAlgebra.NoPivot} = LinearAlgebra.RowMaximum(); check::Bool=true, allowsingular::Bool=false) where {T<:LinearAlgebra.BlasFloat} = LinearAlgebra.lu(A, pivot; check=check, allowsingular=allowsingular)
+
+LinearAlgebra.lu!(A::DMatrix{T}, pivot::Union{LinearAlgebra.RowMaximum,LinearAlgebra.NoPivot} = LinearAlgebra.RowMaximum(); check::Bool=true, allowsingular::Bool=false) where {T<:LinearAlgebra.BlasFloat} = LinearAlgebra.lu(A, pivot; check=check, allowsingular=allowsingular)
+
+function LinearAlgebra.lu(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool = true, allowsingular::Bool = false) where {T<:LinearAlgebra.BlasFloat}
     A_copy = LinearAlgebra._lucopy(A, LinearAlgebra.lutype(T))
     return LinearAlgebra.lu!(A_copy, LinearAlgebra.NoPivot(); check=check)
 end
-function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=true) where T
+function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool = true, allowsingular::Bool = false) where {T<:LinearAlgebra.BlasFloat}
+  
+    check && LinearAlgebra.LAPACK.chkfinite(A)
+
     zone = one(T)
     mzone = -one(T)
     Ac = A.chunks
     mt, nt = size(Ac)
     mb, nb = A.partitioning.blocksize
 
+    info = 0
+
     mb != nb && throw(ArgumentError("Unequal block sizes are not supported: mb = $mb, nb = $nb"))
 
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
-            Dagger.@spawn LinearAlgebra.generic_lufact!(InOut(Ac[k, k]), LinearAlgebra.NoPivot(); check)
+            Dagger.@spawn LinearAlgebra.generic_lufact!(InOut(Ac[k, k]), LinearAlgebra.NoPivot(); check, allowsingular)
             for m in range(k+1, mt)
                 Dagger.@spawn BLAS.trsm!('R', 'U', 'N', 'N', zone, In(Ac[k, k]), InOut(Ac[m, k]))
             end
@@ -28,7 +37,9 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool=t
 
     ipiv = DVector([i for i in 1:min(size(A)...)])
 
-    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)
+    check && LinearAlgebra._check_lu_success(info, allowsingular)
+
+    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, info)
 end
 
 function searchmax_pivot!(piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, A::AbstractMatrix{T}, offset::Int=0) where T
@@ -37,11 +48,13 @@ function searchmax_pivot!(piv_idx::AbstractVector{Int}, piv_val::AbstractVector{
     piv_val[1] = A[max_idx]
 end
 
-function update_ipiv!(ipivl::AbstractVector{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
+function update_ipiv!(ipivl::AbstractVector{Int}, info::Ref{Int}, piv_idx::AbstractVector{Int}, piv_val::AbstractVector{T}, k::Int, nb::Int) where T
     max_piv_idx = LinearAlgebra.BLAS.iamax(piv_val)
     max_piv_val = piv_val[max_piv_idx]
-    abs_max_piv_val = max_piv_val isa Real ? abs(max_piv_val) : abs(real(max_piv_val)) + abs(imag(max_piv_val))
-    isapprox(abs_max_piv_val, zero(T); atol=eps(real(T))) && throw(LinearAlgebra.SingularException(k))
+    abs_max_piv_val = max_piv_val isa Real ? abs(max_piv_val) : abs(real(max_piv_val)) + abs(imag(max_piv_val))    
+    if isapprox(abs_max_piv_val, zero(T); atol=eps(real(T)))
+        info[] = k
+    end
     ipivl[1] = (max_piv_idx+k-2)*nb + piv_idx[max_piv_idx]
 end
 
@@ -69,11 +82,14 @@ function swaprows_trail!(A::AbstractMatrix{T}, M::AbstractMatrix{T}, ipiv::Abstr
     end
 end
 
-function LinearAlgebra.lu(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool=true) where T
+function LinearAlgebra.lu(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool = true, allowsingular::Bool = false) where {T<:LinearAlgebra.BlasFloat}
     A_copy = LinearAlgebra._lucopy(A, LinearAlgebra.lutype(T))
-    return LinearAlgebra.lu!(A_copy, LinearAlgebra.RowMaximum(); check=check)
+    return LinearAlgebra.lu!(A_copy, LinearAlgebra.RowMaximum(); check=check, allowsingular=allowsingular)
 end
-function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool=true) where T
+function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Bool = true, allowsingular::Bool = false) where {T<:LinearAlgebra.BlasFloat}
+    
+    check && LinearAlgebra.LAPACK.chkfinite(A)
+
     zone = one(T)
     mzone = -one(T)
 
@@ -87,6 +103,8 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     ipiv = DVector(collect(1:min(m, n)), Blocks(mb))
     ipivc = ipiv.chunks
 
+    info = Ref(0)
+
     max_piv_idx = zeros(Int,mt)
     max_piv_val = zeros(T, mt)
 
@@ -97,7 +115,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
                 for i in k+1:mt
                     Dagger.@spawn searchmax_pivot!(Out(view(max_piv_idx, i:i)), Out(view(max_piv_val, i:i)), In(view(Ac[i,k],:,p:p)))
                 end
-                Dagger.@spawn update_ipiv!(InOut(view(ipivc[k],p:p)), In(view(max_piv_idx, k:mt)), In(view(max_piv_val, k:mt)), k, nb)
+                Dagger.@spawn update_ipiv!(InOut(view(ipivc[k],p:p)), InOut(info), In(view(max_piv_idx, k:mt)), In(view(max_piv_val, k:mt)), k, nb)
                 for i in k:mt
                     Dagger.@spawn swaprows_panel!(InOut(Ac[k, k]), InOut(Ac[i, k]), In(view(ipivc[k],p:p)), i, p, nb)
                 end
@@ -120,5 +138,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
         end
     end
 
-    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, 0)
+    check && LinearAlgebra._check_lu_success(info[], allowsingular)
+
+    return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, info[])
 end
diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 75bacb17b..871d46169 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -34,4 +34,7 @@
     # Non-square block sizes are not supported
     @test_throws ArgumentError lu(rand(Blocks(64, 32), T, 128, 128), pivot)
     @test_throws ArgumentError lu!(rand(Blocks(64, 32), T, 128, 128), pivot)
+
+        # Singular Values
+    pivot == RowMaximum() || @test_throws LinearAlgebra.SingularException lu(rand(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
 end
\ No newline at end of file

From c1471347e19a1709ac6c231906e7f59526ed6938 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 10:22:45 +0530
Subject: [PATCH 16/22] Fix singular value test to use ones instead of rand for
 RowMaximum pivot

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 871d46169..a8ab9d059 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -36,5 +36,5 @@
     @test_throws ArgumentError lu!(rand(Blocks(64, 32), T, 128, 128), pivot)
 
         # Singular Values
-    pivot == RowMaximum() || @test_throws LinearAlgebra.SingularException lu(rand(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
+    pivot == RowMaximum() || @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
 end
\ No newline at end of file

From d116d51857999f3b30019e995d02442222d24165 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 11:11:00 +0530
Subject: [PATCH 17/22] Refactor LU decomposition functions to handle unequal
 block sizes and improve buffer management

---
 src/array/lu.jl | 27 ++++++++++++++++++++-------
 1 file changed, 20 insertions(+), 7 deletions(-)

diff --git a/src/array/lu.jl b/src/array/lu.jl
index dac2a2ff6..5672a8ab3 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -12,14 +12,21 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool =
 
     zone = one(T)
     mzone = -one(T)
+
+    mb, nb = A.partitioning.blocksize
+
+    if mb != nb 
+        mb = nb = min(mb, nb)
+        A = maybe_copy_buffered(A => Blocks(nb, nb)) do A
+            A
+        end
+    end
+
     Ac = A.chunks
     mt, nt = size(Ac)
-    mb, nb = A.partitioning.blocksize
 
     info = 0
 
-    mb != nb && throw(ArgumentError("Unequal block sizes are not supported: mb = $mb, nb = $nb"))
-
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
             Dagger.@spawn LinearAlgebra.generic_lufact!(InOut(Ac[k, k]), LinearAlgebra.NoPivot(); check, allowsingular)
@@ -93,14 +100,20 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     zone = one(T)
     mzone = -one(T)
 
+    mb, nb = A.partitioning.blocksize
+
+    if mb != nb 
+        mb = nb = min(mb, nb)
+        A = maybe_copy_buffered(A => Blocks(nb, nb)) do A
+            A
+        end
+    end
+
     Ac = A.chunks
     mt, nt = size(Ac)
     m,  n  = size(A)
-    mb, nb = A.partitioning.blocksize
-
-    mb != nb && throw(ArgumentError("Unequal block sizes are not supported: mb = $mb, nb = $nb"))
 
-    ipiv = DVector(collect(1:min(m, n)), Blocks(mb))
+    ipiv = DVector(collect(1:min(m, n)), Blocks(nb))
     ipivc = ipiv.chunks
 
     info = Ref(0)

From ce2480ee0f5a41cd125b98f3d91e01ab916d25ff Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 11:14:22 +0530
Subject: [PATCH 18/22] Update LU tests to support non-square block sizes and
 improve singular value checks

---
 test/array/linalg/lu.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index a8ab9d059..15f72c2bc 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -31,10 +31,10 @@
     @test DA ≈ A
     @test !(B ≈ A)
 
-    # Non-square block sizes are not supported
-    @test_throws ArgumentError lu(rand(Blocks(64, 32), T, 128, 128), pivot)
-    @test_throws ArgumentError lu!(rand(Blocks(64, 32), T, 128, 128), pivot)
+    # Non-square block sizes
+    @test lu(rand(Blocks(64, 32), T, 128, 128), pivot) isa LU{T,DMatrix{T},DVector{Int}}
+    @test lu!(rand(Blocks(64, 32), T, 128, 128), pivot) isa LU{T,DMatrix{T},DVector{Int}}
 
-        # Singular Values
+    # Singular Values
     pivot == RowMaximum() || @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
 end
\ No newline at end of file

From 09b38825750d9a89d0eb5c4be7d4d6cf5609c282 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 11:34:25 +0530
Subject: [PATCH 19/22] Fix logical condition for singular value exception
 check in LU tests

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 15f72c2bc..5ca1832ed 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -36,5 +36,5 @@
     @test lu!(rand(Blocks(64, 32), T, 128, 128), pivot) isa LU{T,DMatrix{T},DVector{Int}}
 
     # Singular Values
-    pivot == RowMaximum() || @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
+    @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
 end
\ No newline at end of file

From d508e965b5586e66949a35f2fd5871100bf0f0fe Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 11:36:36 +0530
Subject: [PATCH 20/22] Fix comment for singular value exception check in LU
 tests

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 5ca1832ed..7f941146e 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -36,5 +36,5 @@
     @test lu!(rand(Blocks(64, 32), T, 128, 128), pivot) isa LU{T,DMatrix{T},DVector{Int}}
 
     # Singular Values
-    @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs Singular Exception Check
+    @test_throws LinearAlgebra.SingularException lu(ones(Blocks(64,64), T, 128, 128)) # FIXME: NoPivot needs to handle info
 end
\ No newline at end of file

From c7175026b8a1915bd553ef1b33a3fdc326dc78df Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Fri, 11 Jul 2025 21:59:13 +0530
Subject: [PATCH 21/22] Fix logical condition for NoPivot check in LU tests

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

diff --git a/test/array/linalg/lu.jl b/test/array/linalg/lu.jl
index 7f941146e..1d9b99ce4 100644
--- a/test/array/linalg/lu.jl
+++ b/test/array/linalg/lu.jl
@@ -7,7 +7,7 @@
     lu_A = lu(A, pivot)
     lu_DA = lu(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end
@@ -21,7 +21,7 @@
     lu_A = lu!(A_copy, pivot)
     lu_DA = lu!(DA, pivot)
     @test lu_DA isa LU{T,DMatrix{T},DVector{Int}}
-    if !(T in (Float32, ComplexF32)) && pivot == NoPivot() # FIXME: NoPivot is unstable for FP32
+    if !(T in (Float32, ComplexF32) && pivot == NoPivot()) # FIXME: NoPivot is unstable for FP32
         @test lu_A.L ≈ lu_DA.L
         @test lu_A.U ≈ lu_DA.U
     end

From 2d2e1224decc5d058353ade8cd0d31e4ed324e90 Mon Sep 17 00:00:00 2001
From: Akhil Akkapelli <akhilakkapelli@gmail.com>
Date: Sat, 12 Jul 2025 08:11:33 +0530
Subject: [PATCH 22/22] Fix argument passing for NoPivot in generic_lufact!
 calls in LU decomposition

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

diff --git a/src/array/lu.jl b/src/array/lu.jl
index 5672a8ab3..b00301e88 100644
--- a/src/array/lu.jl
+++ b/src/array/lu.jl
@@ -29,7 +29,7 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.NoPivot; check::Bool =
 
     Dagger.spawn_datadeps() do
         for k in range(1, min(mt, nt))
-            Dagger.@spawn LinearAlgebra.generic_lufact!(InOut(Ac[k, k]), LinearAlgebra.NoPivot(); check, allowsingular)
+            Dagger.@spawn LinearAlgebra.generic_lufact!(InOut(Ac[k, k]), LinearAlgebra.NoPivot(); check=check, allowsingular=allowsingular)
             for m in range(k+1, mt)
                 Dagger.@spawn BLAS.trsm!('R', 'U', 'N', 'N', zone, In(Ac[k, k]), InOut(Ac[m, k]))
             end
@@ -154,4 +154,4 @@ function LinearAlgebra.lu!(A::DMatrix{T}, ::LinearAlgebra.RowMaximum; check::Boo
     check && LinearAlgebra._check_lu_success(info[], allowsingular)
 
     return LinearAlgebra.LU{T,DMatrix{T},DVector{Int}}(A, ipiv, info[])
-end
+end
\ No newline at end of file