Fallback for weird LAPACK error (#43)

* Fallback for weird LAPACK error

* Fix format

Author: blegat

src/schur.jl

@@ -8,13 +8,26 @@ function condition_number(sf::Schur, I)
     for i in I
         select[i] = 1
     end
-    return LinearAlgebra.LAPACK.trsen!(
-        'E',
-        'N',
-        select,
-        copy(sf.T),
-        copy(sf.Z),
-    )[4]
+    try
+        return LinearAlgebra.LAPACK.trsen!(
+            'E',
+            'N',
+            select,
+            copy(sf.T),
+            copy(sf.Z),
+        )[4]
+    catch err
+        if err isa LinearAlgebra.LAPACKException
+            ε = eps(real(eltype(sf.values)))
+            @warn(
+                "`LAPACK.trsen!` throwed an exception for `$(I)` so using default tolerance `$ε`"
+            )
+            # Not sure why this happens, see `test/schur` for an example
+            return ε
+        else
+            rethrow(err)
+        end
+    end
 end
 
 # Manocha, D. & Demmel, J. Algorithms for intersecting parametric and algebraic curves II: multiple intersections

test/commutativetests.jl

@@ -1,4 +1,5 @@
 include("macro.jl")
 include("groebner.jl")
+include("schur.jl")
 include("solve.jl")
 include("variety.jl")

test/schur.jl

@@ -0,0 +1,77 @@
+module TestSchur
+
+using Test
+
+include("utils.jl")
+
+function test_lapack_exception()
+    A = [
+        -1.96262e-16 0 0 0
+        -1.17757e-16 0 0 0
+        0 0 0 0
+        0 0 0 0
+    ]
+    B = [
+        0 0 0 -3.92523e-17
+        1 0 0 -1.17757e-16
+        0 1 0 0
+        0 0 1 0
+    ]
+    λ = [0.285173013664907, 0.714826986335093]
+    # Throws `LAPACKException(1)`
+    el = SS._solve_multiplication_matrices([A, B], λ, schur_solver)
+    testelements(el, [[0, 0]]; atol = 1e-10)
+    return
+end
+
+function _test_cgt96_e51(solver)
+    ɛ = 1e-4
+    Iɛ = [
+        1-ɛ 0
+        0 1+ɛ
+    ]
+    J = [
+        0 1
+        1 0
+    ]
+    Z = zeros(2, 2)
+    A = [
+        Iɛ Z
+        Z J
+    ]
+    B = [
+        J Z
+        Z Iɛ
+    ]
+    α = 0.219
+    testelements(
+        SS._solve_multiplication_matrices([A, B], [α, 1 - α], solver),
+        [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]];
+        rtol = 1e-7,
+    )
+    return
+end
+
+# Example 5.1 of CGT97
+function test_cgt96_e51()
+    @testset "Schur" begin
+        _test_cgt96_e51(schur_solver)
+    end
+    #    @testset "Newton" begin
+    #        _test_cgt96_e51(newton_solver)
+    #    end
+end
+
+function runtests(args...)
+    for name in names(@__MODULE__; all = true)
+        if startswith("$name", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)(args...)
+            end
+        end
+    end
+end
+
+end
+
+TestSchur.runtests()

test/solve.jl

@@ -30,34 +30,7 @@ end
     @test sort(sort.(clusters); by = first) == [[1, 5, 6], [2]]
 end
 
-function testelements(X, Y; atol = Base.rtoldefault(Float64), kwargs...)
-    @test length(X) == length(Y)
-    for y in Y
-        @test any(x -> isapprox(x, y; atol = atol, kwargs...), X)
-    end
-end
-function testelementstypes(X, Y; kwargs...)
-    testelements(X, Y; kwargs...)
-    for T in [Rational{Int}, Float64]
-        if X isa FixedVariablesSet
-            U = T
-        else
-            U = float(T)
-        end
-        V = similar(X, T)
-        testelements(V, Y; kwargs...)
-        @test eltype(V) == Vector{U}
-        @test collect(V) isa Vector{Vector{U}}
-    end
-end
-
-# We use a fixed RNG in the tests to decrease nondeterminism. There is still nondeterminism in LAPACK though
-using Random
-schur_solver = ReorderedSchurMultiplicationMatricesSolver(
-    sqrt(eps(Float64)),
-    MersenneTwister(0),
-)
-newton_solver = NewtonTypeDiagonalization()
+include("utils.jl")
 
 function zero_dimensional_ideal(solver)
     Mod.@polyvar x y z
@@ -141,46 +114,6 @@ end
     projective_zero_dimensional_ideal(newton_solver)
 end
 
-function cgt96_e51(solver)
-    ɛ = 1e-4
-    Iɛ = [
-        1-ɛ 0
-        0 1+ɛ
-    ]
-    J = [
-        0 1
-        1 0
-    ]
-    Z = zeros(2, 2)
-    A = [
-        Iɛ Z
-        Z J
-    ]
-    B = [
-        J Z
-        Z Iɛ
-    ]
-    α = 0.219
-    return testelements(
-        SemialgebraicSets._solve_multiplication_matrices(
-            [A, B],
-            [α, 1 - α],
-            solver,
-        ),
-        [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]];
-        rtol = 1e-7,
-    )
-end
-
-@testset "Example 5.1 of CGT97" begin
-    @testset "Schur" begin
-        cgt96_e51(schur_solver)
-    end
-    #    @testset "Newton" begin
-    #        cgt96_e51(newton_solver)
-    #    end
-end
-
 function md95_e43(solver)
     Mod.@polyvar x y
     V = @set x^2 + 4y^4 == 4y^2 && 3x^4 + y^2 == 5x^3 solver

test/utils.jl

@@ -0,0 +1,30 @@
+import SemialgebraicSets as SS
+
+function testelements(X, Y; atol = Base.rtoldefault(Float64), kwargs...)
+    @test length(X) == length(Y)
+    for y in Y
+        @test any(x -> isapprox(x, y; atol = atol, kwargs...), X)
+    end
+end
+function testelementstypes(X, Y; kwargs...)
+    testelements(X, Y; kwargs...)
+    for T in [Rational{Int}, Float64]
+        if X isa FixedVariablesSet
+            U = T
+        else
+            U = float(T)
+        end
+        V = similar(X, T)
+        testelements(V, Y; kwargs...)
+        @test eltype(V) == Vector{U}
+        @test collect(V) isa Vector{Vector{U}}
+    end
+end
+
+# We use a fixed RNG in the tests to decrease nondeterminism. There is still nondeterminism in LAPACK though
+using Random
+schur_solver = SS.ReorderedSchurMultiplicationMatricesSolver(
+    sqrt(eps(Float64)),
+    MersenneTwister(0),
+)
+newton_solver = SS.NewtonTypeDiagonalization()