Merge pull request #7 from janbruedigam/add_system_algebra

Add system algebra

Author: janbruedigam

benchmark/adjacency_matrix.jl

@@ -0,0 +1,75 @@
+A = [
+    0 1 0 1 1 0 1 0 1 0
+    1 0 1 0 0 0 0 0 0 0
+    0 1 0 0 0 0 0 0 0 0
+    1 0 0 0 1 0 0 0 0 0
+    1 0 0 1 0 1 0 0 0 0
+    0 0 0 0 1 0 0 0 0 0
+    1 0 0 0 0 0 0 1 0 0
+    0 0 0 0 0 0 1 0 1 1
+    1 0 0 0 0 0 0 1 0 0
+    0 0 0 0 0 0 0 1 0 0
+]
+
+B = [
+    0 1 0 0 0 1 0 0 0
+    1 0 0 0 0 1 0 0 0
+    0 0 0 0 0 1 1 0 1
+    0 0 0 0 0 0 1 1 0
+    0 0 0 0 0 0 0 1 1
+    1 1 1 0 0 0 1 0 0
+    0 0 1 1 0 1 0 1 1
+    0 0 0 1 1 0 1 0 1
+    0 0 1 0 1 0 1 1 0
+]
+
+C = [
+    0 1 1 0 0 0
+    1 0 0 1 0 0
+    1 0 0 1 1 0
+    0 1 1 0 0 1
+    0 0 1 0 0 1
+    0 0 0 1 1 0
+]
+
+D = [
+    0 1 1 0 1 0
+    1 0 0 1 0 1
+    1 0 0 1 0 0
+    0 1 1 0 0 0
+    1 0 0 0 0 1
+    0 1 0 0 1 0]
+
+ZAA = zeros(Int64, 10, 10)
+ZAB = zeros(Int64, 10, 9)
+ZAC = zeros(Int64, 10, 6)
+ZAD = zeros(Int64, 10, 6)
+ZBB = zeros(Int64, 9, 9)
+ZBC = zeros(Int64, 9, 6)
+ZBD = zeros(Int64, 9, 6)
+ZCC = zeros(Int64, 6, 6)
+ZCD = zeros(Int64, 6, 6)
+
+ZBA = ZAB'
+ZCA = ZAC'
+ZDA = ZAD'
+ZCB = ZBC'
+ZDB = ZBD'
+ZDC = ZCD'
+
+
+# Graph 1 is disconnected, 2-3-4-5 are connected, 6 is disconnected, 7 is disconnected
+
+A = [
+    A   ZAA ZAB ZAC ZAA ZAA ZAD
+    ZAA A   ZAB ZAC ZAA ZAA ZAD
+    ZBA ZBA B   ZBC ZBA ZBA ZBD
+    ZCA ZCA ZCB C   ZCA ZCA ZCD
+    ZAA ZAA ZAB ZAC A   ZAA ZAD
+    ZAA ZAA ZAB ZAC ZAA A   ZAD
+    ZDA ZDA ZDB ZDC ZDA ZDA D
+]
+
+A[13,21] = A[21,13] = 1
+A[16,30] = A[30,16] = 1
+A[20,36] = A[36,20] = 1

benchmark/example_benchmark.jl

@@ -1,85 +1,12 @@
 using GraphBasedSystems
 using LinearAlgebra
 
-A = [
-    0 1 0 1 1 0 1 0 1 0
-    1 0 1 0 0 0 0 0 0 0
-    0 1 0 0 0 0 0 0 0 0
-    1 0 0 0 1 0 0 0 0 0
-    1 0 0 1 0 1 0 0 0 0
-    0 0 0 0 1 0 0 0 0 0
-    1 0 0 0 0 0 0 1 0 0
-    0 0 0 0 0 0 1 0 1 1
-    1 0 0 0 0 0 0 1 0 0
-    0 0 0 0 0 0 0 1 0 0
-]
+include("adjacency_matrix.jl")
 
-B = [
-    0 1 0 0 0 1 0 0 0
-    1 0 0 0 0 1 0 0 0
-    0 0 0 0 0 1 1 0 1
-    0 0 0 0 0 0 1 1 0
-    0 0 0 0 0 0 0 1 1
-    1 1 1 0 0 0 1 0 0
-    0 0 1 1 0 1 0 1 1
-    0 0 0 1 1 0 1 0 1
-    0 0 1 0 1 0 1 1 0
-]
 
