tb01pd: API design causes incorrect results with pre-padded arrays #12
Description
Summary
When users pass pre-padded arrays to tb01pd, slicot.c produces incorrect results because it infers m and p from array shapes rather than accepting them as explicit parameters.
Related: python-control/python-control#1200
Root Cause
| Wrapper | API | m, p handling |
|---|---|---|
| slycot | tb01pd(n, m, p, A, B, C, ...) |
User specifies explicitly |
| slicot.c | tb01pd(job, equil, A, B, C, tol) |
Inferred from array shapes |
When a user passes a pre-padded B array with shape (n, max(m,p)) instead of (n, m), our wrapper incorrectly infers m = max(m,p).
Minimal Reproduction
import numpy as np
from slicot import tb01pd as slicot_tb01pd
# System with m=1 input, p=2 outputs
n, m_actual, p_actual = 3, 1, 2
maxmp = max(m_actual, p_actual)
a = np.array([
[1.0, 2.0, 0.0],
[4.0, -1.0, 0.0],
[0.0, 0.0, 1.0]
], dtype=float, order='F')
# CORRECT: Pass actual-sized arrays
b_correct = np.array([[1.0], [0.0], [1.0]], dtype=float, order='F') # shape (3, 1)
c_correct = np.array([
[0.0, 1.0, -1.0],
[0.0, 0.0, 1.0]
], dtype=float, order='F') # shape (2, 3)
_, _, _, nr1, _, info1 = slicot_tb01pd('M', 'N', a.copy(), b_correct.copy(), c_correct.copy(), 0.0)
print(f"Correct arrays: nr={nr1}") # nr=3 ✓
# WRONG: Pre-padded arrays (what python-control was doing)
b_padded = np.zeros((n, maxmp), dtype=float, order='F') # shape (3, 2)
b_padded[:, 0] = [1.0, 0.0, 1.0]
b_padded[:, 1] = [1e100, 1e100, 1e100] # garbage in workspace
c_padded = np.zeros((maxmp, n), dtype=float, order='F') # shape (2, 3)
c_padded[0, :] = [0.0, 1.0, -1.0]
c_padded[1, :] = [0.0, 0.0, 1.0]
_, _, _, nr2, _, info2 = slicot_tb01pd('M', 'N', a.copy(), b_padded.copy(), c_padded.copy(), 0.0)
print(f"Pre-padded arrays: nr={nr2}") # nr=2 ✗ (wrapper infers m=2 from b_padded.shape[1])Output:
Correct arrays: nr=3
Pre-padded arrays: nr=2
Comparison with slycot
slycot works correctly because the user explicitly passes m=1, p=2:
from slycot import tb01pd as slycot_tb01pd
# slycot API takes explicit n, m, p
ar, br, cr, nr = slycot_tb01pd(n, m_actual, p_actual, a.copy(), b_padded.copy(), c_padded.copy())
print(f"slycot with pre-padded: nr={nr}") # nr=3 ✓Investigation Details
- C code is correct - Direct ctypes calls to
tb01pdwith explicit dimensions work correctly - Python wrapper infers dimensions -
py_tb.cline 775-776:m = (i32)b_dims[1]; // infers m from B columns p = (i32)c_dims[0]; // infers p from C rows
- Pre-padding breaks inference - When B has shape
(n, maxmp), wrapper thinksm=maxmp
Proposed Solutions
Option 1: Add optional m, p parameters (recommended)
# New API with backwards compatibility
tb01pd(job, equil, a, b, c, tol, m=None, p=None)
# If m/p not provided, infer from shapes (current behavior)
# If provided, use explicit valuesOption 2: Document the limitation
Document that users should NOT pre-pad arrays - pass actual-sized B(n,m) and C(p,n).
Option 3: Match slycot API exactly
tb01pd(n, m, p, a, b, c, job='M', equil='N', tol=0.0)Related Files
- python-control minreal implementation: https://github.com/python-control/python-control/blob/main/control/statesp.py#L1175
- python-control fix (empty→zeros): https://github.com/python-control/python-control/pull/1200/files
- slycot tb01pd wrapper: https://github.com/python-control/Slycot/blob/master/slycot/transform.py
- slycot f2py interface: https://github.com/python-control/Slycot/blob/master/slycot/src/transform.pyf
Note on the python-control bug
The original python-control bug used np.empty() instead of np.zeros() for workspace padding. While slycot tolerated this due to explicit m/p parameters, our wrapper would fail regardless because of the dimension inference issue. The empty() vs zeros() issue is secondary to the API design difference.