require f2py to only wrap the necessary routines

Author: saidctb

.github/workflows/blas-lapack.yml

@@ -42,13 +42,15 @@ jobs:
             "meson==1.11.2" \
             "ninja==1.13.0" \
             "scipy==1.18.0"
-      - name: Install pinned GFortran
+      - name: Install pinned GFortran and LAPACK link dependencies
         shell: bash
         run: |
+          packages=(libblas-dev liblapack-dev)
           if ! command -v "$PRIK_GFORTRAN_BINARY" >/dev/null 2>&1; then
-            sudo apt-get update
-            sudo apt-get install --yes "$PRIK_GFORTRAN_PACKAGE"
+            packages+=("$PRIK_GFORTRAN_PACKAGE")
           fi
+          sudo apt-get update
+          sudo apt-get install --yes "${packages[@]}"
           compiler_dir="$RUNNER_TEMP/prik-gfortran"
           mkdir -p "$compiler_dir"
           ln -sf "$(command -v "$PRIK_GFORTRAN_BINARY")" "$compiler_dir/gfortran"

examples/lapack/README.md

@@ -35,6 +35,7 @@ same Meson, Ninja, and SciPy versions used by the dedicated BLAS/LAPACK CI lane
 before running this example:
 
 ```console
+sudo apt-get install libblas-dev liblapack-dev
 python3 -m pip install "meson==1.11.2" "ninja==1.13.0" "scipy==1.18.0"
 ```
 
@@ -47,8 +48,9 @@ It performs these operations once:
    cached shared library;
 3. build one PRIK extension from the complete generated contract and shared
    library;
-4. build one raw f2py comparison surface for the reviewed 127 routines, linked
-   to that same complete native library; and
+4. independently build one raw f2py comparison surface from the 125 reviewed
+   f2py-compatible routine sources and their minimal module dependency, linking unselected
+   LAPACK and BLAS helper symbols from the system development libraries; and
 5. reuse the three imported comparison surfaces for every correctness file.
 
 The contract-generation command is:
@@ -63,28 +65,52 @@ f2py command assembled by `conftest.py` is equivalent to:
 
 ```console
 python3 -m numpy.f2py -c -m f2py_reference_lapack_example \
-  <the 127 source files recorded in routine_inventory.py> \
-  <the cached complete libprik_full_lapack.so> \
+  examples/lapack/native/la_constants.f90 \
+  <the 125 f2py-compatible source files recorded in routine_inventory.py> \
+  only: <the 125 f2py-compatible routine names> : \
+  --dep lapack --dep blas \
+  --f2cmap /tmp/prik-lapack/.f2py_f2cmap \
   --build-dir /tmp/prik-lapack/f2py --f77flags=-O0 --f90flags=-O0 --opt=-O0
 ```
 
+`dlartg.f90` imports the `LA_CONSTANTS` module and declares its arguments as
+`REAL(wp)`. The f2py build therefore includes `la_constants.f90` as a build-only
+dependency and supplies a kind map that resolves `wp` to C `double`. The
+`only:` list makes the 125 f2py-compatible routine wrappers explicit; module support metadata
+is not counted as a selected routine. The dependency is compiled only to
+satisfy f2py's standalone selected-source build. f2py does not compile or link
+PRIK's complete 2,062-source native artifact. The reviewed implementations come
+from `examples/lapack/native/`; only their unselected transitive LAPACK and BLAS
+dependencies come from the system libraries.
+
 Build products stay in pytest temporary/cache directories. Failures report the
 compiler identity, command, stdout, and stderr. Neither build dirties the
 repository.
 
 The fixtures import the complete PRIK module, the raw f2py module, and
 `scipy.linalg.lapack`. Character flags are Python `str` through PRIK and
