Build and Validate the Reference BLAS with PRIK¶

This example shows how to take the official Reference BLAS sources and produce two importable Python extension modules:

• one generated by PRIK
• one generated by NumPy’s f2py

It then tests both wrappers against explicit mathematical expectations. The full test suite in the repository applies the same method to all 155 callable routines in the Reference BLAS corpus.

Why this example exists¶

• It demonstrates PRIK on a real Fortran library.
• It makes the differences between the two wrappers visible.
• It verifies numerical behaviour with an independent mathematical check (instead of only checking that the two wrappers agree).

You should already be comfortable with NumPy arrays, basic packaging, and building Fortran extensions.

---

Versions used by the maintained example¶

Component	Version / source
PRIK	current repository checkout (0.1.0)
Reference BLAS	snapshot shipped in Netlib LAPACK 3.12.1
Python	3.12 (dedicated CI job)
NumPy / f2py	NumPy 2.5.1
Meson	1.11.2
Ninja	1.13.0
Fortran compiler	GNU Fortran 13 in CI (any compatible gfortran works locally)

Note: f2py is part of NumPy. On Python 3.12 it uses the Meson backend, which is why Meson and Ninja are required.

For everyday use of this example, prefer the checked-in sources in examples/blas/native/. (See the Source provenance section at the end if you want to verify the upstream archive yourself.)

---

1. Prepare the repository and toolchain¶

Clone PRIK, create a virtual environment, and install the same Python build tools used by the dedicated CI job:

git clone https://github.com/PyNumLab/prik.git
cd prik
python3 -m venv .venv
. .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install -e ".[qa]" \
  "numpy==2.5.1" "meson==1.11.2" "ninja==1.13.0"

Install GNU Fortran separately. On Ubuntu:

sudo apt-get update
sudo apt-get install --yes gfortran
gfortran --version

All remaining commands run from the repository root with the virtual environment active.

The runnable material is self-contained in the repository's examples/ directory. After PRIK and the listed tools are installed, you can copy that directory alone.

---

2. Compile BLAS once and build the PRIK wrapper¶

Run the first build script from the repository root:

export EXAMPLE_WORKSPACE="$PWD"
export BLAS_BUILD_ROOT="$(mktemp -d)"
export BLAS_SHARED_LIBRARY="$(
  python -m examples.native_library blas \
    --compiler "$(command -v gfortran)" \
    --jobs 8
)"

mkdir -p "$BLAS_BUILD_ROOT/prik/generated"
cd "$BLAS_BUILD_ROOT/prik"

python -m prik "$EXAMPLE_WORKSPACE/examples/blas/native" \
  --out prik_reference_blas \
  --out-dir "$BLAS_BUILD_ROOT/prik/generated" \
  --compiler "$(command -v gfortran)" \
  --no-compile-input-sources \
  --native-objects "$BLAS_SHARED_LIBRARY" \
  --jobs 8 \
  --wrapper-fortran-flags="-O0 -g0" \
  --wrapper-c-flags="-O0 -g0"

examples.native_library compiles all 155 implementations and returns the resulting shared-library path. PRIK reads the same source directory to build the Python API, skips native implementation compilation, and links that library.

-O0 keeps the PRIK and f2py correctness builds equivalent and avoids making optimization-dependent claims. This example focuses on correctness, not performance.

---

3. Build f2py against the same native library¶

Run the same f2py build script exercised by the test suite:

cd "$EXAMPLE_WORKSPACE"
export BLAS_F2PY_ROOT="$BLAS_BUILD_ROOT/f2py"
mkdir -p "$BLAS_F2PY_ROOT/generated"
cd "$BLAS_F2PY_ROOT"

export FC="$(command -v gfortran)"
export F77="$FC"
export F90="$FC"
export FFLAGS="-O0"
export F90FLAGS="-O0"
export LDFLAGS="${LDFLAGS:+$LDFLAGS }-Wl,-rpath,$(dirname "$BLAS_SHARED_LIBRARY")"