-C = [
-    0 1 1 0 0 0
-    1 0 0 1 0 0
-    1 0 0 1 1 0
-    0 1 1 0 0 1
-    0 0 1 0 0 1
-    0 0 0 1 1 0
-]
-
-D = [
-    0 1 1 0 1 0
-    1 0 0 1 0 1
-    1 0 0 1 0 0
-    0 1 1 0 0 0
-    1 0 0 0 0 1
-    0 1 0 0 1 0]
-
-ZAA = zeros(Int64,10,10)
-ZAB = zeros(Int64,10,9)
-ZAC = zeros(Int64,10,6)
-ZAD = zeros(Int64,10,6)
-ZBB = zeros(Int64,9,9)
-ZBC = zeros(Int64,9,6)
-ZBD = zeros(Int64,9,6)
-ZCC = zeros(Int64,6,6)
-ZCD = zeros(Int64,6,6)
-
-ZBA = ZAB'
-ZCA = ZAC'
-ZDA = ZAD'
-ZCB = ZBC'
-ZDB = ZBD'
-ZDC = ZCD'
-
-
-# Graph 1 is disconnected, 2-3-4-5 are connected, 6 is disconnected, 7 is disconnected
-
-A = [
-    A   ZAA ZAB ZAC ZAA ZAA ZAD
-    ZAA A   ZAB ZAC ZAA ZAA ZAD
-    ZBA ZBA B   ZBC ZBA ZBA ZBD
-    ZCA ZCA ZCB C   ZCA ZCA ZCD
-    ZAA ZAA ZAB ZAC A   ZAA ZAD
-    ZAA ZAA ZAB ZAC ZAA A   ZAD
-    ZDA ZDA ZDB ZDC ZDA ZDA D
-]
-
-A[13,21] = A[21,13] = 1
-A[16,30] = A[30,16] = 1
-A[20,36] = A[36,20] = 1
-
-function initialize!_posdef!(system)
+function initialize!_posdef!(system::System{N}) where N 
     initialize!(system,rand)
-    for i=1:size(A)[1]
+    for i=1:N
         system.matrix_entries[i,i].value += 1000*I 
     end
 end
@@ -89,20 +16,17 @@ system = System{Float64}(A, ones(Int,size(A)[1])*3)
 systemldlt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
 systemllt = System{Float64}(A, ones(Int,size(A)[1])*3, symmetric=true)
 
-SUITE["ldu"] = @benchmarkable ldu_solve!($system) samples=2 setup=(initialize!($system))
-SUITE["lu"] = @benchmarkable lu_solve!($system) samples=2 setup=(initialize!($system))
-SUITE["ldlt"] = @benchmarkable ldlt_solve!($systemldlt) samples=2 setup=(initialize!($systemldlt))
-SUITE["llt"] = @benchmarkable llt_solve!($systemllt) samples=2 setup=(initialize!_posdef!($systemllt))
-
-# A = [
-#     0 1 1 1 1
-#     1 0 1 1 1
-#     1 1 0 1 1
-#     1 1 1 0 1
-#     1 1 1 1 0
-# ]
-
-
-
-
-
+SUITE["sparse_ldu"] = @benchmarkable ldu_solve!($system) setup=(initialize!($system))
+SUITE["sparse_lu"] = @benchmarkable lu_solve!($system) setup=(initialize!($system))
+SUITE["dense_lu"] = @benchmarkable lu(F)\f setup=(initialize!($system);F=full_matrix($system);f=full_vector($system))
+SUITE["sparse_ldlt"] = @benchmarkable ldlt_solve!($systemldlt) setup=(initialize!($systemldlt))
+SUITE["dense_ldlt"] = @benchmarkable bunchkaufman(F)\f setup=(initialize!($systemldlt);F=full_matrix($systemldlt);f=full_vector($systemldlt))
+SUITE["sparse_llt"] = @benchmarkable llt_solve!($systemllt) setup=(initialize!_posdef!($systemllt))
+SUITE["dense_llt"] = @benchmarkable cholesky(F)\f setup=(initialize!_posdef!($systemllt);F=full_matrix($systemllt);f=full_vector($systemllt))
+
+SUITE["sparse_add"] = @benchmarkable +($system,$system) setup=(initialize!($system))
+SUITE["dense_add"] = @benchmarkable (+(F,F);+(f+f)) setup=(initialize!($system);F=full_matrix($system);f=full_vector($system))
+SUITE["sparse_mul"] = @benchmarkable *($system,$system) setup=(initialize!($system))
+SUITE["dense_mul"] = @benchmarkable *(F,F) setup=(initialize!($system);F=full_matrix($system))
+SUITE["sparse_solve"] = @benchmarkable \($system,Bmat) setup=(initialize!($system);Bmat=deepcopy(system.matrix_entries);initialize!($system))
+SUITE["dense_solve"] = @benchmarkable \(F1,F2) setup=(initialize!($system);F2=full_matrix($system);initialize!($system);F1=full_matrix($system))

src/GraphBasedSystems.jl

@@ -2,6 +2,8 @@ module GraphBasedSystems
 
 using LinearAlgebra
 using SparseArrays
+using SparseArrays: widelength
+using SparseArrays.HigherOrderFns: _sumnnzs, _allocres, _map_zeropres!, _densestructure!
 using StaticArrays
 using Graphs
 
@@ -19,15 +21,19 @@ export System,
     parents,
 
     ldu_solve!,
+    ldu_matrix_solve!,
     ldu_factorization!,
     ldu_backsubstitution!,
     lu_solve!,
+    lu_matrix_solve!,
     lu_factorization!,
     lu_backsubstitution!,
     ldlt_solve!,
+    ldlt_matrix_solve!,
     ldlt_factorization!,
     ldlt_backsubstitution!,
     llt_solve!,