-`bytes` through SciPy/f2py. PRIK and raw f2py retain native argument ordering;
-SciPy projects arrays and optional arguments into its documented Python API.
+`bytes` through SciPy/f2py. PRIK retains the complete native argument order.
+Raw f2py retains the order of the arguments it exposes, but infers and hides
+leading dimensions and a few other shape-only scalar arguments; each visible
+f2py call follows the generated wrapper signature. SciPy projects arrays and
+optional arguments into its documented Python API.
 
 The Reference LAPACK subroutines do not declare Fortran `intent`, so raw f2py
 does not project their scalar writebacks. It still exposes in-place array
-mutation, which provides enough output for the independent oracle in 114 of
-the 123 subroutines; the four LAPACK functions return their values normally. Nine
-subroutines have an essential scalar-only result, or need an unprojected scalar
-to validate their mutated vector, and are explicitly recorded in
-`F2PY_NUMERICAL_LIMITATIONS`. All 127 selected names must still be exported by
-the one f2py comparison module.
+mutation, which provides enough output for the independent oracle in 112 of
+the 121 exported subroutines; the four exported LAPACK functions return their
+values normally. Nine subroutines have an essential scalar-only result, or need
+an unprojected scalar to validate their mutated vector, and are explicitly
+recorded in `F2PY_NUMERICAL_LIMITATIONS`.
+
+`dgees` and `dgges` remain in the 127-routine correctness inventory but are
+recorded in `F2PY_EXPORT_LIMITATIONS`. Their Reference LAPACK interfaces accept
+external selection callbacks. From the unannotated implementation sources,
+NumPy f2py 2.5.1 generates incomplete callback declarations (`select_t` and
+`selctg_t`) and invalid C wrapper code. Those two routines are therefore
+validated through PRIK, SciPy, and independent Schur reconstruction without a
+raw-f2py call; the remaining 125 selected names must be exported by the one
+f2py comparison module.
 
 The nine reviewed limitations are:
 
@@ -109,15 +135,28 @@ pivot ties may choose another valid permutation. Tests therefore use residuals,
 subspaces, factor reconstruction, and orthogonality instead of byte equality
 where the contract permits freedom.
 
-All numerical inputs use `dtype=np.float64`. `assert_allclose_float64` scales
-tolerance with double-precision epsilon and operation length. Residual checks
-also scale by matrix and solution norms. Integer results, `INFO`, shapes, dtypes,
-pivot bounds, and untouched sentinel storage use exact comparisons.
+All numerical arrays use `dtype=np.float64`. PRIK exposes the native scalar
+contract exactly, so its visible calls use `np.int32(...)` for LAPACK `INTEGER`
+arguments, `np.float64(...)` for scalar `DOUBLE PRECISION` arguments, and
+`np.bool_(...)` for scalar `LOGICAL` arguments. SciPy and f2py calls retain their
+own accepted Python conventions beside the PRIK call. `assert_allclose_float64`
+scales tolerance with double-precision epsilon and operation length. Residual
+checks also scale by matrix and solution norms. Integer results, `INFO`, shapes,
+dtypes, pivot bounds, and untouched sentinel storage use exact comparisons.
+
+`DTGSEN` and `DTRSEN` accept default-Fortran `LOGICAL` selection arrays. The
+pinned GFortran ABI stores each such element in four bytes, while PRIK's Python
+surface accepts a NumPy bool buffer. `gfortran_logical_mask` makes that genuine
+ABI difference explicit by placing each truth byte at the start of its native
+four-byte cell. The two affected calls remain visible, independently validate
+the requested reorder, and are recorded in `PRIK_ABI_ADAPTERS` rather than
+silently treating an ordinary two-byte NumPy mask as native storage.
 
 ## Arrays, storage, workspaces, and indexes
 
 Matrices are explicitly Fortran contiguous where the native contract requires
-column-major storage. Leading dimensions remain visible in PRIK/f2py calls.
+column-major storage. Leading dimensions remain visible in PRIK calls and are
+omitted only where the generated f2py signature explicitly infers them.
 Tests use NaNs or sentinel values in unused triangles and padding, then prove
 that those locations remain untouched. Helpers reconstruct general band,
 symmetric band, tridiagonal, packed, triangular, and rectangular-full-packed
