Comply with Dojo

Author: janbruedigam

Project.toml

@@ -5,6 +5,7 @@ version = "0.2.0"
 
 [deps]
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 

src/GraphBasedSystems.jl

@@ -1,23 +1,35 @@
 module GraphBasedSystems
 
+using LinearAlgebra
 using SparseArrays
 using StaticArrays
 using Graphs
 
 
 export System,
+    Entry,
+
     full_matrix,
     full_vector,
     randomize!,
+    reset_inverse_diagonals!,
 
     children,
     connections,
     parents,
 
     ldu_solve!,
+    ldu_factorization!,
+    ldu_backsubstitution!,
     lu_solve!,
+    lu_factorization!,
+    lu_backsubstitution!,
     ldlt_solve!,
-    llt_solve!
+    ldlt_factorization!,
+    ldlt_backsubstitution!,
+    llt_solve!,
+    llt_factorization!,
+    llt_backsubstitution!
 
 
 include(joinpath("util", "custom_static.jl"))
@@ -26,6 +38,9 @@ include(joinpath("system", "entry.jl"))
 include(joinpath("system", "system.jl"))
 include(joinpath("system", "setup_functions.jl"))
 
+include(joinpath("system", "interface.jl"))
+include(joinpath("system", "dense.jl"))
+
 include(joinpath("solvers", "lu.jl"))
 include(joinpath("solvers", "llt.jl"))
 include(joinpath("solvers", "ldlt.jl"))

src/solvers/ldlt.jl

@@ -24,6 +24,8 @@ function ldlt_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             ldlt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
@@ -83,7 +85,6 @@ function ldlt_backsubstitution!(system::System)
 end
 
 function ldlt_solve!(system::System)
-    reset_inverse_diagonals!(system)
     ldlt_factorization!(system)
     ldlt_backsubstitution!(system)
     return

src/solvers/ldu.jl

@@ -24,6 +24,8 @@ function ldu_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         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])
@@ -84,7 +86,6 @@ function ldu_backsubstitution!(system::System)
 end
 
 function ldu_solve!(system::System)
-    reset_inverse_diagonals!(system)
     ldu_factorization!(system)
     ldu_backsubstitution!(system)
     return

src/solvers/llt.jl

@@ -28,6 +28,8 @@ function llt_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             llt_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], diagonal_inverses[c])
@@ -92,7 +94,6 @@ function llt_backsubstitution!(system::System)
 end
 
 function llt_solve!(system::System)
-    reset_inverse_diagonals!(system)
     llt_factorization!(system)
     llt_backsubstitution!(system)
     return

src/solvers/lu.jl

@@ -23,6 +23,8 @@ function lu_factorization!(system::System)
     acyclic_children = system.acyclic_children
     cyclic_children = system.cyclic_children
 
+    reset_inverse_diagonals!(system)
+
     for v in system.dfs_list
         for c in acyclic_children[v]
             lu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
@@ -83,7 +85,6 @@ function lu_backsubstitution!(system::System)
 end
 
 function lu_solve!(system::System)
-    reset_inverse_diagonals!(system)
     lu_factorization!(system)
     lu_backsubstitution!(system)
     return

src/system/dense.jl

@@ -0,0 +1,19 @@
+full_matrix(system::System) = full_matrix(system.matrix_entries, issymmetric(system), system.dims, system.dims)
+# There probably exists a smarter way of getting the dense matrix from the spares one 
+function full_matrix(matrix_entries::SparseMatrixCSC, symmetric::Bool, dimensions_rows, dimensions_cols)
+    range_dict_rows = ranges(dimensions_rows)
+    range_dict_cols = ranges(dimensions_cols)
+    A = zeros(sum(dimensions_rows), sum(dimensions_cols))
+
+    for (i,row) in enumerate(matrix_entries.rowval)
+        col = findfirst(x -> i<x, matrix_entries.colptr)-1
+        A[range_dict_rows[row],range_dict_cols[col]] = matrix_entries[row,col].value
+        if symmetric && col != row
+            A[range_dict_cols[col],range_dict_rows[row]] = matrix_entries[row,col].value'
+        end
+    end
+    return A
+end
+
+full_vector(system::System) = full_vector(system.vector_entries, system.dims)
+full_vector(vector_entries::AbstractVector, dimensions_rows) = vcat([getfield(vector_entries[i],:value) for i=1:size(dimensions_rows)[1]]...)

src/system/entry.jl

@@ -25,7 +25,7 @@ function Base.zero(::Type{Entry{ET}}) where ET
     return Entry{ET.parameters[2]}(dims...)
 end
 function Base.zero(::Type{Entry})
-    return 0
+    return nothing
 end
 
 function randomize!(entry::Entry, rand_function = randn)

src/system/interface.jl

