Skip to content

Commit 9a9cae0

Browse files
janbruedigamjanbruedigam
authored
Merge pull request #3 from janbruedigam/cycle_fixes
Cycle fixes
2 parents 5ef6c4f + c279ac1 commit 9a9cae0

File tree

4 files changed

+71
-44
lines changed

4 files changed

+71
-44
lines changed

‎src/solvers/ldu.jl‎

Lines changed: 11 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -22,19 +22,17 @@ function ldu_factorization!(system)
2222
matrix_entries = system.matrix_entries
2323
diagonal_inverses = system.diagonal_inverses
2424
acyclic_children = system.acyclic_children
25-
cycles = system.cycles
25+
cyclic_children = system.cyclic_children
2626

2727
for v in system.dfs_list
28-
for cyclic_children in cycles[v]
29-
for c in cyclic_children
30-
for cc in cyclic_children
31-
cc == c && break
32-
cc ∉ acyclic_children[c] && continue
33-
ldu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[cc,c])
34-
ldu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,cc], matrix_entries[cc,v])
35-
end
36-
ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
28+
for c in cyclic_children[v]
29+
for cc in cyclic_children[v]
30+
cc == c && break
31+
(cc ∉ children(system,c) && cc ∉ cyclic_children[c]) && continue
32+
ldu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[cc,c])
33+
ldu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,cc], matrix_entries[cc,v])
3734
end
35+
ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
3836
end
3937
for c in acyclic_children[v]
4038
ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
@@ -66,15 +64,13 @@ function ldu_backsubstitution!(system)
6664
diagonal_inverses = system.diagonal_inverses
6765
vector_entries = system.vector_entries
6866
acyclic_children = system.acyclic_children
69-
cycles = system.cycles
67+
cyclic_children = system.cyclic_children
7068
parents = system.parents
7169
dfs_list = system.dfs_list
7270

7371
for v in dfs_list
74-
for cyclic_children in cycles[v]
75-
for c in cyclic_children
76-
ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
77-
end
72+
for c in cyclic_children[v]
73+
ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
7874
end
7975
for c in acyclic_children[v]
8076
ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])

‎src/system/setup_functions.jl‎

Lines changed: 14 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -40,16 +40,6 @@ function split_adjacency(A)
4040
return graphs, roots
4141
end
4242

43-
function reindex_graph(g, ids)
44-
originaledges = collect(edges(g))
45-
newedges = LightGraphs.SimpleEdge{Int64}[]
46-
for edge in originaledges
47-
push!(newedges, LightGraphs.SimpleEdge(ids[edge.src], ids[edge.dst]))
48-
end
49-
50-
return SimpleGraphFromIterator(newedges)
51-
end
52-
5343
function cycle_parent_children(cyclic_members, parents)
5444
parent = -1
5545
for member in cyclic_members
@@ -59,4 +49,18 @@ function cycle_parent_children(cyclic_members, parents)
5949
end
6050

6151
return parent, setdiff(cyclic_members, parent)
52+
end
53+
54+
function lump_cycles!(cycles, cyclic_children)
55+
inds = Int64[]
56+
for (i,cycle) in enumerate(cycles)
57+
if cyclic_children ⊆ cycle
58+
return
59+
elseif cyclic_children ⊇ cycle
60+
append!(inds,i)
61+
end
62+
end
63+
push!(cycles, cyclic_children)
64+
deleteat!(cycles, inds)
65+
return
6266
end

‎src/system/system.jl‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -3,8 +3,8 @@ struct System{N}
33
vector_entries::Vector{Entry}
44
diagonal_inverses::Vector{Entry}
55
acyclic_children::Vector{Vector{Int64}} # Contains direct children that are not part of a cycle
6-
cycles::Vector{Vector{Vector{Int64}}} # Contains direct and indirect children that are part of a cycle (sorted in cycles)
7-
parents::Vector{Vector{Int64}} # Contains direct and cycle-opening parent (2 elemtents at most)
6+
cyclic_children::Vector{Vector{Int64}} # Contains direct and indirect children that are part of a cycle
7+
parents::Vector{Vector{Int64}} # Contains direct and cycle-opening parents
88
dfs_list::SVector{N,Int64}
99
graph::SimpleGraph{Int64}
1010
dfs_graph::SimpleDiGraph{Int64}
@@ -60,21 +60,22 @@ struct System{N}
6060
for cyclic_members in simplecycles(cycle_dfs_graph)
6161
v, cyclic_children = cycle_parent_children(cyclic_members, parents)
6262
cyclic_children = reverse(cyclic_children)
63-
push!(cycles[v], cyclic_children)
63+
lump_cycles!(cycles[v], cyclic_children)
6464

