| title | Build and Validate the Reference BLAS with PRIK |
|---|---|
| audience | users, advanced users |
| prerequisites | arrays, packaging |
| related | lapack-wrapper.md, ../guide/arrays.md |
| status | maintained |
| publication | reviewed |
This example builds the official Reference BLAS sources as two importable Python extension modules:
- one generated by PRIK
- one generated by NumPy’s f2py
It compares both wrappers with independent mathematical results across all 155 callable Reference BLAS routines.
- Build PRIK and f2py wrappers against the same compiled BLAS library.
- Call vector and matrix routines with NumPy arrays.
- Compare numerical results and the Python interfaces produced by each tool.
You should already be comfortable with NumPy arrays, basic packaging, and building Fortran extensions.
| Component | Version / source |
|---|---|
| PRIK | current repository checkout |
| Reference BLAS | snapshot shipped in Netlib LAPACK 3.12.1 |
| Python | 3.12 in the dedicated CI job |
| NumPy / f2py | NumPy 2.5.1 |
| Meson | 1.11.2 |
| Ninja | 1.13.0 |
| Fortran compiler | GNU Fortran 13 in CI; a 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.)
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
python3 -m pip install --upgrade pip
python3 -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 --versionAll 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.
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.
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=-O0The committed blas.pyf defines the f2py
interface. f2py compiles only its wrapper and links it to
BLAS_SHARED_LIBRARY, so both wrappers exercise the same compiled BLAS
implementations.
A few routines expose scalar writebacks differently through the two wrappers; the comparison below shows those return-value differences explicitly.
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_blasPRIK 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
DDOTreturn their numerical result through both wrappers.
Build both wrappers and run the 155-routine suite:
source examples/blas/build_all.sh
python3 -m pytest -q examples/blas/testsThe tests cover vector, matrix, packed, banded, symmetric, Hermitian, and triangular operations. Each routine is called with representative inputs and checked against an independent mathematical result.
Each comparison checks three relationships:
PRIK result == independent mathematical result
f2py result == independent mathematical result
PRIK result == f2py result
The suite also checks mutation, input preservation, dtype, shape, increments, leading dimensions, and unused storage where they are part of a routine's contract. The independent formula or residual remains the primary numerical reference.
The two examples below come directly from the runnable suite and use its small NumPy comparison helpers.
assert_allclose_for_dtypecompares floating-point results with a tolerance matched to their NumPy dtype. Its optionaloperation_sizeis a rounding-error scale: use the number of terms in the calculation, such as3for the three products in the displayed dot product.assert_storage_unchangedrequires an input array to remain exactly equal to its saved value, including its dtype and anyNaNsentinels.
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.
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)After building the wrappers, run a family or one routine:
python3 -m pytest -q examples/blas/tests/test_level1_real.py
python3 -m pytest -q examples/blas/tests/test_level1_real.py::test_daxpy
python3 -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/ - Public routine list →
routine_inventory.py - Routine coverage check →
test_routine_coverage.py
For the copyable build scripts, test commands, and source provenance, see the
examples/blas project README.
-
Confirm that
gfortran,mesonandninjaare on yourPATH. -
On Python 3.12+, do not force the old distutils backend of f2py. Use the pinned Meson and Ninja setup shown above.
-
Run a single failing test with more detail and keep the build directory:
python3 -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.
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.gzOfficial license and provenance: Netlib LAPACK site · LAPACK license