equalize() Parlett-Reinsch balancing diverges to Inf

[jamestjat]

Summary

equalize() in minimal.go:219 (the Parlett-Reinsch state-space balancing routine invoked via Reduce(&ReduceOpts{Equalize: true})) diverges to Inf on at least one realistic input. After one call with the 3×3 closed-loop from issue #26, A[0,1] becomes -Inf, B[0,0] becomes 1e+186, and B[2,0] becomes 5e+302. Any downstream use of the system then panics or produces NaN.

Reproducer

// Build the issue #26 PI+plant closed loop: clToPv has
// A = [[0, -0.56, 0]; [4.19e-9, -4.327e-5, 0]; [4.19e-9, -2.82e-5, -1.51e-5]]
// B = [1; 5.03e-5; 5.03e-5]
// C = [[0, 0, 0.56]]
// D = [0]
mr, _ := clToPv.Reduce(&controlsys.ReduceOpts{Mode: controlsys.ReduceAll, Equalize: true})
// mr.Order == 0; mr.Sys is a bogus 1×1 System before PR #26 fix,
// or a clean empty System after. In both cases, the A/B/C matrices
// that equalize produced internally were garbage (Inf/1e302).

Root cause

Three independent bugs in equalize() vs. reference LAPACK DGEBAL:

1. minimal.go:269 and :278 — c *= sclfac * sclfac should be c *= sclfac (LAPACK updates c and r linearly per iteration, not quadratically in c alone).
2. Missing r update — LAPACK updates r /= sclfac (first loop) and r *= sclfac (second loop) per iteration; the Go loops at minimal.go:264-279 never touch r, so the balance comparison uses stale values.
3. minimal.go:281 — convergence check (c+r)/f >= factor*s divides by f; LAPACK uses (c+r) >= factor*s with c, r already adjusted in the loops. The spurious /f inverts the check once f grows large, falsely signaling "converged" while f keeps exploding through the outer iteration.
4. minimal.go:288-300 — the scaling application direction looks inverted vs. LAPACK: LAPACK does A[i,:] /= f and A[:,i] *= f (shrink row, grow column when r > c); the Go code does A[i,:] *= f and A[:,i] /= f, amplifying the imbalance.

Impact

Low (for now): no code in the repo passes Equalize: true, and MinimalRealization() uses nil opts. This is a latent bug waiting for a user to enable balancing and get silently corrupted results.

Blocks: matching SLICOT AB01ND / python-control minreal behavior, which pre-balance by default for better numerical conditioning.

Proposed fix

Rewrite equalize() to follow LAPACK DGEBAL exactly, with the SLICOT TB01ID extension for including B rows and C columns in the row/column norms. Add a regression test that balances the issue #26 system and verifies the output has finite entries and the same eigenvalues as the input (similarity transform invariance).

Discovered while

Investigating #26 (MinimalRealization over-reduction). Not the root cause of #26 itself — that was a separate fallthrough bug in Reduce fixed in a6cf478.

[jamestjsp]

Fixed in v0.8.0 via 2ec61ea. equalize() rewritten to match LAPACK DGEBAL + SLICOT TB01ID.