IPM/HSD: for BigFloat optimization, computing ξg_ in solve_newton_system! causes a big portion of total allocations

This line causes a lot of allocation, reported by --track-allocation=all:

Line 216 in 2b055fd

@timeit hsd.timer "ξg_" ξg_ = (ξg + ξtk / pt.τ

So there's possibly room for improvement at or around that line.

I tried a couple of different syntaxes, and it looks like, when using BigFloat arithmetic, things like

dot((ξxzl ./ pt.xl) .* dat.lflag, dat.l .* dat.lflag)

allocate a lot of memory indeed.
IIRC, I chose that syntax because then I didn't have to worry about dimensions (everything has the same size, and what needs to be zero is always zero).

Based on your comment, I think a better design choice would be to drop these .* flag and slice the arrays instead.
Something like

dot( (ξxzl ./ pt.xl)[dat.lflag], dat.l[dat.lflag])

A quick thought, maybe useful: I think that part of the reason for the current inefficiency is the bool-float multiplication method from base/bool.jl Base.:*(x::Bool, y::T)::promote_type(Bool,T) where T<:AbstractFloat, which takes special care for the case when the float value is NaN or infinity. I'm not sure how does Tulip handle NaN internally, but an idea that comes to mind is this:
If it's known that y (of floating point type) is finite and not NaN, it would probably be faster to convert x from Bool to UInt8 or UInt before multiplying.

The main reason for using bool-float multiplication was because some of these can be Inf or NaN, yes. They get multiplied by false and get reset to zero, which makes the math work.

I did not realize at the time it could cause big performance issues.

I started optimizing IPM/HSD using MutableArithmetics with good results, however I wonder which arithmetics are supposed to be supported for HSD? It's likely that I'll break some non-BigFloat arithmetics while developing, so I'm wondering how much do the current tests cover this?