diff --git a/src/orbitdata.jl b/src/orbitdata.jl index d7bca32..1e302eb 100644 --- a/src/orbitdata.jl +++ b/src/orbitdata.jl @@ -107,11 +107,7 @@ function sparsify!(M::AbstractArray{T}, eps=eps(T); verbose=false) where T @info("Sparsifying $(size(M))-matrix... ") end - for n in eachindex(M) - if abs(M[n]) < eps - M[n] = zero(T) - end - end + clamp_small!(M, eps) if verbose @info("$(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(count(!iszero, M)) non-zeros)") @@ -120,6 +116,15 @@ function sparsify!(M::AbstractArray{T}, eps=eps(T); verbose=false) where T return sparse(M) end +function clamp_small!(M::AbstractArray{T}, eps=eps(T)) where T + for n in eachindex(M) + if abs(M[n]) < eps + M[n] = zero(T) + end + end + return M +end + function sparsify(U::AbstractArray{T}, tol=eps(T); verbose=false) where T return sparsify!(deepcopy(U), tol, verbose=verbose) end diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index be972ae..22796e9 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -96,8 +96,8 @@ function SOS_problem(X::GroupRingElem, orderunit::GroupRingElem, data::OrbitData end function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M)))) - for π in eachindex(Us) - M[π] = PropertyT.sparsify!(dims[π].*Ust[π]*cnstr*Us[π], eps) + for π in 1:endof(Us) + M[π] = dims[π].*PropertyT.clamp_small!(Ust[π]*cnstr*Us[π], eps) end end