Merge pull request #2 from janbruedigam/custom_dict

Speed up

Author: janbruedigam

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.1"
 
 [deps]
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]

src/GraphBasedSystems.jl

@@ -1,5 +1,6 @@
 module GraphBasedSystems
 
+using SparseArrays
 using StaticArrays
 using LightGraphs
 
@@ -14,15 +15,15 @@ export System,
     ldu_solve!
 
 
-include("custom_static.jl")
+include(joinpath("util", "custom_static.jl"))
 
-include("entry.jl")
-include("system.jl")
-include("setup_functions.jl")
+include(joinpath("system", "entry.jl"))
+include(joinpath("system", "system.jl"))
+include(joinpath("system", "setup_functions.jl"))
 
-# include("lu.jl")
-# include("llt.jl")
-# include("ldlt.jl")
-include("ldu.jl")
+include(joinpath("solvers", "lu.jl"))
+include(joinpath("solvers", "llt.jl"))
+include(joinpath("solvers", "ldlt.jl"))
+include(joinpath("solvers", "ldu.jl"))
 
 end

src/entry.jl

@@ -0,0 +1,95 @@
+# LDU factorization for unsymmetric systems
+
+function ldu_factorization_acyclic!(diagonal_v, offdiagonal_l, diagonal_c, offdiagonal_u, diagonal_inverse_c)
+    if diagonal_inverse_c.isinverted
+        invdiagonal_c = diagonal_inverse_c.value
+    else
+        invdiagonal_c = inv(diagonal_c.value)
+        diagonal_inverse_c.value = invdiagonal_c
+        diagonal_inverse_c.isinverted = true
+    end
+    offdiagonal_l.value = offdiagonal_l.value * invdiagonal_c
+    offdiagonal_u.value = invdiagonal_c * offdiagonal_u.value
+    diagonal_v.value -= offdiagonal_l.value*diagonal_c.value*offdiagonal_u.value
+    return
+end
+function ldu_factorization_cyclic!(entry_lu, offdiagonal_lu, diagonal_c, offdiagonal_ul)
+    entry_lu.value -= offdiagonal_lu.value*diagonal_c.value*offdiagonal_ul.value
+    return
+end
+
+function ldu_factorization!(system)
+    matrix_entries = system.matrix_entries
+    diagonal_inverses = system.diagonal_inverses
+    acyclic_children = system.acyclic_children
+    cycles = system.cycles
+
+    for v in system.dfs_list
+        for cyclic_children in cycles[v]
+            for c in cyclic_children
+                for cc in cyclic_children
+                    cc == c && break 
+                    cc ∉ acyclic_children[c] && continue
+                    ldu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[cc,c])
+                    ldu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,cc], matrix_entries[cc,v])
+                end
+                ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
+            end
+        end
+        for c in acyclic_children[v]
+            ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
+        end
+    end
+    return
+end
+
+function ldu_backsubstitution_l!(vector_v, offdiagonal, vector_c)
+    vector_v.value -= offdiagonal.value*vector_c.value
+    return
+end
+function ldu_backsubstitution_u!(vector_v, offdiagonal, vector_p)
+    vector_v.value -= offdiagonal.value*vector_p.value
+    return
+end
+function ldu_backsubstitution_d!(vector, diagonal, diagonal_inverse)
+    if diagonal_inverse.isinverted
+        vector.value = diagonal_inverse.value*vector.value
+    else
+        vector.value = diagonal.value\vector.value
+    end
+    diagonal_inverse.isinverted = false
+    return
+end
+
+function ldu_backsubstitution!(system)
+    matrix_entries = system.matrix_entries
+    diagonal_inverses = system.diagonal_inverses
+    vector_entries = system.vector_entries
+    acyclic_children = system.acyclic_children
+    cycles = system.cycles
+    parents = system.parents
+    dfs_list = system.dfs_list
+
+    for v in dfs_list
+        for cyclic_children in cycles[v]
+            for c in cyclic_children
+                ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
+            end
+        end
+        for c in acyclic_children[v]
+            ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
+        end
+    end
+    for v in reverse(dfs_list)
+        ldu_backsubstitution_d!(vector_entries[v], matrix_entries[v,v], diagonal_inverses[v])
+        for p in parents[v]
+            ldu_backsubstitution_u!(vector_entries[v], matrix_entries[v,p], vector_entries[p])
+        end
+    end
+end
+
+function ldu_solve!(system)
+    ldu_factorization!(system)
+    ldu_backsubstitution!(system)
+    return
+end