@@ -130,9 +169,11 @@ Invalid argument calls that could route through `XERBLA` are not made in-process
 
 Native LAPACK pivots and positions are one-based. PRIK and the raw f2py surface
 preserve native values. SciPy converts several low-level pivot/index APIs to
-zero-based values; each affected test performs the conversion beside the call.
-The inventory marks every pivot/index-bearing routine so a new conversion
-cannot be hidden in a generic adapter.
+zero-based values, but preserves one-based `JPVT` for `DGELSY`/`DGEQP3`, one-based
+`IPIV` for the general-tridiagonal `DGTTRF`/`DGTTRS` family, and one-based Schur
+reorder positions for `DTGEXC`/`DTREXC`. Each affected test keeps the observed
+convention beside the call. The inventory marks every pivot/index-bearing
+routine so a new conversion cannot be hidden in a generic adapter.
 
 ## Run and diagnose
 
@@ -173,9 +214,11 @@ documented totals.
 | Explicit correctness tests | 127 |
 | PRIK exports required in CI | 127 |
 | SciPy exports used | 127 |
-| f2py exports required in CI | 127 |
-| Routines satisfying the independent oracle through f2py | 118 |
+| f2py exports required in CI | 125 |
+| Routines satisfying the independent oracle through f2py | 116 |
+| Documented raw-f2py export limitations | 2 |
 | Documented raw-f2py numerical projection limitations | 9 |
+| Documented PRIK default-LOGICAL ABI adapters | 2 |
 | Documented unsupported/skipped routines | 0 |
 
 Runtime success, f2py projection limitations, and any PRIK failures are CI
@@ -200,8 +243,10 @@ def test_dgesv_solves_general_system(prik_lapack, scipy_lapack, f2py_lapack):
     prik_piv = np.empty(2, dtype=np.int32)
     f2py_piv = np.empty(2, dtype=np.int32)
 
-    prik_scalars = prik_lapack.dgesv(2, 1, prik_a, 2, prik_piv, prik_b, 2, 0)
-    f2py_result = f2py_lapack.dgesv(2, 1, f2py_a, 2, f2py_piv, f2py_b, 2, 0)
+    prik_scalars = prik_lapack.dgesv(
+        np.int32(2), np.int32(1), prik_a, np.int32(2), prik_piv, prik_b, np.int32(2), np.int32(0)
+    )
+    f2py_result = f2py_lapack.dgesv(2, 1, f2py_a, f2py_piv, f2py_b, 0)
     scipy_lu, scipy_piv, scipy_x, scipy_info = scipy_lapack.dgesv(
         original_a.copy(order="F"), original_b.copy(order="F")
     )
@@ -234,8 +279,8 @@ def test_dgetrf_reconstructs_pivoted_lu(prik_lapack, scipy_lapack, f2py_lapack):
     prik_piv = np.empty(2, dtype=np.int32)
     f2py_piv = np.empty(2, dtype=np.int32)
 
-    prik_scalars = prik_lapack.dgetrf(2, 2, prik_a, 2, prik_piv, 0)
-    f2py_result = f2py_lapack.dgetrf(2, 2, f2py_a, 2, f2py_piv, 0)
+    prik_scalars = prik_lapack.dgetrf(np.int32(2), np.int32(2), prik_a, np.int32(2), prik_piv, np.int32(0))
+    f2py_result = f2py_lapack.dgetrf(2, 2, f2py_a, f2py_piv, 0)
     scipy_lu, scipy_piv, scipy_info = scipy_lapack.dgetrf(original.copy(order="F"))
 
     assert prik_scalars == (2, 2, 2, 0)
@@ -260,8 +305,8 @@ def test_dpotrf_reconstructs_spd_matrix(prik_lapack, scipy_lapack, f2py_lapack):
     stored = np.array([[4.0, np.nan], [1.0, 3.0]], dtype=np.float64, order="F")
     prik_a, f2py_a = column_major(stored), column_major(stored)
 
