fn balance(
operator: &mut DMatrix<Complex>,
row: &mut [f64],
column: &mut [f64],
) -> Result<()>Expand description
Even out the row and column sums of a max-equilibrated operator in place with
power-of-two factors, folded into row and column.
Row and column maxima do not determine the scales: in A = D_r H D_c the largest
entry of a row depends on D_c, so an off-diagonal entry of H can end up as
large as its diagonal, partial pivoting may choose it, and a change of units costs
digits. Unit magnitude sums determine the scaled operator whenever it has total
support (Sinkhorn and Knopp), whatever the units.
Each sweep moves every column factor, then every row factor, to a power of two that
minimizes the convex potential Σ |a_ij| x_i y_j − Σ ln x_i − Σ ln y_j along it. The
sweeps stop after one that moves no factor by more than two, or after
BALANCE_SWEEPS. A sum that is not normal, or factors beyond 2^±256, keep the
max scales.