Skip to content

Commit 7c9cbd9

Browse files
eeshan9815eeshan9815
committed
Add Gauss Elimination with preconditioning
1 parent 1388397 commit 7c9cbd9

File tree

2 files changed

+96
-141
lines changed

2 files changed

+96
-141
lines changed

‎perf/linear_eq.jl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,3 +40,24 @@ function benchmark(max=10)
4040
println(dfnew)
4141
dfnew
4242
end
43+
44+
function benchmark_elimination(max=10)
45+
df = DataFrame()
46+
df[:Method] = ["Gauss Elimination", "Base.\\"]
47+
for n in 1:max
48+
A, mA, sA = randMat(n)
49+
b, mb, sb = randVec(n)
50+
t1 = @belapsed gauss_elimination_interval($A, $b)
51+
t2 = @belapsed gauss_elimination_interval1($A, $b)
52+
df[Symbol("$n")] = [t1, t2]
53+
end
54+
a = []
55+
for i in 1:max
56+
push!(a, Symbol("$i"))
57+
end
58+
df1 = stack(df, a)
59+
dfnew = unstack(df1, :variable, :Method, :value)
60+
dfnew = rename(dfnew, :variable => :n)
61+
println(dfnew)
62+
dfnew
63+
end

‎src/linear_eq.jl

Lines changed: 75 additions & 141 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
using IntervalArithmetic
12
"""
23
Preconditions the matrix A and b with the inverse of mid(A)
34
"""
@@ -15,10 +16,6 @@ function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition
1516
n = size(A, 1)
1617
x = fill(-1e16..1e16, n)
1718

18-
if (typeof(b) <: SArray || typeof(b) <: MArray)
19-
x = MVector{n}(x)
20-
end
21-
2219
gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter)
2320
return x
2421
end
@@ -47,148 +44,11 @@ function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::Abstract
4744
end
4845
x
4946
end
50-
#
51-
# function preconditioner_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}})
52-
#
53-
# Aᶜ = mid.(A)
54-
# B = inv(Aᶜ)
55-
#
56-
# return B*A, B*b
57-
#
58-
# end
59-
#
60-
#
61-
# function gauss_seidel_interval_static{T, N}(A::MMatrix{N, N, Interval{T}}, b::MVector{N, Interval{T}}; precondition=true, maxiter=100)
62-
#
63-
# n = size(A, 1)
64-
# x = @MVector fill(-1e16..1e16, n)
65-
# gauss_seidel_interval_static!(x, A, b, precondition=precondition, maxiter=maxiter)
66-
# return x
67-
# end
68-
#
69-
# """
70-
# Iteratively solves the system of interval linear
71-
# equations and returns the solution set. Uses the
72-
# Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
73-
# Keyword `precondition` to turn preconditioning off.
74-
# Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
75-
# """
76-
# 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)
77-
#
78-
# precondition && ((A, b) = preconditioner_static(A, b))
79-
#
80-
# n = size(A, 1)
81-
#
82-
# for iter in 1:maxiter
83-
# for i in 1:n
84-
# Y = b[i]
85-
# for j in 1:n
86-
# (i == j) || (Y -= A[i, j] * x[j])
87-
# end
88-
# Z = extended_div(Y, A[i, i])
89-
# x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
90-
# end
91-
# end
92-
# x
93-
# end
94-
#
95-
# function preconditioner_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}})
96-
#
97-
# Aᶜ = mid.(A)
98-
# B = inv(Aᶜ)
99-
#
100-
# return B*A, B*b
101-
#
102-
# end
103-
#
104-
#
105-
# function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
106-
#
107-
# n = size(A, 1)
108-
# x = @MVector fill(-1e16..1e16, n)
109-
# gauss_seidel_interval_static1!(x, A, b, precondition=precondition, maxiter=maxiter)
110-
# return x
111-
# end
112-
#
113-
# """
114-
# Iteratively solves the system of interval linear
115-
# equations and returns the solution set. Uses the
116-
# Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
117-
# Keyword `precondition` to turn preconditioning off.
118-
# Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
119-
# """
120-
# function gauss_seidel_interval_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
121-
#
122-
# precondition && ((A, b) = preconditioner_static1(A, b))
123-
#
124-
# n = size(A, 1)
125-
#
126-
# for iter in 1:maxiter
127-
# for i in 1:n
128-
# Y = b[i]
129-
# for j in 1:n
130-
# (i == j) || (Y -= A[i, j] * x[j])
131-
# end
132-
# Z = extended_div(Y, A[i, i])
133-
# x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
134-
# end
135-
# end
136-
# x
137-
# end
138-
#
139-
# function preconditioner_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}})
140-
#
141-
# Aᶜ = mid.(A)
142-
# B = inv(Aᶜ)
143-
#
144-
# return B*A, B*b
145-
#
146-
# end
147-
#
148-
#
149-
# function gauss_seidel_interval_static2{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
150-
#
151-
# n = size(A, 1)
152-
# x = @SVector fill(-1e16..1e16, n)
153-
# x = gauss_seidel_interval_static2!(x, A, b, precondition=precondition, maxiter=maxiter)
154-
# return x
155-
# end
156-
#
157-
# """
158-
# Iteratively solves the system of interval linear
159-
# equations and returns the solution set. Uses the
160-
# Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
161-
# Keyword `precondition` to turn preconditioning off.
162-
# Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
163-
# """
164-
# function gauss_seidel_interval_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
165-
#
166-
# precondition && ((A, b) = preconditioner_static2(A, b))
167-
#
168-
# n = size(A, 1)
169-
#
170-
# for iter in 1:maxiter
171-
# for i in 1:n
172-
# Y = b[i]
173-
# for j in 1:n
174-
# (i == j) || (Y -= A[i, j] * x[j])
175-
# end
176-
# Z = extended_div(Y, A[i, i])
177-
# x = setindex(x, hull((x[i] ∩ Z[1]), x[i] ∩ Z[2]), i)
178-
# end
179-
# end
180-
# x
181-
# end
18247

