Use StaticArrays and add benchmarking

src/benchmark.jl

@@ -0,0 +1,17 @@
+using IntervalArithmetic, StaticArrays, BenchmarkTools, Compat
+
+include("linear_eq.jl")
+
+function randVec(n::Int)
+           a = randn(n)
+           A = Interval.(a)
+           sA = MVector{n}(A)
+           return A, sA
+       end
+
+function randMat(n::Int)
+           a = randn(n, n)
+           A = Interval.(a)
+           sA = MMatrix{n, n}(A)
+           return A, sA
+       end

src/linear_eq.jl

@@ -1,21 +1,29 @@
+using IntervalArithmetic
 """
 Preconditions the matrix A and b with the inverse of mid(A)
 """
 function preconditioner{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}})
 
     Aᶜ = mid.(A)
-    B = pinv(Aᶜ)
+    B = inv(Aᶜ)
 
     return B*A, B*b
 
 end
 
+function gauss_seidel_interval{T}(A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100)
 
+    n = size(A, 1)
+    x = fill(-1e16..1e16, n)
+    gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter)
+    return x
+end
 """
 Iteratively solves the system of interval linear
 equations and returns the solution set. Uses the
 Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
 Keyword `precondition` to turn preconditioning off.
+Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
 """
 function gauss_seidel_interval!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}, b::Array{Interval{T}}; precondition=true, maxiter=100)
 
@@ -29,8 +37,52 @@ function gauss_seidel_interval!{T}(x::Array{Interval{T}}, A::Matrix{Interval{T}}
             for j in 1:n
                 (i == j) || (Y -= M[i, j] * x[j])
             end
-            Y = extended_div(Y, M[i, i])
-            x[i] = hull((x[i] ∩ Y[1]), x[i] ∩ Y[2])
+            Z = extended_div(Y, M[i, i])
+            x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
+        end
+    end
+    x
+end
+
+function preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}})
+
+    Aᶜ = mid.(A)
+    B = inv(Aᶜ)
+
+    return B*A, B*b
+
+end
+
+
+function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100)
+
+    n = size(A, 1)
+    x = @MVector fill(-1e16..1e16, n)
+    gauss_seidel_interval_static!(x, A, b, precondition=precondition, maxiter=maxiter)
+    return x
+end
+
+"""
+Iteratively solves the system of interval linear
+equations and returns the solution set. Uses the
+Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
+Keyword `precondition` to turn preconditioning off.
+Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
+"""
+function gauss_seidel_interval_static!{T, N}(x::MVector{N, Interval{T}}, A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100)
+
+    precondition && ((M, r) = preconditioner_static(A, b))
+
+    n = size(A, 1)
+
+    for iter in 1:maxiter
+        for i in 1:n
+            Y = r[i]
+            for j in 1:n
+                (i == j) || (Y -= M[i, j] * x[j])
+            end
+            Z = extended_div(Y, M[i, i])
+            x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
         end
     end
     x