+    llt_matrix_solve!,
     llt_factorization!,
     llt_backsubstitution!
 
@@ -36,11 +42,12 @@ include(joinpath("util", "custom_static.jl"))
 
 include(joinpath("system", "entry.jl"))
 include(joinpath("system", "system.jl"))
-include(joinpath("system", "setup_functions.jl"))
+include(joinpath("system", "graph_functions.jl"))
 
 include(joinpath("system", "interface.jl"))
 include(joinpath("system", "dense.jl"))
 
+include(joinpath("solvers", "matrix.jl"))
 include(joinpath("solvers", "lu.jl"))
 include(joinpath("solvers", "llt.jl"))
 include(joinpath("solvers", "ldlt.jl"))

src/solvers/ldlt.jl

@@ -89,3 +89,12 @@ function ldlt_solve!(system::System)
     ldlt_backsubstitution!(system)
     return
 end  
+
+function ldlt_matrix_solve!(system::System, matrix::SparseMatrixCSC{Entry, Int64}; keep_vector = true)
+    keep_vector && (vector_entries = deepcopy(system.vector_entries))
+    ldlt_factorization!(system)
+    C = matrix_backsubsitution!(system, matrix, ldlt_backsubstitution!)
+    keep_vector && (system.vector_entries .= vector_entries)
+
+    return C
+end

src/solvers/ldu.jl

@@ -90,3 +90,12 @@ function ldu_solve!(system::System)
     ldu_backsubstitution!(system)
     return
 end  
+
+function ldu_matrix_solve!(system::System, matrix::SparseMatrixCSC{Entry, Int64}; keep_vector = true)
+    keep_vector && (vector_entries = deepcopy(system.vector_entries))
+    ldu_factorization!(system)
+    C = matrix_backsubsitution!(system, matrix, ldu_backsubstitution!)
+    keep_vector && (system.vector_entries .= vector_entries)
+
+    return C
+end

src/solvers/llt.jl

@@ -59,10 +59,7 @@ function llt_backsubstitution_d!(vector, diagonal, diagonal_inverse)
     if diagonal_inverse.isinverted
         vector.value = diagonal_inverse.value * vector.value
     else
-        invdiagonal = inv(diagonal.value)
-        diagonal_inverse.value = invdiagonal
-        diagonal_inverse.isinverted = true
-        vector.value = diagonal_inverse.value * vector.value
+        vector.value = diagonal.value \ vector.value
     end
     return
 end
@@ -98,3 +95,12 @@ function llt_solve!(system::System)
     llt_backsubstitution!(system)
     return
 end  
+
+function llt_matrix_solve!(system::System, matrix::SparseMatrixCSC{Entry, Int64}; keep_vector = true)
+    keep_vector && (vector_entries = deepcopy(system.vector_entries))
+    llt_factorization!(system)
+    C = matrix_backsubsitution!(system, matrix, llt_backsubstitution!)
+    keep_vector && (system.vector_entries .= vector_entries)
+
+    return C
+end

src/solvers/lu.jl

@@ -89,3 +89,12 @@ function lu_solve!(system::System)
     lu_backsubstitution!(system)
     return
 end
+
+function lu_matrix_solve!(system::System, matrix::SparseMatrixCSC{Entry, Int64}; keep_vector = true)
+    keep_vector && (vector_entries = deepcopy(system.vector_entries))
+    lu_factorization!(system)
+    C = matrix_backsubsitution!(system, matrix, lu_backsubstitution!)
+    keep_vector && (system.vector_entries .= vector_entries)
+
+    return C
+end

src/solvers/matrix.jl

@@ -0,0 +1,22 @@
+# Assumes quadratic matrices of same size
+function matrix_backsubsitution!(system::System{N}, matrix::SparseMatrixCSC{Entry, Int64}, backsubsitution!) where N
+    vector_entries = system.vector_entries
+    dims = system.dims
+
+    maxnnzC = Int(widelength(system.matrix_entries))
+    C = _allocres((N,N), Int64,  Entry, maxnnzC)
+    # _densestructure!(C)
+
+    for i = 1:N
+        for j = 1:N
+            Bji = matrix[j,i]
+            Bji isa Entry{Nothing} ? vector_entries[j] = Entry(dims[j], dims[i]) : vector_entries[j] = Bji
+        end
+        backsubsitution!(system)
+        for j = 1:N
+            C[j,i] = vector_entries[j]
+        end
+    end
+
+    return C
+end

src/system/dense.jl

@@ -1,19 +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)
+function full_matrix(matrix::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
+    for (i,row) in enumerate(matrix.rowval)
+        col = findfirst(x -> i<x, matrix.colptr)-1
+        A[range_dict_rows[row],range_dict_cols[col]] = matrix[row,col].value
         if symmetric && col != row
-            A[range_dict_cols[col],range_dict_rows[row]] = matrix_entries[row,col].value'
+            A[range_dict_cols[col],range_dict_rows[row]] = matrix[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]]...)
+full_vector(vector_entries::AbstractVector, dimensions_rows::SVector{N}) where N = vcat([getfield(vector_entries[i], :value) for i=1:N]...)