Merge branch '40-bug-lu-may-have-0-pivots-if-the-whole-matrix-is-abstol' into 'main'

Resolve "Bug: lu may have 0 pivots if the whole matrix is < abstol"

Closes #40 and #41

See merge request tensors4fields/TensorCrossInterpolation.jl!77

Author: Marc Ritter

Project.toml

@@ -10,8 +10,8 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 
 [compat]
 EllipsisNotation = "1"
-julia = "1.6"
 QuadGK = "2.9"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

src/matrixlu.jl

@@ -136,7 +136,8 @@ function _optimizerrlu!(
         k = lu.npivot + 1
         newpivot = submatrixargmax(abs2, A, k)
         lu.error = abs(A[newpivot...])
-        if abs(lu.error) < reltol * maxerror || abs(lu.error) < abstol
+        # Add at least 1 pivot to get a well-defined L * U
+        if (abs(lu.error) < reltol * maxerror || abs(lu.error) < abstol) && lu.npivot > 0
             break
         end
         maxerror = max(maxerror, lu.error)

test/test_matrixlu.jl

@@ -192,4 +192,20 @@ using LinearAlgebra
         lu3 = TCI.rrlu(A, abstol=0.0)
         @test TCI.lastpivoterror(lu3) == 0.0
     end
+
+    @testset "LU for matrices with very small absolute values" begin
+        A = 1e-13 * [
+            0.585383 0.124568 0.352426 0.573507
+            0.865875 0.600153 0.727443 0.902388
+            0.913477 0.954081 0.116965 0.817
+            0.985918 0.516114 0.600366 0.0200085
+        ]
+
+        lu = TCI.rrlu(A, abstol=1e-3)
+        @test TCI.npivots(lu) == 1
+        @test length(TCI.pivoterrors(lu)) > 0
+        @test TCI.lastpivoterror(lu) > 0
+        @test size(lu) == size(A)
+        @test maximum(abs.(TCI.left(lu) * TCI.right(lu) .- A)) < 1e-3
+    end
 end