@@ -0,0 +1,54 @@
+@inline children(system::System, v) = outneighbors(system.dfs_graph, v) # all direct children of v
+@inline connections(system::System, v) = neighbors(system.graph, v)     # all connected nodes of v
+@inline parents(system::System, v) = inneighbors(system.dfs_graph, v)   # same elements as system.parents[v], but potentially different order
+
+LinearAlgebra.issymmetric(::System{N, Symmetric}) where N = true
+LinearAlgebra.issymmetric(::System{N, Unsymmetric}) where N = false
+
+function ranges(system::System{N}) where N
+    dims = system.dims
+    range_dict = Dict(1=>1:dims[1])
+    for i=2:N
+        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
+    end
+
+    return range_dict
+end
+function ranges(dims::AbstractVector)
+    range_dict = Dict(1=>1:dims[1])
+    for i=2:size(dims)[1]
+        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
+    end
+
+    return range_dict
+end
+
+function randomize!(system::System, rand_function = randn)
+    for entry in system.matrix_entries.nzval
+        randomize!(entry, rand_function)
+    end
+    for entry in system.vector_entries
+        randomize!(entry, rand_function)
+    end
+end
+
+function randomize!(system::System{N,<:Symmetric}, rand_function = randn) where N
+    matrix_entries = system.matrix_entries
+
+    for entry in matrix_entries.nzval
+        randomize!(entry, rand_function)
+    end
+    for i=1:N
+        matrix_entries[i,i].value += matrix_entries[i,i].value'
+    end
+
+    for entry in system.vector_entries
+        randomize!(entry, rand_function)
+    end
+end
+
+function reset_inverse_diagonals!(system::System)
+    for entry in system.diagonal_inverses
+        entry.isinverted = false
+    end
+end

src/system/system.jl

@@ -11,6 +11,7 @@ struct System{N,T}
     dfs_list::SVector{N,Int64}      # depth-first search list of nodes [last-found node, ..., first-found node]
     graph::SimpleGraph{Int64}       # the graph built from the adjacency matrix
     dfs_graph::SimpleDiGraph{Int64} # the directed graph built from the depth-first search
+    dims::Vector{Int64} # Dimensions of the matrix entries
 
     function System{T}(A, dims; force_static = false, symmetric = false) where T
         N = length(dims)
@@ -79,7 +80,7 @@ struct System{N,T}
         cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]
         cyclic_children = [intersect(dfs_list, cyclic_children[i]) for i=1:N]
 
-        new{N,S}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph)
+        new{N,S}(matrix_entries, vector_entries, diagonal_inverses, acyclic_children, cyclic_children, parents, dfs_list, full_graph, full_dfs_graph, dims)
     end
 
     System(A, dims; force_static = false, symmetric = false) = System{Float64}(A, dims; force_static, symmetric)
@@ -94,79 +95,3 @@ function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N,S}
     println(io, "System with "*string(N)*" nodes.")
     SparseArrays._show_with_braille_patterns(io, system.matrix_entries)
 end
-
-
-@inline children(system::System, v) = outneighbors(system.dfs_graph, v) # all direct children of v
-@inline connections(system::System, v) = neighbors(system.graph, v)     # all connected nodes of v
-@inline parents(system::System, v) = inneighbors(system.dfs_graph, v)   # same elements as system.parents[v], but potentially different order
-
-dimensions(system::System{N}) where N = [size(system.vector_entries[i].value)[1] for i=1:N]
-function ranges(system::System{N}) where N
-    dims = dimensions(system)
-    range_dict = Dict(1=>1:dims[1])
-    for i=2:N
-        range_dict[i] = last(range_dict[i-1])+1:sum(dims[1:i])
-    end
-
-    return range_dict
-end
-
-# There probably exists a smarter way of getting the dense matrix from the spares one 
-function full_matrix(system::System{N}) where N
-    dims = dimensions(system)
-    range_dict = ranges(system)
-    A = zeros(sum(dims), sum(dims))
-
-    for (i,row) in enumerate(system.matrix_entries.rowval)
-        col = findfirst(x -> i<x, system.matrix_entries.colptr)-1
-        A[range_dict[row],range_dict[col]] = system.matrix_entries[row,col].value
-    end
-    return A
-end
-
-function full_matrix(system::System{N,<:Symmetric}) where N
-    dims = dimensions(system)
-    range_dict = ranges(system)
-    A = zeros(sum(dims),sum(dims))
-
-    for (i,row) in enumerate(system.matrix_entries.rowval)
-        col = findfirst(x -> i<x, system.matrix_entries.colptr)-1
-        A[range_dict[row],range_dict[col]] = system.matrix_entries[row,col].value
-        if col != row
-            A[range_dict[col],range_dict[row]] = system.matrix_entries[row,col].value'
-        end
-    end
-    return A
-end
-
-full_vector(system::System) = vcat(getfield.(system.vector_entries,:value)...)
-
-function randomize!(system::System, rand_function = randn)
-    for entry in system.matrix_entries.nzval
-        randomize!(entry, rand_function)
-    end
-    for entry in system.vector_entries
-        randomize!(entry, rand_function)
-    end
-end
-
-function randomize!(system::System{N,<:Symmetric}, rand_function = randn) where N
-    matrix_entries = system.matrix_entries
-
-    for entry in matrix_entries.nzval
-        randomize!(entry, rand_function)
-    end
-    for i=1:N
-        matrix_entries[i,i].value += matrix_entries[i,i].value'
-    end
-
-    for entry in system.vector_entries
-        randomize!(entry, rand_function)
-    end
-end
-
-function reset_inverse_diagonals!(system::System)
-    for entry in system.diagonal_inverses
-        entry.isinverted = false
-    end
-end