Add linear hull for solving system of linear equations

Author: eeshan9815

perf/linear_eq.jl

@@ -43,13 +43,14 @@ end
 
 function benchmark_elimination(max=10)
     df = DataFrame()
-    df[:Method] = ["Gauss Elimination", "Base.\\"]
+    df[:Method] = ["Gauss Elimination", "Base.\\", "Linear Hull"]
     for n in 1:max
         A, mA, sA = randMat(n)
         b, mb, sb = randVec(n)
         t1 = @belapsed gauss_elimination_interval($A, $b)
         t2 = @belapsed gauss_elimination_interval1($A, $b)
-        df[Symbol("$n")] = [t1, t2]
+        t3 = @belapsed linear_hull($A, $b)
+        df[Symbol("$n")] = [t1, t2, t3]
     end
     a = []
     for i in 1:max

perf/linear_eq_tabular.txt

@@ -39,11 +39,15 @@ Gauss Seidel
 
 Gauss Elimination
 
-5×3 DataFrames.DataFrame
-│ Row │ n │ Base.\\    │ Gauss Elimination │
-├─────┼───┼────────────┼───────────────────┤
-│ 1   │ 1 │ 1.84e-6    │ 8.127e-6          │
-│ 2   │ 2 │ 3.1615e-6  │ 1.0029e-5         │
-│ 3   │ 3 │ 4.53986e-6 │ 1.1211e-5         │
-│ 4   │ 4 │ 6.8655e-6  │ 1.4342e-5         │
-│ 5   │ 5 │ 1.0773e-5  │ 1.7725e-5         │
+│ Row │ n  │ Base.\\    │ Gauss Elimination │ Linear Hull │
+├─────┼────┼────────────┼───────────────────┼─────────────┤
+│ 1   │ 1  │ 1.9415e-6  │ 8.40867e-6        │ 3.94457e-6  │
+│ 2   │ 10 │ 7.1184e-5  │ 7.8812e-5         │ 4.6999e-5   │
+│ 3   │ 2  │ 3.34387e-6 │ 1.012e-5          │ 5.30567e-6  │
+│ 4   │ 3  │ 4.548e-6   │ 1.0955e-5         │ 6.5832e-6   │
+│ 5   │ 4  │ 7.256e-6   │ 1.3845e-5         │ 1.1115e-5   │
+│ 6   │ 5  │ 1.0878e-5  │ 1.7847e-5         │ 1.4564e-5   │
+│ 7   │ 6  │ 1.5747e-5  │ 2.1526e-5         │ 1.9468e-5   │
+│ 8   │ 7  │ 2.296e-5   │ 3.109e-5          │ 2.7597e-5   │
+│ 9   │ 8  │ 3.4512e-5  │ 4.0778e-5         │ 3.7466e-5   │
+│ 10  │ 9  │ 4.8627e-5  │ 5.7028e-5         │ 5.3766e-5   │

src/IntervalRootFinding.jl

@@ -19,7 +19,7 @@ export
     Root, is_unique,
     roots, find_roots,
     bisect, newton1d, slope,
-    slope_newton1d,
+    slope_newton1d, linear_hull,
     gauss_seidel_interval, gauss_seidel_interval!,
     gauss_seidel_contractor, gauss_seidel_contractor!,
     gauss_elimination_interval, gauss_elimination_interval!

src/linear_eq.jl

@@ -143,3 +143,33 @@ function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::Ab
 
     a \ b
 end
+
+function linear_hull(M::AbstractMatrix, r::AbstractArray)
+
+	n = size(M, 1)
+
+    ((M, r) = preconditioner(M, r))
+    M_lo = inf.(M)
+    M_hi = sup.(M)
+    # M = interval.((2*eye(n) - M_hi), M_hi)  #is it necessary?
+    if all(.≤(M_lo, zero(M_lo)))
+        return M \ r
+    end
+    P = inv(M_lo)
+    if all(.≤(eye(n), (P)))
+        H1 = P * sup.(r)
+        C = 1 ./ (2diag(P) - 1)
+        Z = ((2mid.(r)) .* diag(P)) - H1
+        H2 = C .* Z
+        for i in 1:n
+            if Z[i] < 0
+                H2[i] = Z[i]
+            end
+        end
+        H = interval.(min.(H1, H2), max.(H1, H2))
+
+        return H
+    else
+        return M \ r
+    end
+end

test/linear_eq.jl

@@ -9,7 +9,7 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding
 
     for i in 1:length(A)
         for precondition in (false, true)
-            for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval)
+            for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, linear_hull)
                 @test all(x[i] .⊆ f(A[i], b[i]))
             end
         end