Skip to content

Latest commit

 

History

History
332 lines (248 loc) · 10.7 KB

File metadata and controls

332 lines (248 loc) · 10.7 KB
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

Build and Validate the Reference BLAS with PRIK

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.

What this example shows

  • 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.


Versions used

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.)


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
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 --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 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_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 complete test suite

Build both wrappers and run the 155-routine suite:

source examples/blas/build_all.sh
python3 -m pytest -q examples/blas/tests

The 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.


5. See how results are validated

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.

Test helper conventions

  • assert_allclose_for_dtype compares floating-point results with a tolerance matched to their NumPy dtype. Its optional operation_size is a rounding-error scale: use the number of terms in the calculation, such as 3 for the three products in the displayed dot product.
  • assert_storage_unchanged requires an input array to remain exactly equal to its saved value, including its dtype and any NaN sentinels.

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 focused examples

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

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 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.


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