-    prik_scalars = prik_lapack.dpotrf("L", 2, prik_a, 2, 0)
-    f2py_result = f2py_lapack.dpotrf(b"L", 2, f2py_a, 2, 0)
+    prik_scalars = prik_lapack.dpotrf("L", np.int32(2), prik_a, np.int32(2), np.int32(0))
+    f2py_result = f2py_lapack.dpotrf(b"L", 2, f2py_a, 0)
     scipy_factor, scipy_info = scipy_lapack.dpotrf(stored.copy(order="F"), lower=1, clean=0)
 
     assert prik_scalars == (2, 2, 0)
@@ -289,8 +334,10 @@ def test_dgeqrf_reconstructs_qr_factorization(prik_lapack, scipy_lapack, f2py_la
     prik_tau = np.empty(2, dtype=np.float64)
     f2py_tau = np.empty(2, dtype=np.float64)
 
-    prik_scalars = prik_lapack.dgeqrf(3, 2, prik_a, 3, prik_tau, np.empty(16), 16, 0)
-    f2py_result = f2py_lapack.dgeqrf(3, 2, f2py_a, 3, f2py_tau, np.empty(16), 16, 0)
+    prik_scalars = prik_lapack.dgeqrf(
+        np.int32(3), np.int32(2), prik_a, np.int32(3), prik_tau, np.empty(16), np.int32(16), np.int32(0)
+    )
+    f2py_result = f2py_lapack.dgeqrf(3, 2, f2py_a, f2py_tau, np.empty(16), 16, 0)
     scipy_qr, scipy_tau, _scipy_work, scipy_info = scipy_lapack.dgeqrf(matrix.copy(order="F"), lwork=16)
 
     assert prik_scalars == (3, 2, 3, 16, 0)
@@ -318,8 +365,10 @@ def test_dsyev_returns_orthonormal_eigenvectors(prik_lapack, scipy_lapack, f2py_
     prik_w = np.empty(2, dtype=np.float64)
     f2py_w = np.empty(2, dtype=np.float64)
 
-    prik_scalars = prik_lapack.dsyev("V", "U", 2, prik_vectors, 2, prik_w, np.empty(16), 16, 0)
-    f2py_result = f2py_lapack.dsyev(b"V", b"U", 2, f2py_vectors, 2, f2py_w, np.empty(16), 16, 0)
+    prik_scalars = prik_lapack.dsyev(
+        "V", "U", np.int32(2), prik_vectors, np.int32(2), prik_w, np.empty(16), np.int32(16), np.int32(0)
+    )
+    f2py_result = f2py_lapack.dsyev(b"V", b"U", 2, f2py_vectors, f2py_w, np.empty(16), 16, 0)
     scipy_w, scipy_vectors, scipy_info = scipy_lapack.dsyev(matrix.copy(order="F"), compute_v=1, lower=0, lwork=16)
 
     assert prik_scalars == (2, 2, 16, 0)
@@ -346,8 +395,23 @@ def test_dgesvd_reconstructs_matrix(prik_lapack, scipy_lapack, f2py_lapack):
     prik_u, f2py_u = column_major(np.zeros((3, 3))), column_major(np.zeros((3, 3)))
     prik_vt, f2py_vt = column_major(np.zeros((2, 2))), column_major(np.zeros((2, 2)))
 
-    prik_scalars = prik_lapack.dgesvd("A", "A", 3, 2, prik_a, 3, prik_s, prik_u, 3, prik_vt, 2, np.empty(32), 32, 0)
-    f2py_result = f2py_lapack.dgesvd(b"A", b"A", 3, 2, f2py_a, 3, f2py_s, f2py_u, 3, f2py_vt, 2, np.empty(32), 32, 0)
+    prik_scalars = prik_lapack.dgesvd(
+        "A",
+        "A",
+        np.int32(3),
+        np.int32(2),
+        prik_a,
+        np.int32(3),
+        prik_s,
+        prik_u,
+        np.int32(3),
+        prik_vt,
+        np.int32(2),
+        np.empty(32),
+        np.int32(32),
+        np.int32(0),
+    )
+    f2py_result = f2py_lapack.dgesvd(b"A", b"A", 3, 2, f2py_a, f2py_s, f2py_u, f2py_vt, np.empty(32), 32, 0)
     scipy_u, scipy_s, scipy_vt, scipy_info = scipy_lapack.dgesvd(
         matrix.copy(order="F"), compute_uv=1, full_matrices=1, lwork=32
     )
@@ -379,8 +443,10 @@ def test_dtbtrs_solves_triangular_band_system(prik_lapack, scipy_lapack, f2py_la
     expected = np.array([[1.0], [2.0]], dtype=np.float64)
     prik_b, f2py_b = rhs.copy(order="F"), rhs.copy(order="F")
 
-    prik_scalars = prik_lapack.dtbtrs("U", "N", "N", 2, 1, 1, band, 2, prik_b, 2, 0)
-    f2py_result = f2py_lapack.dtbtrs(b"U", b"N", b"N", 2, 1, 1, band, 2, f2py_b, 2, 0)
+    prik_scalars = prik_lapack.dtbtrs(
+        "U", "N", "N", np.int32(2), np.int32(1), np.int32(1), band, np.int32(2), prik_b, np.int32(2), np.int32(0)
+    )
+    f2py_result = f2py_lapack.dtbtrs(b"U", b"N", b"N", 2, 1, 1, band, f2py_b, 0)
     scipy_x, scipy_info = scipy_lapack.dtbtrs(band, rhs.copy(order="F"), uplo=b"U", trans=b"N", diag=b"N")
 
     assert f2py_result is None
@@ -400,10 +466,18 @@ def test_dgecon_estimates_reciprocal_condition(prik_lapack, scipy_lapack, f2py_l
     expected = 1.0
 
     prik_scalars = prik_lapack.dgecon(
-        "1", 1, factor.copy(order="F"), 1, 4.0, 0.0, np.empty(4), np.empty(1, dtype=np.int32), 0
+        "1",
+        np.int32(1),
+        factor.copy(order="F"),
+        np.int32(1),
+        np.float64(4.0),
+        np.float64(0.0),
+        np.empty(4),
+        np.empty(1, dtype=np.int32),
+        np.int32(0),
     )
     f2py_result = f2py_lapack.dgecon(
-        b"1", 1, factor.copy(order="F"), 1, 4.0, 0.0, np.empty(4), np.empty(1, dtype=np.int32), 0
+        b"1", 1, factor.copy(order="F"), 4.0, 0.0, np.empty(4), np.empty(1, dtype=np.int32), 0
     )
     scipy_rcond, scipy_info = scipy_lapack.dgecon(factor.copy(order="F"), 4.0, norm=b"1")
 

examples/lapack/conftest.py

@@ -21,6 +21,9 @@
 NATIVE_ROOT = EXAMPLE_ROOT / "native"
 BUILD_FLAGS = "-O0"
 FORTRAN_SUFFIXES = (".f", ".f90", ".f95", ".f03", ".f08", ".for", ".f77", ".ftn")
+F2PY_BUILD_DEPENDENCIES = ("la_constants.f90",)
+F2PY_KIND_MAP = "{'real': {'wp': 'double'}}\n"
+F2PY_LINK_DEPENDENCIES = ("lapack", "blas")
 
 
 @dataclass(frozen=True)
@@ -91,6 +94,45 @@ def _selected_source(routine: str) -> Path:
     return matches[0]
 
 
+def _f2py_source_plan() -> tuple[Path, ...]:
+    """Return reviewed implementations plus their minimal compile dependency."""
+    dependencies = tuple(NATIVE_ROOT / name for name in F2PY_BUILD_DEPENDENCIES)
+    missing_dependencies = [str(path) for path in dependencies if not path.is_file()]
+    if missing_dependencies:
+        pytest.fail(f"missing f2py build dependencies: {missing_dependencies}")
+    selected = tuple(_selected_source(name) for name in ROUTINES if name not in F2PY_EXPORT_LIMITATIONS)
+    return dependencies + selected
+
+
+def _f2py_build_command(workdir: Path) -> tuple[str, ...]:
+    """Build only reviewed implementations and link external helper symbols."""
+    module_name = "f2py_reference_lapack_example"
+    f2cmap = workdir / ".f2py_f2cmap"
+    f2cmap.write_text(F2PY_KIND_MAP, encoding="utf-8")
+    selected_routines = tuple(name for name in ROUTINES if name not in F2PY_EXPORT_LIMITATIONS)
+    link_dependencies = tuple(item for dependency in F2PY_LINK_DEPENDENCIES for item in ("--dep", dependency))
+    return (
+        sys.executable,
+        "-m",
+        "numpy.f2py",
+        "-c",
+        "-m",
+        module_name,
+        *(str(source) for source in _f2py_source_plan()),
+        "only:",
+        *selected_routines,
+        ":",
+        *link_dependencies,
+        "--f2cmap",
+        str(f2cmap),
+        "--build-dir",
+        str(workdir / "generated"),
+        f"--f77flags={BUILD_FLAGS}",
+        f"--f90flags={BUILD_FLAGS}",
+        f"--opt={BUILD_FLAGS}",
+    )
+
+
 @pytest.fixture(scope="session")
 def prik_build(tmp_path_factory: pytest.TempPathFactory) -> BuiltLapack:
     """Build the one complete PRIK LAPACK wrapper through the established path."""
@@ -134,27 +176,11 @@ def prik_lapack(prik_build: BuiltLapack):
 
 @pytest.fixture(scope="session")
 def f2py_build(tmp_path_factory: pytest.TempPathFactory) -> BuiltLapack:
-    """Build one raw f2py comparison surface linked to the complete library."""
+    """Build one raw f2py surface from only the reviewed implementations."""
     compiler = full._compiler()
     workdir = tmp_path_factory.mktemp("f2py-reference-lapack-example")
-    shared = full._cached_native_shared_library("lapack")
     module_name = "f2py_reference_lapack_example"
-    selected_sources = tuple(_selected_source(name) for name in ROUTINES if name not in F2PY_EXPORT_LIMITATIONS)
-    command = (
-        sys.executable,
-        "-m",
-        "numpy.f2py",
-        "-c",
-        "-m",
-        module_name,
-        *(str(source) for source in selected_sources),
-        str(shared),
-        "--build-dir",
-        str(workdir / "generated"),
-        f"--f77flags={BUILD_FLAGS}",
-        f"--f90flags={BUILD_FLAGS}",
-        f"--opt={BUILD_FLAGS}",
-    )
+    command = _f2py_build_command(workdir)
     result = _run_build(command, workdir, compiler)
     return BuiltLapack(
         module=_import_built_module(module_name, workdir),

examples/lapack/helpers.py

@@ -53,6 +53,13 @@ def native_pivots(scipy_pivots: np.ndarray) -> np.ndarray:
     return np.asarray(scipy_pivots, dtype=np.int32) + np.int32(1)
 
 
+def gfortran_logical_mask(values) -> np.ndarray:
+    """Represent a default-GFortran LOGICAL vector through PRIK's bool buffer ABI."""
+    logical_bytes = np.zeros(len(values) * np.dtype(np.int32).itemsize, dtype=np.bool_)
+    logical_bytes[:: np.dtype(np.int32).itemsize] = np.asarray(values, dtype=np.bool_)
+    return logical_bytes
+
+
 def pivot_matrix(pivots: np.ndarray, size: int, *, one_based: bool) -> np.ndarray:
     """Build the row permutation represented by sequential LAPACK pivots."""
     permutation = np.eye(size, dtype=np.float64)