IPM/HSD: use logical slicing instead of elementwise multiplication

We now do basically this for a logical vector l: `dot(a[l], b[l])`, instead of `dot(a .* l, b .* l) as before. This is as suggested here: ds4dm#122 (comment) Helps with the performance of the BigFloat arithmetic. This is the output with the script from the previous commit: ``` exp = 6: 27.391877 seconds (142.35 M allocations: 21.646 GiB, 11.38% gc time, 23.70% compilation time) elapsed time (ns): 27391877044 gc time (ns): 3116697350 bytes allocated: 23242616545 pool allocs: 142136686 non-pool GC allocs: 9144 malloc() calls: 118007 realloc() calls: 88649 free() calls: 116981 minor collections: 506 full collections: 0 exp = 9: 289.050404 seconds (1.12 G allocations: 175.440 GiB, 29.73% gc time) elapsed time (ns): 289050404133 gc time (ns): 85940338982 bytes allocated: 188377360520 pool allocs: 1117528611 non-pool GC allocs: 75494 malloc() calls: 728937 realloc() calls: 554273 free() calls: 743507 minor collections: 4186 full collections: 3 ``` Fixes ds4dm#122
nsajko · Jan 10, 2023 · b41ea93 · b41ea93
1 parent 4273f5c
commit b41ea93
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 12 deletions.
diff --git a/src/IPM/HSD/HSD.jl b/src/IPM/HSD/HSD.jl
@@ -203,8 +203,8 @@ function compute_residuals!(hsd::HSD{T}
         (
             (dat.c, pt.x),
             (dat.b, pt.y),
-            (dat.l .* dat.lflag, pt.zl),
-            (dat.u .* dat.uflag, pt.zu)),
+            (dat.l[dat.lflag], pt.zl[dat.lflag]),
+            (dat.u[dat.uflag], pt.zu[dat.uflag])),
         (
             1, -1, -1, 1))
     )
@@ -222,8 +222,8 @@ function compute_residuals!(hsd::HSD{T}
         dot_buf,
         (
             (dat.b, pt.y),
-            (dat.l .* dat.lflag, pt.zl),
-            (dat.u .* dat.uflag, pt.zu)),
+            (dat.l[dat.lflag], pt.zl[dat.lflag]),
+            (dat.u[dat.uflag], pt.zu[dat.uflag])),
         (
             1, 1, -1)) / pt.τ + dat.c0
 
@@ -295,8 +295,8 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
         dot_buf,
         (
             (dat.b, pt.y),
-            (dat.l .* dat.lflag, pt.zl),
-            (dat.u .* dat.uflag, pt.zu)),
+            (dat.l[dat.lflag], pt.zl[dat.lflag]),
+            (dat.u[dat.uflag], pt.zu[dat.uflag])),
         (
             1, 1, -1)) * ϵi
 

diff --git a/src/IPM/HSD/step.jl b/src/IPM/HSD/step.jl
@@ -69,8 +69,8 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
     h0 = buffered_dot_weighted_sum!!(
         dot_buf,
         (
-            (dat.l .* dat.lflag, (dat.l .* θl) .* dat.lflag),
-            (dat.u .* dat.uflag, (dat.u .* θu) .* dat.uflag),
+            (dat.l[dat.lflag], (dat.l .* θl)[dat.lflag]),
+            (dat.u[dat.uflag], (dat.u .* θu)[dat.uflag]),
             ((@. (c + (θl * dat.l) * dat.lflag + (θu * dat.u) * dat.uflag)), hx),
             (b, hy)),
         (
@@ -224,10 +224,10 @@ function solve_newton_system!(Δ::Point{T, Tv},
         buffered_dot_weighted_sum!!(
             dot_buf,
             (
-                ((ξxzl ./ pt.xl) .* dat.lflag, dat.l .* dat.lflag),            # l'(Xl)^-1 * ξxzl
-                ((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag),
-                (((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag),
-                (((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag)),
+                ((ξxzl ./ pt.xl)[dat.lflag], dat.l[dat.lflag]),            # l'(Xl)^-1 * ξxzl
+                ((ξxzu ./ pt.xu)[dat.uflag], dat.u[dat.uflag]),
+                (((pt.zl ./ pt.xl) .* ξl)[dat.lflag], dat.l[dat.lflag]),
+                (((pt.zu ./ pt.xu) .* ξu)[dat.uflag], dat.u[dat.uflag])),
             (
                 -1, 1, -1, -1))