Merge pull request #7 from kevin-tracy/master

Removing allocations in `update_values!` and adding tests

src/QDLDL.jl

@@ -49,7 +49,9 @@ struct QDLDLWorkspace{Tf<:AbstractFloat,Ti<:Integer}
     regularize_delta::Tf
 
     #number of regularized entries in D
-    regularize_count::Ref{Ti}
+    #length 1 vector instead of ref to avoid allocations
+    #while maintaining immutability
+    regularize_count::Vector{Ti}
 
 end
 
@@ -88,7 +90,7 @@ function QDLDLWorkspace(triuA::SparseMatrixCSC{Tf,Ti},
     positive_inertia = Ref{Ti}(-1)
 
     #number of regularized entries in D. None to start
-    regularize_count = Ref{Ti}(0)
+    regularize_count = zeros(Ti,1)
 
     QDLDLWorkspace(etree,Lnz,iwork,bwork,fwork,
                    Ln,Lp,Li,Lx,D,Dinv,positive_inertia,triuA,
@@ -199,7 +201,7 @@ function positive_inertia(F::QDLDLFactorisation)
 end
 
 function regularized_entries(F::QDLDLFactorisation)
-    F.workspace.regularize_count[]
+    F.workspace.regularize_count[1]
 end
 
 
@@ -378,14 +380,29 @@ end
 
 
 function QDLDL_factor!(
-        n,Ap,Ai,Ax,Lp,Li,Lx,
-        D,Dinv,Lnz,etree,bwork,iwork,fwork,
-        logicalFactor::Bool,Dsigns,
-        regularize_eps,regularize_delta,regularize_count
+        n,
+        Ap,
+        Ai,
+        Ax,
+        Lp,
+        Li,
+        Lx,
+        D,
+        Dinv,
+        Lnz,
+        etree,
+        bwork,
+        iwork,
+        fwork,
+        logicalFactor::Bool,
+        Dsigns,
+        regularize_eps,
+        regularize_delta,
+        regularize_count
 )
 
     positiveValuesInD  = 0
-    regularize_count[] = 0
+    regularize_count[1] = 0
 
     #partition working memory into pieces
     yMarkers        = bwork
@@ -415,7 +432,7 @@ function QDLDL_factor!(
         D[1]     = Ax[1]
         if(Dsigns !== nothing && Dsigns[1]*D[1] < regularize_eps)
             D[1] = regularize_delta * Dsigns[1]
-            regularize_count[] += 1
+            regularize_count[1] += 1
         end
 
         if(D[1] == 0.0) return -1 end
@@ -535,7 +552,7 @@ function QDLDL_factor!(
         #vector has been specified.
         if(Dsigns !== nothing && Dsigns[k]*D[k] < regularize_eps)
             D[k] = regularize_delta * Dsigns[k]
-            regularize_count[] += 1
+            regularize_count[1] += 1
         end
 
         #Maintain a count of the positive entries

test/UnitTests/update_values.jl

@@ -0,0 +1,36 @@
+using Test, LinearAlgebra, SparseArrays, Random
+using QDLDL
+rng = Random.MersenneTwister(1312)
+
+@testset "update values" begin
+
+    # create random KKT system
+    nz = 100
+    nc = 70
+    H = sprand(nz,nz,0.05);
+    H = H'*H + I
+    b = randn(nz + nc)
+    A1 = sprand(nc,nz,0.8)
+    K1 = [H A1';A1 -1e-3*I]
+
+    # get triu of K1 and create QDLDL struct
+    triuK1 = triu(K1);
+    F = qdldl(triuK1)
+
+    # compare with backslash
+    @test norm(K1\b - F\b,Inf) < 1e-12
+
+    # create a new KKT system with the same sparsity pattern
+    A2 = copy(A1)
+    A2.nzval .= randn(length(A2.nzval))
+    K2 = [H A2';A2 -1e-7*I]
+    triuK2 = triu(K2)
+
+    # update factorization of F in place (non allocating)
+    update_values!(F,1:length(triuK2.nzval),triuK2.nzval)
+    refactor!(F)
+
+    # compare with backslash using the refactored F
+    @test norm(K2\b - F\b,Inf)<1e-12
+
+end

test/runtests.jl

@@ -16,7 +16,7 @@ end
   include("./UnitTests/non-quasidef.jl")
   include("./UnitTests/inertia.jl")
   include("./UnitTests/regularization.jl")
-
+  include("./UnitTests/update_values.jl")
 
 end
 nothing