python -m numpy.f2py -c \
  "$EXAMPLE_WORKSPACE/examples/blas/blas.pyf" \
  "-L$(dirname "$BLAS_SHARED_LIBRARY")" \
  -lprik_full_blas \
  --build-dir "$BLAS_F2PY_ROOT/generated" \
  --f77flags=-O0 \
  --f90flags=-O0 \
  --opt=-O0

The committed blas.pyf is the reviewed f2py interface. f2py compiles only its wrapper and links it to BLAS_SHARED_LIBRARY, so both wrappers exercise the same compiled BLAS implementations.

Six rotation routines document scalar writebacks without declaring Fortran intent. Their reviewed intent(inout) declarations live directly in blas.pyf, and the tests pass typed 0-D arrays. PRIK needs neither: it returns unannotated scalar writebacks directly, with ordinary scalar arguments.

Import both modules:

import os
import sys

build_root = os.environ["BLAS_BUILD_ROOT"]
sys.path.insert(0, f"{build_root}/prik")
sys.path.insert(0, f"{build_root}/f2py")

import f2py_reference_blas
import prik_reference_blas

Important wrapper difference¶

PRIK deliberately follows the native scalar contract:

• For a subroutine such as DAXPY, arrays are mutated in place and PRIK returns the visible input scalars because they are treated as inout.
• The f2py comparison module also mutates the output array but returns None.
• Function routines such as DDOT return their numerical result through both wrappers.

---

4. Run the correctness tests¶

The names prik_blas and f2py_blas in the tests are session-scoped pytest fixtures from conftest.py. After the build scripts finish, they import the modules and reuse them for every test under examples/blas/tests/.

The tests import small, explicit helpers from tests/helpers.py. The two helpers used below are intentionally narrow:

• assert_allclose_for_dtype chooses a tolerance from the result dtype and the number of accumulated operations. It delegates the final comparison to numpy.testing.assert_allclose.
• assert_storage_unchanged uses exact array equality, including NaNs, to prove that an input-only array or padding was not modified.

Their complete implementations are short enough to show here. Both helper definitions and the test functions below assume:

import numpy as np

def assert_allclose_for_dtype(actual, expected, *, operation_size: int = 1) -> None:
    """Compare floating values with dtype- and accumulation-aware tolerances."""
    actual_array = np.asarray(actual)
    expected_array = np.asarray(expected)
    dtype = np.result_type(actual_array.dtype, expected_array.dtype)
    real_dtype = np.empty((), dtype=dtype).real.dtype
    epsilon = np.finfo(real_dtype).eps
    scale = max(1, operation_size)
    magnitude = max(1.0, float(np.max(np.abs(expected_array), initial=0.0)))
    np.testing.assert_allclose(
        actual_array,
        expected_array,
        rtol=epsilon * 8 * scale,
        atol=epsilon * 8 * scale * magnitude,
    )

def assert_storage_unchanged(actual: np.ndarray, original: np.ndarray) -> None:
    """Require exact preservation, including NaNs and sentinel padding."""
    np.testing.assert_array_equal(actual, original, strict=True)

---

5. Validate behaviour, not just compilation¶

A solid numerical test checks three relationships:

PRIK result      == independent mathematical result
f2py result      == independent mathematical result
PRIK result      == f2py result

It should also verify mutation, input preservation, dtype, shape, increments, leading dimensions, and unused storage when those properties are part of the routine contract.

f2py provides useful differential evidence, but the independent formula or residual is the primary oracle.

The two examples below are taken verbatim from the runnable suite. A documentation test compares these blocks with the Python AST, so the page and the real tests stay in sync.

DAXPY – in-place vector update¶