18348
function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
18449

18550
n = size(A, 1)
18651
x = fill(-1e16..1e16, n)
187-
188-
if (typeof(b) <: SArray || typeof(b) <: MArray)
189-
x = MVector{n}(x)
190-
end
191-
19252
x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter)
19353
return x
19454
end
@@ -215,3 +75,77 @@ function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::Abstra
21575
end
21676
x
21777
end
78+
79+
function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true)
80+
81+
n = size(A, 1)
82+
x = fill(-1e16..1e16, n)
83+
84+
x = gauss_seidel_interval!(x, A, b, precondition=precondition)
85+
86+
return x
87+
end
88+
"""
89+
Iteratively solves the system of interval linear
90+
equations and returns the solution set. Uses the
91+
Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
92+
Keyword `precondition` to turn preconditioning off.
93+
Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
94+
"""
95+
function gauss_elimination_interval!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
96+
97+
precondition && ((a, b) = preconditioner(a, b))
98+
99+
n = size(a, 1)
100+
101+
p = zeros(x)
102+
103+
for i in 1:(n-1)
104+
if 0 ∈ a[i, i] # diagonal matrix is not invertible
105+
p .= entireinterval(b[1])
106+
return p .∩ x
107+
end
108+
109+
for j in (i+1):n
110+
α = a[j, i] / a[i, i]
111+
b[j] -= α * b[i]
112+
113+
for k in (i+1):n
114+
a[j, k] -= α * a[i, k]
115+
end
116+
end
117+
end
118+
119+
for i in n:-1:1
120+
sum = 0
121+
for j in (i+1):n
122+
sum += a[i, j] * p[j]
123+
end
124+
p[i] = (b[i] - sum) / a[i, i]
125+
end
126+
127+
p .∩ x
128+
end
129+
130+
function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true)
131+
132+
n = size(A, 1)
133+
x = fill(-1e16..1e16, n)
134+
135+
x = gauss_seidel_interval!(x, A, b, precondition=precondition)
136+
137+
return x
138+
end
139+
"""
140+
Iteratively solves the system of interval linear
141+
equations and returns the solution set. Uses the
142+
Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
143+
Keyword `precondition` to turn preconditioning off.
144+
Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
145+
"""
146+
function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
147+
148+
precondition && ((a, b) = preconditioner(a, b))
149+
150+
a \ b
151+
end

0 commit comments

Comments
 (0)