Add symmetric system constructor

Author: janbruedigam

src/system/system.jl

@@ -1,4 +1,8 @@
-struct System{N}
+abstract type Symmetric end                     # symmetric = true
+abstract type BlockSymmetric <: Symmetric end   # symmetric = true
+abstract type SparsitySymmetric end             # symmetric = false, Any matrix is sparsity symmetric since zero entries can be assumed nonzero
+
+struct System{N,T}
     matrix_entries::SparseMatrixCSC{GraphBasedSystems.Entry, Int64}
     vector_entries::Vector{Entry}
     diagonal_inverses::Vector{Entry}
@@ -9,23 +13,21 @@ struct System{N}
     graph::SimpleGraph{Int64}
     dfs_graph::SimpleDiGraph{Int64}
 
-    function System{T}(A, dims; force_static = false) where T
+    function System{T}(A, dims; force_static = false, symmetric = false) where T
         N = length(dims)
+        if symmetric
+            all(dims .== 1) ? S = Symmetric : S = BlockSymmetric
+        else
+            S = SparsitySymmetric
+        end
         static = force_static || all(dims.<=10)
 
         full_graph = Graph(A)
 
         matrix_entries = spzeros(Entry,N,N)
 
         for (i,dimi) in enumerate(dims)
-            for (j,dimj) in enumerate(dims)
-                if i == j
-                    matrix_entries[i,j] = Entry{T}(dimi, dimj, static = static)
-                elseif j ∈ all_neighbors(full_graph,i)
-                    matrix_entries[i,j] = Entry{T}(dimi, dimj, static = static)
-                    matrix_entries[j,i] = Entry{T}(dimj, dimi, static = static)
-                end
-            end
+                matrix_entries[i,i] = Entry{T}(dimi, dimi, static = static)
         end
 
         vector_entries = [Entry{T}(dim, static = static) for dim in dims]
@@ -65,23 +67,36 @@ struct System{N}
                 acyclic_children[v] = setdiff(acyclic_children[v], cyclic_children)
                 for c in cyclic_children
                     matrix_entries[v,c] = Entry{T}(dims[v], dims[c], static = static)
-                    matrix_entries[c,v] = Entry{T}(dims[c], dims[v], static = static)
+                    !symmetric && (matrix_entries[c,v] = Entry{T}(dims[c], dims[v], static = static))
 
                     v ∉ parents[c] && push!(parents[c],v)
                 end
-            end            
+            end  
+            for v in sub_dfs_list
+                for c in acyclic_children[v]
+                    matrix_entries[v,c] = Entry{T}(dims[v], dims[c], static = static)
+                    !symmetric && (matrix_entries[c,v] = Entry{T}(dims[c], dims[v], static = static))
+                end
+            end
         end
 
         full_dfs_graph = SimpleDiGraph(edgelist)
         cyclic_children = [unique(vcat(cycles[i]...)) for i=1:N]
 
-        new{N}(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)
     end
 
-    System(A, dims; force_static = false) = System{Float64}(A, dims; force_static = force_static)
+    System(A, dims; force_static = false, symmetric = false) = System{Float64}(A, dims; force_static, symmetric)
 end
 
-function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N}) where {N}
+function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, system::System{N,S}) where {N,S}
+    if S <: BlockSymmetric
+        print(io, "BlockSymmetric ")
+    elseif S <: Symmetric
+        print(io, "Symmetric ")
+    elseif S <: SparsitySymmetric
+        print(io, "SparsitySymmetric ")
+    end
     println(io, "System with "*string(N)*" nodes.")
     SparseArrays._show_with_braille_patterns(io, system.matrix_entries)
 end
@@ -109,4 +124,24 @@ function full_matrix(system::System{N}) where N
     return A
 end
 
+# function full_matrix(system::System{N,<:Symmetric}) where N
+#     dims = [length(system.vector_entries[i].value) for i=1:N]
+
+#     range = [1:dims[1]]
+
+#     for (i,dim) in enumerate(collect(Iterators.rest(dims, 2)))
+#         push!(range,sum(dims[1:i])+1:sum(dims[1:i])+dim)
+#     end
+#     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[row],range[col]] = system.matrix_entries[row,col].value
+#         if col != row
+#             A[range[col],range[row]] = system.matrix_entries[row,col].value'
+#         end
+#     end
+#     return A
+# end
+
 full_vector(system) = vcat(getfield.(system.vector_entries,:value)...)