Close translation parity gaps

.github/workflows/test.yml

@@ -24,6 +24,8 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
+        with:
+          submodules: recursive
 
       - name: Set up Python
         uses: actions/setup-python@v5

.gitignore

@@ -65,3 +65,4 @@ _codeql_detected_source_root
 .agent/rules/beads-workflow-context.md
 *.txt
 .cache/
+.ergo/lock

docs/incomplete_translations.md

@@ -1,162 +1,61 @@
-# Incomplete SLICOT C Translations
+# Current SLICOT Translation Gap Audit
 
-Systematic audit of ~620 C source files against Fortran77 originals.
-Generated 2025-02-07.
+Refreshed 2026-05-10 against the current C sources and
+`SLICOT-Reference/src`. The older 2025 line-ratio table is no longer reliable:
+several files it called stubs now contain substantial translated logic.
 
 ## Summary
 
-| Category | Count | Action |
-|----------|-------|--------|
-| INCOMPLETE | 8 | Must fix — wrong/degraded results |
-| STUB | 2 | Must fix — mostly unimplemented |
-| DEGRADED | 5 | Should fix — correct but missing optimizations |
-| FALSE_POSITIVE | ~30 | No action needed |
-
----
-
-## INCOMPLETE — Major algorithmic sections missing
-
-### ab13hd.c (4% ratio, 141/3439 lines)
-**Missing:** 33 subroutine calls. Entire JOBE='G'/'C' (general/companion) code paths.
-Only handles JOBE='I' (identity). Missing: DGEBAL, DGEHRD, DHSEQR, DGESVD, DGETRF,
-DGETRS, DGGBAK, DGGBAL, DHGEQZ, DLASRT, DORMHR, DORMQR, DSYRK, DTGSEN, DTRCON,
-DTRSM, MA02AD, MB01SD, MB02RD, MB02SD, MB02TD, MB03XD, MB04BP, SB04OD, TB01ID,
-TG01AD, TG01BD, ZGESVD. The file explicitly notes "simplified implementation".
-**Impact:** Returns wrong results for non-identity E matrices.
-
-### ab09iy.c (34% ratio, 290/844 lines)
-**Missing:** MB01WD, DSYEV (spectral factorization path for frequency-weighted
-Gramians when alpha < 1). The basic weighted/unweighted Lyapunov paths are present
-but the spectral reweighting for stability-enforcing alpha parameters is absent.
-**Impact:** Wrong Gramians when ALPHAC < 1 or ALPHAO < 1 with non-trivial weights.
-
-### ab13dd.c (66% ratio, 1241/1853 lines)
-**Missing:** DLASRT, MB01SD. Missing eigenvalue sorting and symmetric matrix update
-in the continuous-time structured singular value computation.
-**Impact:** May produce incorrect mu upper bounds for continuous-time systems.
-
-### mb01vd.c (7% ratio, 133/1678 lines)
-**Missing:** 28-branch sparse optimization using DCOPY-based operations for all
-UPLO/TRANS/TYPE combinations. Only basic dense path implemented.
-**Impact:** Correct for dense matrices but may produce wrong results for sparse
-structured operations where the Fortran has specialized handling.
-
-### mb02pd.c (22% ratio, 124/540 lines)
-**Missing:** DGEEQU (equilibration), DGERFS (iterative refinement), DLAQGE
-(row/column scaling). Error bounds (ferr/berr) are zeroed instead of computed.
-Equilibration (FACT='E') is silently ignored.
-**Impact:** Reduced numerical accuracy. No iterative refinement means larger
-backward errors for ill-conditioned systems.
-
-### mb03bd.c (68% ratio, 1335/1954 lines)
-**Missing:** Zero-chasing algorithm for periodic Schur decomposition (TODO marker at
-line 1246). The main eigenvalue reordering logic exists but the ILO/IHI zero-chasing
-sweep path is unimplemented.
-**Impact:** May fail or produce wrong results for certain eigenvalue reordering cases.
-
-### mb04rb.c (23% ratio, 80/347 lines)
-**Missing:** MB01KD, MB04PA, DGEMM. Only handles basic case, missing the blocked
-symplectic URV decomposition update algorithm.
-**Impact:** Wrong results for non-trivial block operations.
-
-### mb01uy.c (34% ratio, 158/457 lines)
-**Missing:** DGEMM, DGEQRF. The Fortran has optimized paths using DGEMM for large
-matrices and DGEQRF for QR-based symmetric rank-k update. C version only has
-element-wise loops.
-**Impact:** Correct for small matrices but missing the DGEMM/DGEQRF optimized path
-that handles the actual computation differently for larger cases. Needs verification
-that element-wise path covers all cases.
-
----
-
-## STUB — Only parameter validation, most logic missing
-
-### mb04pb.c (7% ratio, 27/341 lines)
-**Missing:** DGEHRD, DGEMM, DSYR2K, MB04PA. Only has parameter validation.
-Requires UE01MD (not yet translated) for blocked symplectic factorization.
-**Impact:** Function does nothing. Blocked on UE01MD dependency.
-
-### mb04tb.c (12% ratio, 89/687 lines)
-**Missing:** DGEMM, DGEQRF, MB03XU. Delegates to mb04ts for unblocked case only.
-The entire blocked algorithm (the primary purpose of the routine) is missing.
-**Impact:** Only works for small matrices where unblocked path suffices.
-
----
-
-## DEGRADED — Correct but missing optimizations
-
-### mb03xp.c (16% ratio, 105/633 lines)
-**Missing:** MB03YA, DLARNV, DLARFG, DLARFX, DGEMV, DTRMV (multishift QZ bulge
-chasing). Delegates entirely to mb03yd (single-shift QZ).
-**Impact:** Correct results but O(n) slower for large matrices. The multishift
-optimization is ~300 lines of Fortran.
-
-### mb02cu.c (71% ratio, 713/1001 lines)
-**Missing:** DGELQ2, DGEQR2 (unblocked QR/LQ). Uses DGEQRF/DGELQF (blocked
-versions) instead.
-**Impact:** Correct. Blocked versions are actually preferred for large matrices.
-Minor: unblocked versions may be slightly faster for very small matrices.
-
-### mb03bb.c (120% ratio, 495/411 lines)
-**Missing:** DLADIV, ZLARTG, ZROT. These are inlined as direct computations.
-C version is actually longer than Fortran (complex arithmetic expanded).
-**Impact:** Correct. The complex operations are correctly implemented inline.
-
-### mc01td.c (46% ratio, 136/290 lines)
-**Missing:** DRSCL (reciprocal scaling). Uses inline `1.0/scale` multiplication
-instead.
-**Impact:** Correct. DRSCL is a trivial BLAS-like operation.
-
-### mb03bf.c (64% ratio, 84/130 lines)
-**Missing:** DROT (Givens rotation application). Uses inline rotation code.
-**Impact:** Correct. DROT is trivially inlined.
-
----
-
-## FALSE_POSITIVE — No issues found
-
-| File | Reason |
-|------|--------|
-| ab05rd.c | Full translation, low ratio due to verbose Fortran comments |
-| ab09ed.c | Full translation, calls all required subroutines |
-| ab08nz.c | Full translation with all complex variants |
-| ab08mz.c | Full translation with all complex variants |
-| dg01nd.c | Legitimate wrapper delegating to dg01md + dg01ny |
-| fd01ad.c | Full translation with DLARTG inlined via SLC_DLARTG |
-| ib01ad.c | Dispatcher calling ib01md/ib01od, all paths present |
-| ib01px.c | Full translation, low ratio due to Fortran comment blocks |
-| ib01py.c | Full translation, low ratio due to Fortran comment blocks |
-| ib03ad.c | Full translation, "simplified" comment refers to workspace calc only |
-| mb01ru.c | Full translation with DGEMM/DSYR2K inlined |
-| mb03md.c | Full translation (bisection algorithm fully implemented) |
-| mb03nd.c | Full translation (small routine, QR iteration) |
-| mb03pd.c | Full translation (product eigenvalue computation) |
-| mb03vd.c | Full translation (sorting eigenvalues of product) |
-| mb04nd.c | Full translation (RQ/QR factorization update) |
-| mb04od.c | Full translation (RQ factorization update) |
-| mb04ty.c | Full translation (symplectic butterfly) |
-| mc03nd.c | Full translation (polynomial GCD computation) |
-| nf01bq.c | Full translation, parameter validation is complete |
-| sb02ru.c | Calls mb02pd (which itself is incomplete, but sb02ru's translation is faithful) |
-| sb03ou.c | Full translation (Cholesky factor of Lyapunov solution) |
-| sb16cd.c | Full translation (coprime factorization controller) |
-| sg03bd.c | Uses DGGES (modern replacement for deprecated DGEGS) |
-| tb01ux.c | Full translation delegating to tb01ud + tb01xd |
-| tb04ad.c | Full translation (transfer matrix computation) |
-| tf01qd.c | Full translation (output response, simple loop) |
-| tf01rd.c | Full translation (output response, simple loop) |
-| tg01hd.c | Full translation (descriptor system reduction) |
-| tg01jy.c | Uses tg01hx instead of tg01hy (equivalent replacement) |
-| tg01nd.c | Full translation (descriptor system reduction) |
-
----
-
-## Notes
-
-- **Method:** Automated triage via subroutine call comparison + line-count ratio,
-  followed by manual verification of each candidate (reading both C and Fortran sources).
-- **Filters applied:** Inlined BLAS (DCOPY, DSCAL, DSWAP, DAXPY, DLASET, DLACPY,
-  DLASCL, DROT, DLADIV, DRSCL), legitimate delegation patterns, modern LAPACK
-  replacements (DGEGS→DGGES).
-- **mb01vd** sparse optimization is the largest single gap (~1300 lines of Fortran).
-- **ab13hd** is the most severe — 96% of Fortran logic missing.
+| Classification | Count | Meaning |
+| --- | ---: | --- |
+| Active parity task | 0 | No current behavior gap remains open in the refreshed audit. |
+| Needs differential proof | 0 | The risky paths identified by this audit now have focused proof or a task result documenting why no fix was needed. |
+| Completed / stale entry | 13 | The missing calls, stub claims, or risky paths from the old audit are now closed by implementation, differential coverage, or branch-specific proof. |
+| Deliberate replacement / inlining | 5 | The C code intentionally uses equivalent wrappers, expanded arithmetic, or inline loops. |
+| Performance-only degradation | 1 | Behavior is expected to be correct; remaining gap is optimization. |
+
+## Closed parity tasks
+
+| Routine | Current status | Tracking |
+| --- | --- | --- |
+| `ab13dd.c` | The continuous-time candidate-frequency ordering path has branch coverage and Fortran-reference proof. | `LZ6EVQ` |
+| `mb03bd.c` | The `DLAS2` singularity-test concern is covered by focused parity proof. | `KJDQQC` |
+| `mb01uy.c` | Workspace-query semantics now match the Fortran `DGEQRF` query contract. | `XDMHPE` |
+| `mb02cu.c` | Blocked QR/LQ paths and blocked `B` reflector-scalar side effects now have parity evidence. | `2ZMHMW`, `VERVRD` |
+
+## Differential proof added
+
+| Routine | Current finding |
+| --- | --- |
+| `ab13hd.c` | The old "4% identity-only" claim is stale; current C contains the broader descriptor/generalized code paths and `DLASRT` calls, with targeted coverage and Fortran differential tests for risky branches. |
+| `ab09iy.c` | The old missing `MB01WD`/`DSYEV` claim is stale; current C calls both and now has focused weighted-Gramian differential coverage. |
+| `mb01vd.c` | Sparse/dense Kronecker product branches now have branch-matrix differential coverage across transpose combinations. |
+| `mb04rb.c` | The blocked symplectic update concern has been audited and covered by a narrow parity task. |
+| `mb04tb.c` | The blocked path has Fortran differential coverage; the old stub/delegation claim is stale. |
+
+## Completed or stale entries
+
+| Routine | Current finding |
+| --- | --- |
+| `mb02pd.c` | Current C calls `DGEEQU`, `DLAQGE`, and `DGERFS`; the old equilibration/refinement gap is stale. |
+| `mb04pb.c` | Current C calls `DGEHRD`, `MB04PA`, `DGEMM`, `DSYR2K`, and `MB04PU`; the old validation-only stub claim is stale. |
+| `ab09iy.c` | No longer missing `MB01WD` or `DSYEV`; retained above only because differential proof is still useful. |
+| `mb04tb.c` | No longer delegates entirely to `mb04ts`; retained above only because blocked-path proof is still useful. |
+
+## Deliberate replacement / inlining
+
+| Routine | Current finding |
+| --- | --- |
+| `mb03bb.c` | `DLADIV`, `ZLARTG`, and `ZROT` are expanded inline as complex arithmetic. |
+| `mc01td.c` | `DRSCL` is replaced with explicit reciprocal scaling. |
+| `mb03bf.c` | `DROT` is implemented inline. |
+| `sg03bd.c` | Uses modern `DGGES` instead of deprecated `DGEGS`. |
+| `tg01jy.c` | Uses the current helper path instead of the older Fortran helper name. |
+
+## Current gates
+
+- `scripts/check_translation_smells.py`
+- `scripts/check_translation_coverage.py --report-only`
+- `python3 -m py_compile scripts/check_translation_smells.py scripts/check_translation_coverage.py`
+
+Unresolved questions: none

include/slicot_blas.h

@@ -957,6 +957,9 @@ void SLC_FC_FUNC(dtgsen, DTGSEN)(const int* ijob, const int* wantq, const int* w
 void SLC_FC_FUNC(dgeqr2, DGEQR2)(const int* m, const int* n, f64* a,
                                   const int* lda, f64* tau, f64* work, int* info);
 
+void SLC_FC_FUNC(dgelq2, DGELQ2)(const int* m, const int* n, f64* a,
+                                  const int* lda, f64* tau, f64* work, int* info);
+
 void SLC_FC_FUNC(dorg2r, DORG2R)(const int* m, const int* n, const int* k,
                                   f64* a, const int* lda, const f64* tau,
                                   f64* work, int* info);
@@ -1068,6 +1071,7 @@ void SLC_FC_FUNC(dlasrt, DLASRT)(const char* id, const int* n, f64* d, int* info
 #define SLC_DSYTRS   SLC_FC_FUNC(dsytrs, DSYTRS)
 #define SLC_DLANSY   SLC_FC_FUNC(dlansy, DLANSY)
 #define SLC_DGEQRF   SLC_FC_FUNC(dgeqrf, DGEQRF)
+#define SLC_DGELQ2   SLC_FC_FUNC(dgelq2, DGELQ2)
 #define SLC_DGELQF   SLC_FC_FUNC(dgelqf, DGELQF)
 #define SLC_DGEQP3   SLC_FC_FUNC(dgeqp3, DGEQP3)
 #define SLC_DGERQF   SLC_FC_FUNC(dgerqf, DGERQF)

python/wrappers/py_mb01.c

@@ -714,8 +714,10 @@ PyObject* py_mb01uy(PyObject* self, PyObject* args) {
     PyObject *t_obj, *a_obj;
     PyArrayObject *t_array, *a_array;
     i32 info;
+    i32 ldwork = 0;
+    Py_ssize_t nargs = PyTuple_Size(args);
 
-    if (!PyArg_ParseTuple(args, "sssiidOO", &side, &uplo, &trans, &m, &n, &alpha, &t_obj, &a_obj)) {
+    if (!PyArg_ParseTuple(args, "sssiidOO|i", &side, &uplo, &trans, &m, &n, &alpha, &t_obj, &a_obj, &ldwork)) {
         return NULL;
     }
 
@@ -753,8 +755,12 @@ PyObject* py_mb01uy(PyObject* self, PyObject* args) {
         }
     }
 
-    i32 ldwork = m * n;
-    f64 *dwork = (f64*)malloc(ldwork * sizeof(f64));
+    if (nargs <= 8) {
+        ldwork = m * n;
+    }
+    i32 dwork_len = (ldwork == -1) ? 1 : ldwork;
+    if (dwork_len < 1) dwork_len = 1;
+    f64 *dwork = (f64*)malloc((size_t)dwork_len * sizeof(f64));
     if (dwork == NULL) {
         Py_DECREF(t_array);
         Py_DECREF(a_array);
@@ -763,9 +769,18 @@ PyObject* py_mb01uy(PyObject* self, PyObject* args) {
     }
 
     mb01uy(side, uplo, trans, m, n, alpha, t_data, ldt, a_data, lda, dwork, ldwork, &info);
+    f64 work = dwork[0];
 
     free(dwork);
 
+    if (ldwork == -1) {
+        PyObject *result = Py_BuildValue("Oid", Py_None, info, work);
+        Py_DECREF(t_array);
+        Py_DECREF(a_array);
+        Py_DECREF(t_out);
+        return result;
+    }
+
     npy_intp result_dims[2] = {m, n};
     npy_intp result_strides[2] = {sizeof(f64), m * sizeof(f64)};
     PyObject *result_array = PyArray_New(&PyArray_Type, 2, result_dims, NPY_DOUBLE, result_strides,
@@ -784,7 +799,12 @@ PyObject* py_mb01uy(PyObject* self, PyObject* args) {
         }
     }
 
-    PyObject *result = Py_BuildValue("Oi", result_array, info);
+    PyObject *result;
+    if (nargs <= 8) {
+        result = Py_BuildValue("Oi", result_array, info);
+    } else {
+        result = Py_BuildValue("Oid", result_array, info, work);
+    }
     Py_DECREF(t_array);
     Py_DECREF(a_array);
     Py_DECREF(t_out);

python/wrappers/py_mb02.c

@@ -1337,15 +1337,15 @@ PyObject* py_mb02cu(PyObject* self, PyObject* args, PyObject* kwargs) {
         }
     } else if (lcol) {
         if (col2 > 0) {
-            lcs = 5 * k + (k < q ? k : q);
+            lcs = 6 * k;
         } else {
-            lcs = 2 * k + (k < q ? k : q);
+            lcs = 4 * k;
         }
     } else {
         if (col2 > 0) {
-            lcs = 5 * k + (k < q ? k : q);
+            lcs = 6 * k;
         } else {
-            lcs = 2 * k + (k < q ? k : q);
+            lcs = 4 * k;
         }
     }
     if (lcs < 1) lcs = 1;

scripts/check_translation_coverage.py

@@ -24,18 +24,18 @@
 DEFAULT_TARGETS = [
     ("src/AB/ab13hd.c", 1842, "AB13HD C-row access branch"),
     ("src/MB/mb04hd.c", 228, "MB04HD reorder/BWORK call"),
-    ("src/MB/mb04cd.c", 758, "MB04CD first-stage A column save"),
-    ("src/MB/mb04cd.c", 760, "MB04CD first-stage A column update"),
-    ("src/MB/mb04cd.c", 788, "MB04CD first-stage B column save"),
-    ("src/MB/mb04cd.c", 790, "MB04CD first-stage B column update"),
-    ("src/MB/mb04cd.c", 817, "MB04CD first-stage D column save"),
-    ("src/MB/mb04cd.c", 819, "MB04CD first-stage D column update"),
-    ("src/MB/mb04cd.c", 844, "MB04CD first-stage Q1 column save"),
-    ("src/MB/mb04cd.c", 846, "MB04CD first-stage Q1 column update"),
-    ("src/MB/mb04cd.c", 1559, "MB04CD second-stage A column update"),
-    ("src/MB/mb04cd.c", 1592, "MB04CD second-stage B column update"),
-    ("src/MB/mb04cd.c", 1625, "MB04CD second-stage D column update"),
-    ("src/MB/mb04cd.c", 1659, "MB04CD second-stage Q1 column update"),
+    ("src/MB/mb04cd.c", 774, "MB04CD first-stage A column save"),
+    ("src/MB/mb04cd.c", 776, "MB04CD first-stage A column update"),
+    ("src/MB/mb04cd.c", 804, "MB04CD first-stage B column save"),
+    ("src/MB/mb04cd.c", 806, "MB04CD first-stage B column update"),
+    ("src/MB/mb04cd.c", 833, "MB04CD first-stage D column save"),
+    ("src/MB/mb04cd.c", 835, "MB04CD first-stage D column update"),
+    ("src/MB/mb04cd.c", 860, "MB04CD first-stage Q1 column save"),
+    ("src/MB/mb04cd.c", 862, "MB04CD first-stage Q1 column update"),
+    ("src/MB/mb04cd.c", 1589, "MB04CD second-stage A column update"),
+    ("src/MB/mb04cd.c", 1622, "MB04CD second-stage B column update"),
+    ("src/MB/mb04cd.c", 1655, "MB04CD second-stage D column update"),
+    ("src/MB/mb04cd.c", 1676, "MB04CD second-stage Q1 column update"),
 ]