Fix matrix inverse for varying dims

Author: janbruedigam

src/GraphBasedSystems.jl

@@ -47,6 +47,7 @@ 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/matrix.jl

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

src/system/system.jl

@@ -139,27 +139,4 @@ function Base.:/(S1::System, r::Real)
     return S2
 end
 Base.:\(r::Real, S::System) = S/r
-Base.:\(A::System, B::SparseMatrixCSC{Entry, Int64}) = ldu_matrix_solve!(A, B; keep_vector = true)
-
-# Assumes quadratic matrices of same size
-function matrix_backsubsitution!(systemA::System{N}, B::SparseMatrixCSC{Entry, Int64}, backsubsitution!) where N
-    vector_entries = systemA.vector_entries
-    dims = systemA.dims
-
-    maxnnzC = Int(widelength(systemA.matrix_entries))
-    C = _allocres((N,N), Int64,  Entry, maxnnzC)
-    _densestructure!(C)
-
-    for i = 1:N
-        for j = 1:N
-            Bji = B[j,i]
-            Bji isa Entry{Nothing} ? vector_entries[j] = Entry(dims[i], dims[j]) : vector_entries[j] = Bji
-        end
-        backsubsitution!(systemA)
-        for j = 1:N
-            C[j,i] = vector_entries[j]
-        end
-    end
-
-    return C
-end
+Base.:\(A::System, B::SparseMatrixCSC{Entry, Int64}) = lu_matrix_solve!(A, B; keep_vector = true)

test/ldlt_test.jl

@@ -15,7 +15,7 @@ for i=1:5
 end
 
 for i=1:5
-    system = System{Float64}(A, ones(Int64, size(A)[1]), symmetric=true)
+    system = System{Float64}(A, rand(0:3, size(A)[1]), symmetric=true)
     initialize!(system)
     Bmat = deepcopy(system.matrix_entries)
     F2 = full_matrix(Bmat,false,system.dims,system.dims)

test/ldu_test.jl

@@ -15,7 +15,7 @@ for i=1:5
 end
 
 for i=1:5
-    system = System{Float64}(A, ones(Int64, size(A)[1]))
+    system = System{Float64}(A, rand(0:3, size(A)[1]))
     initialize!(system)
     Bmat = deepcopy(system.matrix_entries)
     F2 = full_matrix(Bmat,false,system.dims,system.dims)

test/llt_test.jl

@@ -22,7 +22,7 @@ for i=1:5
 end
 
 for i=1:5
-    system = System{Float64}(A, ones(Int64, size(A)[1]), symmetric=true)
+    system = System{Float64}(A, rand(0:3, size(A)[1]), symmetric=true)
     initialize!(system)
     Bmat = deepcopy(system.matrix_entries)
     F2 = full_matrix(Bmat,false,system.dims,system.dims)

test/lu_test.jl

@@ -15,7 +15,7 @@ for i=1:5
 end
 
 for i=1:5
-    system = System{Float64}(A, ones(Int64, size(A)[1]))
+    system = System{Float64}(A, rand(0:3, size(A)[1]))
     initialize!(system)
     Bmat = deepcopy(system.matrix_entries)
     F2 = full_matrix(Bmat,false,system.dims,system.dims)
@@ -26,6 +26,18 @@ for i=1:5
     @test maximum(abs.(full_matrix(Cmat,false,system.dims,system.dims)-F1\F2)) < 1e-3
 end
 
+for i=1:5
+    system = System{Float64}(A, rand(0:3, size(A)[1]))
+    initialize!(system)
+    Bmat = deepcopy(system.matrix_entries)
+    F2 = full_matrix(Bmat,false,system.dims,system.dims)
+    initialize!(system)
+    F1 = full_matrix(system)
+
+    Cmat = system \ Bmat
+    @test maximum(abs.(full_matrix(Cmat,false,system.dims,system.dims)-F1\F2)) < 1e-3
+end
+
 system = System{Float64}(A, rand(0:3, size(A)[1]))
 
 display(system)