dsymm_bug.py

gistfile1.txt

"""Minimal reproducer: OpenBLAS 0.3.20 DSYMM is wrong on its Cooperlake kernel.

C = B @ A with A symmetric, column-major, via
    cblas_dsymm(CblasColMajor, CblasRight, CblasUpper,
                m, n, 1.0, A(n x n), n, B(m x n), m, 0.0, C, m)
compared against numpy. (This is the call OpenMEEG makes in
Matrix::operator*(const SymMatrix&).)

Needs only numpy and a libopenblas. To exercise a kernel your CPU does not
have, use Intel SDE -- `-cpx` sets CPUID to Cooper Lake -- together with
OpenBLAS's own dispatch override:

    python3 dsymm_bug.py                                        # native kernel
    OPENBLAS_CORETYPE=Cooperlake sde64 -cpx -- python3 dsymm_bug.py

Observed on an i7-7700K (no AVX-512), so every Cooperlake row below is emulated:

    Ubuntu 22.04, libopenblas 0.3.20+ds-1
      kernel Haswell      max rel err 1.0e-15   ok
      kernel Cooperlake   max rel err 1.4e+00   WRONG, for every shape tried

    Ubuntu 24.04, libopenblas 0.3.26+ds-1ubuntu0.1
      kernel Haswell      max rel err 9.7e-16   ok
      kernel Cooperlake   max rel err 1.1e-15   ok

    OpenBLAS 0.3.32
      kernel Cooperlake   max rel err 1.1e-15   ok

So it is fixed somewhere between 0.3.20 and 0.3.26. Impact: any runner whose
CPU dispatches 0.3.20 to Cooperlake silently returns a badly wrong DSYMM.

Exits non-zero if a wrong result is detected.
"""

import ctypes
import os
import sys

import numpy as np

CBLAS_COL_MAJOR, CBLAS_RIGHT, CBLAS_UPPER = 102, 142, 121


def find_openblas():
    for cand in (
        "/usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so",
        "/usr/lib/x86_64-linux-gnu/libopenblas.so.0",
        "libopenblas.so.0",
    ):
        try:
            return ctypes.CDLL(cand), cand
        except OSError:
            continue
    raise SystemExit("no libopenblas found")


def main():
    dll, path = find_openblas()
    try:
        dll.openblas_get_corename.restype = ctypes.c_char_p
        core = dll.openblas_get_corename().decode()
    except AttributeError:
        core = "?"
    print(f"library : {os.path.basename(path)}")
    print(f"kernel  : {core}")

    rng = np.random.default_rng(0)
    worst = 0.0
    for m, n in ((32, 1126), (32, 512), (32, 200), (8, 1126), (1, 1126)):
        A = rng.standard_normal((n, n))
        A = (A + A.T) / 2.0  # symmetric
        B = rng.standard_normal((m, n))
        Af = np.asfortranarray(A)
        Bf = np.asfortranarray(B)
        C = np.zeros((m, n), order="F")
        dll.cblas_dsymm(
            ctypes.c_int(CBLAS_COL_MAJOR), ctypes.c_int(CBLAS_RIGHT),
            ctypes.c_int(CBLAS_UPPER),
            ctypes.c_int(m), ctypes.c_int(n), ctypes.c_double(1.0),
            Af.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(n),
            Bf.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(m),
            ctypes.c_double(0.0),
            C.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(m),
        )
        # Reference: only the upper triangle of A is read by DSYMM
        Aup = np.triu(A) + np.triu(A, 1).T
        ref = B @ Aup
        rel = np.abs(C - ref).max() / (np.abs(ref).max() or 1.0)
        worst = max(worst, rel)
        verdict = "ok" if rel < 1e-12 else "*** WRONG ***"
        print(f"  m={m:<4} n={n:<6} max rel err = {rel:.3e}  {verdict}")
    print(f"worst   : {worst:.3e}")
    return 1 if worst > 1e-12 else 0


if __name__ == "__main__":
    sys.exit(main())