From 0de6a346daec6facc9aae36cce3308f95483756f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Wed, 29 May 2024 14:30:56 +0200
Subject: [PATCH] Reinstate support SVector

---
 src/linear_eq.jl | 10 ++++++----
 src/region.jl    |  8 ++++++++
 test/roots.jl    |  2 --
 test/runtests.jl |  3 +++
 test/svectors.jl | 19 +++++++++++++++++++
 5 files changed, 36 insertions(+), 6 deletions(-)
 create mode 100644 test/svectors.jl

diff --git a/src/linear_eq.jl b/src/linear_eq.jl
index 790330b..26c43fd 100644
--- a/src/linear_eq.jl
+++ b/src/linear_eq.jl
@@ -75,13 +75,15 @@ function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::Abstra
     x
 end
 
-function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray ; precondition=true)
+function gauss_elimination_interval(A0::AbstractMatrix, b0::AbstractArray ; precondition=true)
     if precondition
-        A, b = preconditioner(A, b)
+        A0, b0 = preconditioner(A0, b0)
     end
 
-    A = copy(A)
-    b = copy(b)
+    A = similar(A0)
+    A .= A0
+    b = similar(b0)
+    b .= b0
     n = size(A, 1)
 
     p = similar(b)
diff --git a/src/region.jl b/src/region.jl
index b166032..31a3c9e 100644
--- a/src/region.jl
+++ b/src/region.jl
@@ -38,5 +38,13 @@ function bisect_region(X::AbstractVector, α)
     return X1, X2
 end
 
+function bisect_region(X::SVector{N, T}, α) where {N, T}
+    i = argmax(diam.(X))
+    x1, x2 = bisect_region(X[i], α)
+    X1 = SVector{N, T}(k == i ? x1 : X[k] for k in 1:N)
+    X2 = SVector{N, T}(k == i ? x2 : X[k] for k in 1:N)
+    return X1, X2
+end
+
 istrivial(X::Interval) = decoration(X) <= trv
 istrivial(X::AbstractVector) = any(istrivial.(X))
\ No newline at end of file
diff --git a/test/roots.jl b/test/roots.jl
index d95741c..913a078 100644
--- a/test/roots.jl
+++ b/test/roots.jl
@@ -29,8 +29,6 @@ function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10)
     @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol
 end
 
-newtonlike_methods = [Newton, Krawczyk]
-
 @testset "1D roots" begin
     # Default
     rts = roots(sin, interval(-5, 5))
diff --git a/test/runtests.jl b/test/runtests.jl
index c9c0aff..fdddf0e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,10 @@ using IntervalRootFinding
 using IntervalArithmetic.Symbols
 using Test
 
+newtonlike_methods = [Newton, Krawczyk]
+
 include("roots.jl")
+include("svectors.jl")
 include("test_smiley.jl")
 include("linear_eq.jl")
 
diff --git a/test/svectors.jl b/test/svectors.jl
new file mode 100644
index 0000000..5702a2c
--- /dev/null
+++ b/test/svectors.jl
@@ -0,0 +1,19 @@
+@testset "Static vectors" begin
+    f(xx) = [xx[1]^2 - 1, xx[2]^2 - 2]
+    g(xx) = SVector(xx[1]^2 - 1, xx[2]^2 - 2)
+
+    X = [interval(-5, 5), interval(-5, 5)]
+    S = SVector(interval(-5, 5), interval(-5, 5))
+
+    @testset for contractor in newtonlike_methods
+        rts = roots(f, X ; contractor)
+        rts2 = roots(g, S ; contractor)
+
+        for (rt, rt2) in zip(rts, rts2)
+            @test isequal_interval(
+                root_region(rt),
+                root_region(rt2)
+            )
+        end
+    end
+end
\ No newline at end of file