6565
acyclic_children[v] = setdiff(acyclic_children[v], cyclic_children)
6666
for c in cyclic_children
6767
matrix_entries[v,c] = Entry{T}(dims[v], dims[c], static = static)
6868
matrix_entries[c,v] = Entry{T}(dims[c], dims[v], static = static)
6969

70-
v != parents[c][1] && push!(parents[c],v)
70+
v ∉ parents[c] && push!(parents[c],v)
7171
end
7272
end
7373
end
7474

7575
full_dfs_graph = SimpleDiGraph(edgelist)
76+
cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]
7677

77-
new{N}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cycles, parents, dfs_list, full_graph, full_dfs_graph)
78+
new{N}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph)
7879
end
7980

8081
System(A, dims; force_static = false) = System{Float64}(A, dims; force_static = force_static)

‎test/ldu_test.jl‎

Lines changed: 40 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@ using GraphBasedSystems
22
using LinearAlgebra
33
GBS = GraphBasedSystems
44

5-
Z = zeros(Int64,10,10)
65
A = [
76
0 1 0 1 1 0 1 0 1 0
87
1 0 1 0 0 0 0 0 0 0
@@ -16,24 +15,51 @@ A = [
1615
0 0 0 0 0 0 0 1 0 0
1716
]
1817

18+
B = [
19+
0 1 0 0 0 1 0 0 0
20+
1 0 0 0 0 1 0 0 0
21+
0 0 0 0 0 1 1 0 1
22+
0 0 0 0 0 0 1 1 0
23+
0 0 0 0 0 0 0 1 1
24+
1 1 1 0 0 0 1 0 0
25+
0 0 1 1 0 1 0 1 1
26+
0 0 0 1 1 0 1 0 1
27+
0 0 1 0 1 0 1 1 0
28+
]
29+
30+
C = [
31+
0 1 1 0 0 0
32+
1 0 0 1 0 0
33+
1 0 0 1 1 0
34+
0 1 1 0 0 1
35+
0 0 1 0 0 1
36+
0 0 0 1 1 0
37+
]
38+
39+
ZAA = zeros(Int64,10,10)
40+
ZAB = zeros(Int64,10,9)
41+
ZBA = ZAB'
42+
ZAC = zeros(Int64,10,6)
43+
ZCA = ZAC'
44+
ZBC = zeros(Int64,9,6)
45+
ZCB = ZBC'
46+
1947
# Graph 1 is disconnected, 2-3-4-5 are connected, 6 is disconnected
2048

2149
A = [
22-
A Z Z Z Z Z
23-
Z A Z Z Z Z
24-
Z Z A Z Z Z
25-
Z Z Z A Z Z
26-
Z Z Z Z A Z
27-
Z Z Z Z Z A
50+
A ZAA ZAB ZAC ZAA ZAA
51+
ZAA A ZAB ZAC ZAA ZAA
52+
ZBA ZBA B ZBC ZBA ZBA
53+
ZCA ZCA ZCB C ZCA ZCA
54+
ZAA ZAA ZAB ZAC A ZAA
55+
ZAA ZAA ZAB ZAC ZAA A
2856
]
2957

3058
A[13,21] = A[21,13] = 1
31-
A[16,31] = A[31,16] = 1
32-
A[20,41] = A[41,20] = 1
59+
A[16,30] = A[30,16] = 1
60+
A[20,36] = A[36,20] = 1
3361

3462

35-
ids = rand(1:1000,1000)
36-
ids = unique(ids)
3763
system = System{Float64}(A, rand(0:3,size(A)[1]))
3864

3965
for i=1:10
@@ -44,9 +70,9 @@ for i=1:10
4470
GBS.randomize!(entry)
4571
end
4672

47-
B = full_matrix(system)
48-
b = full_vector(system)
73+
F = full_matrix(system)
74+
f = full_vector(system)
4975
ldu_solve!(system)
50-
@test maximum(abs.(full_vector(system)-B\b)) < 1e-3
76+
@test maximum(abs.(full_vector(system)-F\f)) < 1e-3
5177
end
5278

0 commit comments

Comments
 (0)