Skip to content

tb01pd: API design causes incorrect results with pre-padded arrays #12

Closed as not planned

Description

@jamestjsp

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

  1. C code is correct - Direct ctypes calls to tb01pd with explicit dimensions work correctly
  2. Python wrapper infers dimensions - py_tb.c line 775-776:
    m = (i32)b_dims[1];  // infers m from B columns
    p = (i32)c_dims[0];  // infers p from C rows
  3. Pre-padding breaks inference - When B has shape (n, maxmp), wrapper thinks m=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 values

Option 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

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.

Activity

  1. added a commit that references this issue on Feb 3, 2026
    2443fb9
  2. jamestjsp commented on Feb 6, 2026

    @jamestjsp
    OwnerAuthor

    Closing as won't-fix. Our API is correct by design — we infer dimensions from array shapes, which is the Pythonic approach. Pass correctly-shaped arrays (e.g., B with shape (n, m), C with shape (p, n)) and the wrapper works correctly. Pre-padding arrays to (n, max(m,p)) is a caller bug, not a wrapper bug. The python-control fix (python-control/python-control#1200) addresses this on their side.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions