Skip to content

Commit 4fbf8db

Browse files
janbruedigamjanbruedigam
authored
Merge pull request #5 from janbruedigam/updates
Clean up
2 parents 8df3a63 + b5c7771 commit 4fbf8db

File tree

18 files changed

+470
-258
lines changed

18 files changed

+470
-258
lines changed

‎.gitignore‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,3 @@
11
/Manifest.toml
2-
/testscripts
2+
/testscripts
3+
/benchmark/tune.json

‎README.md‎

Lines changed: 4 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,9 @@
22
[![Build Status](https://github.com/janbruedigam/GraphBasedSystems.jl/workflows/CI/badge.svg)](https://github.com/janbruedigam/GraphBasedSystems.jl/actions?query=workflow%3ACI)
33
[![codecov](https://codecov.io/gh/janbruedigam/GraphBasedSystems.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/janbruedigam/GraphBasedSystems.jl)
44

5-
A package for efficiently solving linear systems of equations based on graph properties. Any matrix representing an underlying real-world system can be represented by a graph, for example a mechanical system, checmical molecules, neural networks, ... The code of the package is currently used for the robotics simulator [Dojo.jl](https://github.com/dojo-sim/Dojo.jl)
5+
A package for efficiently solving linear systems of equations based on graph properties. Any matrix representing an underlying real-world system can be represented by a graph, for example a mechanical system, chemical molecules, neural networks, ... The code of the package is currently used for the robotics simulator [Dojo.jl](https://github.com/dojo-sim/Dojo.jl)
66

7-
By providing the adjacency matrix of a graph-based system, the `GraphBasedSystems` package can automatically exploit the existing sparsity when solving the linear system to speed up calculations. Currently, the LU and LDU decomposition/backsubstitution are implemented.
7+
By providing the adjacency matrix of a graph-based system, the `GraphBasedSystems` package can automatically exploit the existing sparsity when solving the linear system to speed up calculations. Currently, the LDU, LU, LDLt, and LLt decomposition/backsubstitution are implemented. LLt is currently slow due to memory allocations.
88

99
```julia
1010
using GraphBasedSystems
@@ -18,16 +18,10 @@ graph_matrix = [ # The adjacency matrix for the underlying graph
1818

1919
dimensions = [2; 3; 0; 1] # The dimension of row/column
2020

21-
system = System{Float64}(graph_matrix, dimensions) # The resulting linear system
21+
system = System{Float64}(graph_matrix, dimensions; symmetric=false) # The resulting linear system. Set symmetric=true for symmetric systems
2222

23-
for entry in system.matrix_entries.nzval # Randomize all matrix entries
24-
GraphBasedSystems.randomize!(entry)
25-
end
23+
randomize!(system) # Randomize all system entries
2624
system.matrix_entries[1,2].value = rand(2,3) # Directly set the value of a matrix entry
27-
28-
for entry in system.vector_entries # Randomize all vector entries
29-
GraphBasedSystems.randomize!(entry)
30-
end
3125
system.vector_entries[4].value = rand(1) # Directly set the value of a vector entry
3226

3327
A = full_matrix(system) # Inspect the matrix

‎benchmark/example_benchmark.jl‎

Lines changed: 28 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,6 @@
11
using GraphBasedSystems
2+
using LinearAlgebra
3+
24
A = [
35
0 1 0 1 1 0 1 0 1 0
46
1 0 1 0 0 0 0 0 0 0
@@ -75,17 +77,32 @@ A[13,21] = A[21,13] = 1
7577
A[16,30] = A[30,16] = 1
7678
A[20,36] = A[36,20] = 1
7779

78-
79-
system = System{Float64}(A, rand(0:3,size(A)[1]))
80-
81-
function init!(system)
82-
for entry in system.matrix_entries.nzval
83-
GraphBasedSystems.randomize!(entry)
84-
end
85-
for entry in system.vector_entries
86-
GraphBasedSystems.randomize!(entry)
80+
function randomize_posdef!(system)
81+
randomize!(system,rand)
82+
for i=1:size(A)[1]
83+
system.matrix_entries[i,i].value += 1000*I
8784
end
8885
end
8986

90-
SUITE["ldu"] = @benchmarkable ldu_solve!($system) samples=2 setup=(init!($system))
91-
SUITE["lu"] = @benchmarkable lu_solve!($system) samples=2 setup=(init!($system))
87+
88+
system = System{Float64}(A, ones(Int,size(A)[1])*3)
89+
systemldlt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
90+
systemllt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
91+
92+
SUITE["ldu"] = @benchmarkable ldu_solve!($system) samples=2 setup=(randomize!($system))
93+
SUITE["lu"] = @benchmarkable lu_solve!($system) samples=2 setup=(randomize!($system))
94+
SUITE["ldlt"] = @benchmarkable ldlt_solve!($systemldlt) samples=2 setup=(randomize!($systemldlt))
95+
SUITE["llt"] = @benchmarkable llt_solve!($systemllt) samples=2 setup=(randomize_posdef!($systemllt))
96+
97+
# A = [
98+
# 0 1 1 1 1
99+
# 1 0 1 1 1
100+
# 1 1 0 1 1
101+
# 1 1 1 0 1
102+
# 1 1 1 1 0
103+
# ]
104+
105+
106+
107+
108+

‎src/GraphBasedSystems.jl‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,14 +8,16 @@ using Graphs
88
export System,
99
full_matrix,
1010
full_vector,
11-
reordered_matrix,
12-
reordered_vector,
11+
randomize!,
12+
1313
children,
1414
connections,
1515
parents,
1616

1717
ldu_solve!,
18-
lu_solve!
18+
lu_solve!,
19+
ldlt_solve!,
20+
llt_solve!
1921

2022

2123
include(joinpath("util", "custom_static.jl"))

‎src/solvers/ldlt.jl‎

Lines changed: 90 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,90 @@
1-
# LDLᵀ factorization for block-symmetric systems
1+
# LDLᵀ (Cholesky) factorization for symmetric systems
2+
3+
function ldlt_factorization_acyclic!(diagonal_v, offdiagonal_l, diagonal_c, diagonal_inverse_c)
4+
if diagonal_inverse_c.isinverted
5+
invdiagonal_c = diagonal_inverse_c.value
6+
else
7+
invdiagonal_c = inv(diagonal_c.value)
8+
diagonal_inverse_c.value = invdiagonal_c
9+
diagonal_inverse_c.isinverted = true
10+
end
11+
offdiagonal_l.value = offdiagonal_l.value * invdiagonal_c
12+
diagonal_v.value -= offdiagonal_l.value * diagonal_c.value * offdiagonal_l.value'
13+
return
14+
end
15+
function ldlt_factorization_cyclic!(entry_lu, offdiagonal_lu, diagonal_c, offdiagonal_ul)
16+
entry_lu.value -= offdiagonal_lu.value * diagonal_c.value * offdiagonal_ul.value'
17+
return
18+
end
19+
20+
21+
function ldlt_factorization!(system::System)
22+
matrix_entries = system.matrix_entries
23+
diagonal_inverses = system.diagonal_inverses
24+
acyclic_children = system.acyclic_children
25+
cyclic_children = system.cyclic_children
26+
27+
for v in system.dfs_list
28+
for c in acyclic_children[v]
29+
ldlt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
30+
end
31+
for c in cyclic_children[v]
32+
for cc in cyclic_children[v]
33+
cc == c && break
34+
(cc ∉ acyclic_children[c] && cc ∉ cyclic_children[c]) && continue
35+
ldlt_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[c,cc])
36+
end
37+
ldlt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
38+
end
39+
end
40+
return
41+
end
42+
43+
function ldlt_backsubstitution_l!(vector_v, offdiagonal, vector_c)
44+
vector_v.value -= offdiagonal.value * vector_c.value
45+
return
46+
end
47+
function ldlt_backsubstitution_lt!(vector_v, offdiagonal, vector_p)
48+
vector_v.value -= offdiagonal.value' * vector_p.value
49+
return
50+
end
51+
function ldlt_backsubstitution_d!(vector, diagonal, diagonal_inverse)
52+
if diagonal_inverse.isinverted
53+
vector.value = diagonal_inverse.value * vector.value
54+
else
55+
vector.value = diagonal.value \ vector.value
56+
end
57+
return
58+
end
59+
60+
function ldlt_backsubstitution!(system::System)
61+
matrix_entries = system.matrix_entries
62+
diagonal_inverses = system.diagonal_inverses
63+
vector_entries = system.vector_entries
64+
acyclic_children = system.acyclic_children
65+
cyclic_children = system.cyclic_children
66+
parents = system.parents
67+
dfs_list = system.dfs_list
68+
69+
for v in dfs_list
70+
for c in cyclic_children[v]
71+
ldlt_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
72+
end
73+
for c in acyclic_children[v]
74+
ldlt_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
75+
end
76+
end
77+
for v in reverse(dfs_list)
78+
ldlt_backsubstitution_d!(vector_entries[v], matrix_entries[v,v], diagonal_inverses[v])
79+
for p in parents[v]
80+
ldlt_backsubstitution_lt!(vector_entries[v], matrix_entries[p,v], vector_entries[p])
81+
end
82+
end
83+
end
84+
85+
function ldlt_solve!(system::System)
86+
reset_inverse_diagonals!(system)
87+
ldlt_factorization!(system)
88+
ldlt_backsubstitution!(system)
89+
return
90+
end

‎src/solvers/ldu.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ function ldu_factorization_cyclic!(entry_lu, offdiagonal_lu, diagonal_c, offdiag
1818
return
1919
end
2020

21-
function ldu_factorization!(system)
21+
function ldu_factorization!(system::System)
2222
matrix_entries = system.matrix_entries
2323
diagonal_inverses = system.diagonal_inverses
2424
acyclic_children = system.acyclic_children
@@ -31,7 +31,7 @@ function ldu_factorization!(system)
3131
for c in cyclic_children[v]
3232
for cc in cyclic_children[v]
3333
cc == c && break
34-
(cc ∉ children(system,c) && cc ∉ cyclic_children[c]) && continue
34+
(cc ∉ acyclic_children[c] && cc ∉ cyclic_children[c]) && continue
3535
ldu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[cc,c])
3636
ldu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,cc], matrix_entries[cc,v])
3737
end
@@ -55,11 +55,10 @@ function ldu_backsubstitution_d!(vector, diagonal, diagonal_inverse)
5555
else
5656
vector.value = diagonal.value \ vector.value
5757
end
58-
diagonal_inverse.isinverted = false
5958
return
6059
end
6160

62-
function ldu_backsubstitution!(system)
61+
function ldu_backsubstitution!(system::System)
6362
matrix_entries = system.matrix_entries
6463
diagonal_inverses = system.diagonal_inverses
6564
vector_entries = system.vector_entries
@@ -84,7 +83,8 @@ function ldu_backsubstitution!(system)
8483
end
8584
end
8685

87-
function ldu_solve!(system)
86+
function ldu_solve!(system::System)
87+
reset_inverse_diagonals!(system)
8888
ldu_factorization!(system)
8989
ldu_backsubstitution!(system)
9090
return

‎src/solvers/llt.jl‎

Lines changed: 99 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,99 @@
1-
# LLᵀ (Choleski) factorization for symmetric systems
1+
# LLᵀ (Cholesky) factorization for symmetric systems
2+
3+
function llt_factorization_acyclic!(diagonal_v, offdiagonal_l, diagonal_c, diagonal_inverse_c)
4+
if diagonal_inverse_c.isinverted
5+
invdiagonal_c = diagonal_inverse_c.value
6+
else
7+
invdiagonal_c = inv(diagonal_c.value)
8+
diagonal_inverse_c.value = invdiagonal_c
9+
diagonal_inverse_c.isinverted = true
10+
end
11+
offdiagonal_l.value = offdiagonal_l.value * invdiagonal_c
12+
diagonal_v.value -= offdiagonal_l.value * offdiagonal_l.value'
13+
return
14+
end
15+
function llt_factorization_cyclic!(entry_lu, offdiagonal_lu, offdiagonal_ul)
16+
entry_lu.value -= offdiagonal_lu.value * offdiagonal_ul.value'
17+
return
18+
end
19+
# TODO currently allocates memory
20+
function llt_diagonal_root!(diagonal)
21+
diagonal.value = sqrt(diagonal.value)
22+
end
23+
24+
25+
function llt_factorization!(system::System)
26+
matrix_entries = system.matrix_entries
27+
diagonal_inverses = system.diagonal_inverses
28+
acyclic_children = system.acyclic_children
29+
cyclic_children = system.cyclic_children
30+
31+
for v in system.dfs_list
32+
for c in acyclic_children[v]
33+
llt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
34+
end
35+
for c in cyclic_children[v]
36+
for cc in cyclic_children[v]
37+
cc == c && break
38+
(cc ∉ acyclic_children[c] && cc ∉ cyclic_children[c]) && continue
39+
llt_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[c,cc])
40+
end
41+
llt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
42+
end
43+
llt_diagonal_root!(matrix_entries[v,v])
44+
end
45+
return
46+
end
47+
48+
function llt_backsubstitution_l!(vector_v, offdiagonal, vector_c)
49+
vector_v.value -= offdiagonal.value * vector_c.value
50+
return
51+
end
52+
function llt_backsubstitution_lt!(vector_v, offdiagonal, vector_p)
53+
vector_v.value -= offdiagonal.value' * vector_p.value
54+
return
55+
end
56+
function llt_backsubstitution_d!(vector, diagonal, diagonal_inverse)
57+
if diagonal_inverse.isinverted
58+
vector.value = diagonal_inverse.value * vector.value
59+
else
60+
invdiagonal = inv(diagonal.value)
61+
diagonal_inverse.value = invdiagonal
62+
diagonal_inverse.isinverted = true
63+
vector.value = diagonal_inverse.value * vector.value
64+
end
65+
return
66+
end
67+
68+
function llt_backsubstitution!(system::System)
69+
matrix_entries = system.matrix_entries
70+
diagonal_inverses = system.diagonal_inverses
71+
vector_entries = system.vector_entries
72+
acyclic_children = system.acyclic_children
73+
cyclic_children = system.cyclic_children
74+
parents = system.parents
75+
dfs_list = system.dfs_list
76+
77+
for v in dfs_list
78+
for c in cyclic_children[v]
79+
llt_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
80+
end
81+
for c in acyclic_children[v]
82+
llt_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
83+
end
84+
llt_backsubstitution_d!(vector_entries[v], matrix_entries[v,v], diagonal_inverses[v])
85+
end
86+
for v in reverse(dfs_list)
87+
for p in parents[v]
88+
llt_backsubstitution_lt!(vector_entries[v], matrix_entries[p,v], vector_entries[p])
89+
end
90+
llt_backsubstitution_d!(vector_entries[v], matrix_entries[v,v], diagonal_inverses[v])
91+
end
92+
end
93+
94+
function llt_solve!(system::System)
95+
reset_inverse_diagonals!(system)
96+
llt_factorization!(system)
97+
llt_backsubstitution!(system)
98+
return
99+
end

‎src/solvers/lu.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ function lu_factorization_cyclic!(entry_lu, offdiagonal_lu, offdiagonal_ul)
1717
return
1818
end
1919

20-
function lu_factorization!(system)
20+
function lu_factorization!(system::System)
2121
matrix_entries = system.matrix_entries
2222
diagonal_inverses = system.diagonal_inverses
2323
acyclic_children = system.acyclic_children
@@ -30,7 +30,7 @@ function lu_factorization!(system)
3030
for c in cyclic_children[v]
3131
for cc in cyclic_children[v]
3232
cc == c && break
33-
(cc ∉ children(system,c) && cc ∉ cyclic_children[c]) && continue
33+
(cc ∉ acyclic_children[c] && cc ∉ cyclic_children[c]) && continue
3434
lu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,c])
3535
lu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,v])
3636
end
@@ -54,11 +54,10 @@ function lu_backsubstitution_d!(vector, diagonal, diagonal_inverse)
5454
else
5555
vector.value = diagonal.value \ vector.value
5656
end
57-
diagonal_inverse.isinverted = false
5857
return
5958
end
6059

61-
function lu_backsubstitution!(system)
60+
function lu_backsubstitution!(system::System)
6261
matrix_entries = system.matrix_entries
6362
diagonal_inverses = system.diagonal_inverses
6463
vector_entries = system.vector_entries
@@ -83,7 +82,8 @@ function lu_backsubstitution!(system)
8382
end
8483
end
8584

86-
function lu_solve!(system)
85+
function lu_solve!(system::System)
86+
reset_inverse_diagonals!(system)
8787
lu_factorization!(system)
8888
lu_backsubstitution!(system)
8989
return

0 commit comments

Comments
 (0)