def test_daxpy(prik_blas, f2py_blas):
    alpha = np.float64(-1.5)
    x = np.array([2.0, -4.0, 1.0], dtype=np.float64)
    original_y = np.array([3.0, 5.0, -2.0], dtype=np.float64)
    prik_x, f2py_x = x.copy(), x.copy()
    prik_y, f2py_y = original_y.copy(), original_y.copy()

    prik_scalars = prik_blas.daxpy(np.int32(3), alpha, prik_x, np.int32(1), prik_y, np.int32(1))
    f2py_result = f2py_blas.daxpy(np.int32(3), alpha, f2py_x, np.int32(1), f2py_y, np.int32(1))

    expected_y = alpha * x + original_y
    assert_allclose_for_dtype(prik_y, expected_y)
    assert_allclose_for_dtype(f2py_y, expected_y)
    assert_allclose_for_dtype(prik_y, f2py_y)
    assert prik_scalars == (np.int32(3), alpha, np.int32(1), np.int32(1))
    assert f2py_result is None
    assert_storage_unchanged(prik_x, x)
    assert_storage_unchanged(f2py_x, x)

Both wrappers must mutate y to the expected value. The input-only array x must remain unchanged.

DDOT – scalar function result¶

def test_ddot(prik_blas, f2py_blas):
    x = np.array([1.0, -2.0, 4.0], dtype=np.float64)
    y = np.array([3.0, 5.0, -1.0], dtype=np.float64)
    prik_x, f2py_x = x.copy(), x.copy()
    prik_y, f2py_y = y.copy(), y.copy()

    prik_value, n, incx, incy = prik_blas.ddot(np.int32(3), prik_x, np.int32(1), prik_y, np.int32(1))
    f2py_value = f2py_blas.ddot(np.int32(3), f2py_x, np.int32(1), f2py_y, np.int32(1))

    expected = np.float64(1.0 * 3.0 + (-2.0) * 5.0 + 4.0 * (-1.0))
    assert_allclose_for_dtype(prik_value, expected, operation_size=3)
    assert_allclose_for_dtype(f2py_value, expected, operation_size=3)
    assert_allclose_for_dtype(prik_value, f2py_value, operation_size=3)
    assert (n, incx, incy) == (np.int32(3), np.int32(1), np.int32(1))
    assert_storage_unchanged(prik_x, x)
    assert_storage_unchanged(f2py_x, x)
    assert_storage_unchanged(prik_y, y)
    assert_storage_unchanged(f2py_y, y)

---

6. Run the maintained example¶

Build both wrappers once, then run any user-facing test selection:

cd "$REPOSITORY_ROOT"
source examples/blas/build_all.sh
python -m pytest -q examples/blas/tests

Focused commands for quick debugging:

python -m pytest -q examples/blas/tests/test_level1_real.py
python -m pytest -q examples/blas/tests/test_level1_real.py::test_daxpy
python -m pytest -q examples/blas/tests -k dgemm

• Complete Level-1 examples → test_level1_real.py
• Matrix / packed / banded / symmetric / Hermitian / triangular examples → files under examples/blas/tests/
• Authoritative classification → routine_inventory.py
• Coverage guard → test_routine_coverage.py

For the copyable build scripts, test commands, and source provenance, see the examples/blas project README.

---

Troubleshooting¶

• Confirm that gfortran, meson and ninja are on your PATH.
• On Python 3.12+, do not force the old distutils backend of f2py. Use the pinned Meson + Ninja setup shown above.
• Run a single failing test with more detail and keep the build directory:

bash python -m pytest -vv -s --basetemp=/tmp/prik-blas-debug examples/blas/tests/test_level1_real.py::test_daxpy

• Read the compiler output from build_all.sh.
• Keep correctness tests and benchmarking completely separate. This suite uses small deterministic inputs and makes no performance claims.

---

Source provenance¶

The files under examples/blas/native/ are byte-for-byte copies of the 155 files in BLAS/SRC/ from the official LAPACK 3.12.1 archive.

If you want to reconstruct the upstream sources yourself:

curl --location --output lapack-3.12.1.tar.gz \
  https://www.netlib.org/lapack/lapack-3.12.1.tar.gz

printf '%s  %s\n' \
  37b00c90947488521f475b5a187fff4da4a5cfe61b525efcacf7a97f39a45ec6 \
  lapack-3.12.1.tar.gz | sha256sum --check -

tar -xzf lapack-3.12.1.tar.gz

Official license and provenance: Netlib LAPACK site · LAPACK license