diff --git a/src/QDLDL.jl b/src/QDLDL.jl index e80f09c..0ed0d9c 100644 --- a/src/QDLDL.jl +++ b/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 diff --git a/test/UnitTests/update_values.jl b/test/UnitTests/update_values.jl new file mode 100644 index 0000000..d1bb498 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 96209ac..d59824b 100644 --- a/test/runtests.jl +++ b/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