From e4bee6b8fe54a0c8020b78600d76fae76058bb0d Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Sat, 10 Oct 2026 11:31:55 +0200 Subject: [PATCH 1/3] madmatrix GPU: running flavour couplings, no good helicity, stream order, unset selections Four GPU-only bugs of the template (the CPU and SIMD backends were right), found by the review that followed the HIP cut crash (#254): - Running (alpha_s-dependent) flavour couplings, nDPF > 0 (MSSM SUSY-QCD: squark-quark-gluino GC_106/GC_110, e.g. every subprocess of p p > go go): the gather into the per-thread dpf_value and the vertex routines reading it used CD_ACCESS, which on the device adds the offset of the event in the all-events coupling buffer: every thread past the first event page read and wrote outside dpf_value. The vertex calls on flvCOUPs_dep are now instantiated with CDPF_ACCESS: CD_ACCESS on the CPU/SIMD backends (unchanged), the trivial one-event access on the GPU, where dpf_value holds this event's values. The generator refuses a vertex mixing such a coupling with another running one (one access type per routine). - No good helicity (a subprocess whose |M|^2 vanishes identically, e.g. FCNC with zero Wilson coefficients): the helicity loop was launched with gridDim.y = 0, an invalid configuration, and checkGpu aborted. Skipped now, |M|^2 = 0 as on the CPU. - The selected helicity and colour were left unset for a vanishing or NaN |M|^2 (umami then returned stale pool memory); now 0 (none, -1 out of umami) unless one is chosen, as the zero-initialised CPU buffers give. - Stream order: in async (umami) mode the memsets and the helicity, output and colour kernels went to the default stream while everything else ran on the stream umami is given. Correct only because madspace never passes its stream (a separate madspace bug) and its streams are blocking; with a non-blocking stream (a torch side stream) it races. All of sigmaKin now runs on the given stream in async mode, and umami waits for its input copies before the helicity filtering (default stream) of the first call. The memset sizes are computed in size_t. GPU CI: new job madmatrix (madmatrix_gpu.sh, madmatrix_checks.py) comparing the GPU backend with the scalar CPU backend through umami, each check in its own process so that an abort fails only that check: - mssm_gpu_vs_cpu: p p > go go event by event with random alpha_s and flavour; - zero_me: b b~ > ta+ ta- / z a with ymtau = 0 (no good helicity); - nonblocking_stream: g g > t t~ g on a non-blocking stream vs the default one. The umami ctypes harness moves to umami_harness.py, shared with interference_checks.py; its CPU buffers now keep the host arrays alive (they were freed before the call when passed as temporaries). --- .github/gpu_ci/README.md | 1 + .github/gpu_ci/interference_checks.py | 89 +----- .github/gpu_ci/madmatrix_checks.py | 257 ++++++++++++++++++ .github/gpu_ci/madmatrix_gpu.sh | 46 ++++ .github/gpu_ci/umami_harness.py | 183 +++++++++++++ .github/workflows/gpu_runner_ci.yml | 50 +++- .../madmatrix/backend/cpu/SigmaKin.cc | 5 +- .../madmatrix/backend/gpu/GpuAbstraction.h | 2 + .../madmatrix/backend/gpu/SigmaKin.cc | 63 +++-- .../madmatrix/backend/gpu/umami.cc | 4 + .../madmatrix/backend/simd/SigmaKin.cc | 5 +- madmatrix/model_handling.py | 13 +- 12 files changed, 598 insertions(+), 120 deletions(-) create mode 100644 .github/gpu_ci/madmatrix_checks.py create mode 100755 .github/gpu_ci/madmatrix_gpu.sh create mode 100644 .github/gpu_ci/umami_harness.py diff --git a/.github/gpu_ci/README.md b/.github/gpu_ci/README.md index 91857ded4..fb0ab447d 100644 --- a/.github/gpu_ci/README.md +++ b/.github/gpu_ci/README.md @@ -24,6 +24,7 @@ allocation per workflow run, and every CI job of that run executes in it: | `pp_jj` | cluster | the same with `generate p p > j j` (`PROCESS`/`TAG` of `pp_ttx_mg7.sh`): its jet cuts act on several objects, which `p p > t t~` has none of | | `madspace_tests` | cluster | the madspace unit tests, `madspace/tests`, with `MADSPACE_TEST_DEVICE=cuda` or `hip` (`madspace_tests_gpu.sh`): the tests put their torch tensors on the GPU, so madspace runs them on its GPU runtime. Host-side tests (marker `cpu_only`) are skipped. pytest and scipy are installed in the cache when the modules lack them, and so is a CUDA torch wheel when the PyTorch module is a CPU build (manneback). The log and `junit.xml` are uploaded as an artifact | | `interference` | cluster | squared split orders on the GPU backend of madmatrix (`interference_gpu.sh`, `interference_checks.py`): standalone `u u~ > u u~` constraints, the helicity choice of an interference matrix element through umami (exact and summed), and `p p > u u~ QCD^2==2` with mg7 against its CPU reference. Input `interference`, on by default and on push | +| `madmatrix` | cluster | the GPU backend of madmatrix against its scalar CPU backend through umami (`madmatrix_gpu.sh`, `madmatrix_checks.py`): MSSM `p p > go go` event by event (running flavour couplings), a subprocess with no good helicity (`b b~ > ta+ ta- / z a`, ymtau = 0), and a non-blocking stream handed to umami (`g g > t t~ g`). Each check runs in its own process | | `stop_runner` | cluster | clean up, then create the stop file: `runner_batch.sh` stops the runner and the allocation ends | | `watchdog` | vmadgraph | from the start of the runner to the end of the allocation: if the allocation ends before `stop_runner` completed, cancel the run (`ceci-gate watch-runner`, `remote_watch_runner.sh`) | diff --git a/.github/gpu_ci/interference_checks.py b/.github/gpu_ci/interference_checks.py index 3656ed907..28f4ce338 100755 --- a/.github/gpu_ci/interference_checks.py +++ b/.github/gpu_ci/interference_checks.py @@ -30,7 +30,6 @@ import argparse import ctypes -import ctypes.util import glob import gzip import json @@ -42,7 +41,7 @@ import sys import time -GPU_BACKENDS = ('cuda', 'hip') +from umami_harness import GPU_BACKENDS, Memory, rambo # u u~ > u u~ QED^2==2 at the check_sa RAMBO point, from the Fortran split-order # standalone (madmatrix CPU: -5.5828746494657258e-02) @@ -136,92 +135,6 @@ def check_standalone(repo, backend): # 2. helicity choice through umami # --------------------------------------------------------------------------- -class Memory: - """Buffers for umami: host memory for a CPU library, device memory for a GPU one - (through the CUDA/HIP runtime the library itself is linked with).""" - - def __init__(self, backend, library): - self.gpu = backend in GPU_BACKENDS - if not self.gpu: - return - name, prefix = ('cudart', 'cuda') if backend == 'cuda' else ('amdhip64', 'hip') - path = None - ldd = subprocess.run(['ldd', library], capture_output=True, text=True).stdout - for line in ldd.splitlines(): - if 'lib%s' % name in line and '=>' in line: - path = line.split('=>')[1].split()[0] - self.rt = ctypes.CDLL(path or ctypes.util.find_library(name) or 'lib%s.so' % name) - self.malloc = getattr(self.rt, prefix + 'Malloc') - self.malloc.argtypes = [ctypes.POINTER(ctypes.c_void_p), ctypes.c_size_t] - self.memcpy = getattr(self.rt, prefix + 'Memcpy') - self.memcpy.argtypes = [ctypes.c_void_p, ctypes.c_void_p, ctypes.c_size_t, ctypes.c_int] - self.free = getattr(self.rt, prefix + 'Free') - self.free.argtypes = [ctypes.c_void_p] - self.sync = getattr(self.rt, prefix + 'DeviceSynchronize') - self.pointers = [] - - def buffer(self, host): - """A buffer for umami holding (a copy of) the ctypes array host.""" - if not self.gpu: - return ctypes.addressof(host) - pointer = ctypes.c_void_p() - assert self.malloc(ctypes.byref(pointer), ctypes.sizeof(host)) == 0, 'device malloc' - # cudaMemcpyHostToDevice == hipMemcpyHostToDevice == 1 - assert self.memcpy(pointer, ctypes.addressof(host), ctypes.sizeof(host), 1) == 0 - self.pointers.append(pointer) - return pointer.value - - def fetch(self, pointer, host): - """Copy a umami output back into the ctypes array host.""" - if not self.gpu: - return - assert self.sync() == 0, 'device synchronize' - # cudaMemcpyDeviceToHost == hipMemcpyDeviceToHost == 2 - assert self.memcpy(ctypes.addressof(host), pointer, ctypes.sizeof(host), 2) == 0 - - def release(self): - if self.gpu: - for pointer in self.pointers: - self.free(pointer) - self.pointers = [] - - -def rambo(masses_out, roots, rng): - """One phase-space point: two massless beams along z, then RAMBO with masses.""" - n = len(masses_out) - q = [] - for _ in range(n): - c, f = 2 * rng.random() - 1, 2 * math.pi * rng.random() - e = -math.log(rng.random() * rng.random()) - s = math.sqrt(1 - c * c) - q.append([e, e * s * math.cos(f), e * s * math.sin(f), e * c]) - big_q = [sum(qi[k] for qi in q) for k in range(4)] - mass = math.sqrt(big_q[0] ** 2 - sum(big_q[k] ** 2 for k in (1, 2, 3))) - b = [-big_q[k] / mass for k in (1, 2, 3)] - g, x = big_q[0] / mass, roots / mass - a = 1 / (1 + g) - p = [] - for qi in q: - bq = sum(b[k] * qi[k + 1] for k in range(3)) - p.append([x * (g * qi[0] + bq)] + - [x * (qi[k + 1] + b[k] * qi[0] + a * bq * b[k]) for k in range(3)]) - # rescale the three-momenta so that the energies add up to roots with the masses - xi = 1. - for _ in range(100): - energies = [math.sqrt((xi * pi[0]) ** 2 + m * m) for pi, m in zip(p, masses_out)] - f = sum(energies) - roots - df = sum(xi * pi[0] ** 2 / e for pi, e in zip(p, energies)) - xi -= f / df - if abs(f) < 1e-12 * roots: - break - out = [] - for pi, m in zip(p, masses_out): - vec = [xi * v for v in pi[1:]] - out.append([math.sqrt(sum(v * v for v in vec) + m * m)] + vec) - half = roots / 2 - return [[half, 0., 0., half], [half, 0., 0., -half]] + out - - def check_helicity_choice(repo, backend, npoints=16, ngrid=4000): section('helicity choice of an interference |M|^2 through umami (%s)' % backend) proc_dir = standalone(repo, backend, 'uuxttxg_qed2_d', 'u u~ > t t~ g QED^2==2', 'd') diff --git a/.github/gpu_ci/madmatrix_checks.py b/.github/gpu_ci/madmatrix_checks.py new file mode 100644 index 000000000..8eec47646 --- /dev/null +++ b/.github/gpu_ci/madmatrix_checks.py @@ -0,0 +1,257 @@ +#!/usr/bin/env python3 +"""The GPU backend of madmatrix against its scalar CPU backend, through umami. + +Run by madmatrix_gpu.sh (job madmatrix of gpu_runner_ci.yml) in its run directory. + +1. mssm_gpu_vs_cpu: MSSM ``p p > go go`` (both subprocesses), whose squark-quark-gluino + couplings are running flavour couplings (nDPF > 0, gathered per event into + dpf_value), on random points with a random alpha_s and flavour per event: the GPU + |M|^2 equal the CPU ones event by event (the GPU gather used to address its one-event + dpf_value with the offset of the event in the all-events coupling buffer, out of + bounds past the first event page), and scale as alpha_s^2. +2. zero_me: ``b b~ > ta+ ta- / z a`` (s-channel Higgs only) with the tau Yukawa coupling + set to 0 in the param card, so that no helicity survives the helicity filtering: the + GPU library returns |M|^2 = 0 and no helicity (-1), and for the colour either none + (-1) or the only colour flow (0, what a -ffast-math build makes of 0/0, as the CPU + one), never stale memory; it used to abort on the invalid launch configuration + (gridDim.y = 0) of the helicity loop. +3. nonblocking_stream: ``g g > t t~ g`` with umami handed a non-blocking stream (as a + torch side stream, which does not synchronise with the default stream) gives the same + |M|^2 and helicities as on the default stream, call after call: every memset, kernel + and copy of the matrix element must go to the stream umami is given. + +Each check runs in a process of its own (--only NAME): an abort of the GPU runtime +fails that check, not the others. Writes summary.txt (key=value lines) in the current +directory; exit status 1 if a check fails. +""" + +import argparse +import glob +import json +import os +import random +import re +import subprocess +import sys +import time + +from umami_harness import (GPU_BACKENDS, OUT_COLOR_INDEX, OUT_HELICITY_INDEX, + OUT_MATRIX_ELEMENT, Umami, rambo) + +CPU_BACKEND = 'scalar' +CHECKS = ('mssm_gpu_vs_cpu', 'zero_me', 'nonblocking_stream') + +results = {} + + +def section(title): + print('\n=================== %s (%s)' % (title, time.strftime('%H:%M:%S')), flush=True) + + +def check(name, ok, detail): + print('%s %s: %s' % ('PASS' if ok else 'FAIL', name, detail), flush=True) + results[name] = (ok, detail) + + +def run(cmd, cwd=None, log=None): + with open(log or os.devnull, 'w') as out: + return subprocess.call(cmd, cwd=cwd, stdout=out, stderr=subprocess.STDOUT) + + +def output(repo, name, model, process, backends): + """output standalone + make for these backends in every subprocess: the output + directory and its subprocess names.""" + with open(name + '.mg5', 'w') as f: + f.write('import model %s\ngenerate %s\noutput standalone %s -f\n' % (model, process, name)) + if run([sys.executable, os.path.join(repo, 'bin', 'madgraph'), name + '.mg5'], + log=name + '.log') != 0: + raise RuntimeError('madgraph failed, see %s.log' % name) + subprocesses = sorted(os.path.basename(d) + for d in glob.glob(os.path.join(name, 'SubProcesses', 'P*'))) + if not subprocesses: + raise RuntimeError('no subprocess directory for %s' % process) + jobs = os.environ.get('SLURM_CPUS_PER_TASK', '4') + for sub in subprocesses: + for backend in backends: + log = '%s_%s_%s_make.log' % (name, sub, backend) + if run(['make', '-j', jobs, 'BACKEND=' + backend, 'FPTYPE=d'], + cwd=os.path.join(name, 'SubProcesses', sub), log=log) != 0: + raise RuntimeError('make failed, see %s' % log) + return subprocesses + + +def library(name, sub, backend): + return os.path.join(name, 'lib', 'libmadmatrix_%s_%s.so' % (sub, backend)) + + +def card(name): + return os.path.join(name, 'Cards', 'param_card.dat') + + +# --------------------------------------------------------------------------- +# 1. running flavour couplings (MSSM p p > go go): GPU against CPU +# --------------------------------------------------------------------------- + +def check_mssm_gpu_vs_cpu(repo, backend, npoints=1024): + section('MSSM p p > go go, running flavour couplings: %s against %s' % (backend, CPU_BACKEND)) + name = 'mssm_gogo' + subprocesses = output(repo, name, 'MSSM_SLHA2', 'p p > go go', (CPU_BACKEND, backend)) + rng = random.Random(11) + problems, details = [], [] + for sub in subprocesses: + data = open(os.path.join(name, 'SubProcesses', sub, 'ProcessData.h')).read() + ndpf = int(re.search(r'constexpr int nDPF = (\d+);', data).group(1)) + nflavor = int(re.search(r'constexpr int nmaxflavor = (\d+);', data).group(1)) + if ndpf == 0: + problems.append('%s: nDPF = 0, no running flavour coupling to check' % sub) + cpu = Umami(CPU_BACKEND, library(name, sub, CPU_BACKEND), card(name)) + gpu = Umami(backend, library(name, sub, backend), card(name)) + points = [rambo(gpu.masses[2:], 3000., rng, gpu.masses[:2]) for _ in range(npoints)] + # the first two events: one point at two values of alpha_s + points[1] = points[0] + alpha_s = [0.10, 0.12] + [rng.uniform(0.09, 0.14) for _ in range(npoints - 2)] + flavor = [rng.randrange(nflavor) for _ in range(npoints)] + (me_cpu,) = cpu(points, alpha_s=alpha_s, flavor=flavor) + (me_gpu,) = gpu(points, alpha_s=alpha_s, flavor=flavor) + worst, iworst = 0., -1 + for ievt, (a, b) in enumerate(zip(me_cpu, me_gpu)): + rel = abs(a - b) / abs(a) if a != 0 else (0. if b == 0 else float('inf')) + if not rel <= worst: # also catches nan + worst, iworst = rel, ievt + if not worst < 1e-9: + problems.append('%s: event %d %s %.10e vs %s %.10e (rel %.1e)' + % (sub, iworst, backend, me_gpu[iworst], CPU_BACKEND, + me_cpu[iworst], worst)) + if not any(v != 0 for v in me_cpu): + problems.append('%s: all |M|^2 vanish' % sub) + scaling = me_gpu[1] / me_gpu[0] / (0.12 / 0.10) ** 2 - 1 if me_gpu[0] else float('inf') + if not abs(scaling) < 1e-9: + problems.append('%s: |M|^2 not ~ alpha_s^2 (ratio - 1 = %.1e)' % (sub, scaling)) + details.append('%s nDPF %d, %d flavours, max rel diff %.1e' % (sub, ndpf, nflavor, worst)) + print(' ' + details[-1]) + cpu.free() + gpu.free() + check('mssm_gpu_vs_cpu', not problems, '; '.join(problems[:3]) or '; '.join(details)) + + +# --------------------------------------------------------------------------- +# 2. a subprocess with no good helicity +# --------------------------------------------------------------------------- + +def check_zero_me(repo, backend, npoints=64): + section('no good helicity: b b~ > ta+ ta- / z a with ymtau = 0 (%s)' % backend) + name = 'zero_me' + (sub,) = output(repo, name, 'sm', 'b b~ > ta+ ta- / z a', (CPU_BACKEND, backend)) + text = open(card(name)).read() + zero_card, count = re.subn(r'(?mi)^(\s*15\s+)\S+(\s*#\s*ymtau)', r'\g<1>0.000000e+00\2', text) + if count != 1: + raise RuntimeError('no ymtau in %s' % card(name)) + zero_path = os.path.join(name, 'Cards', 'param_card_ymtau0.dat') + open(zero_path, 'w').write(zero_card) + rng = random.Random(5) + # the helicity filtering runs once per library (and process): the GPU library only + # ever sees the card with ymtau = 0, the CPU library the default one (|M|^2 != 0) + gpu = Umami(backend, library(name, sub, backend), zero_path) + cpu = Umami(CPU_BACKEND, library(name, sub, CPU_BACKEND), card(name)) + points = [rambo(gpu.masses[2:], 500., rng, gpu.masses[:2]) for _ in range(npoints)] + rnd = [rng.random() for _ in range(npoints)] + outputs = (OUT_MATRIX_ELEMENT, OUT_HELICITY_INDEX, OUT_COLOR_INDEX) + me, hel, col = gpu(points, random_helicity=rnd, random_color=rnd, outputs=outputs) + (me_default,) = cpu(points) + problems = [] + if any(v != 0 for v in me): + problems.append('|M|^2 not 0: %s' % me[:4]) + if set(hel) != {-1} or not set(col) <= {-1, 0}: + problems.append('helicities %s, colours %s (expected -1, and -1 or 0)' + % (sorted(set(hel))[:4], sorted(set(col))[:4])) + if not all(v > 0 for v in me_default): + problems.append('|M|^2 with the default card not > 0: %s' % me_default[:4]) + gpu.free() + cpu.free() + check('zero_me', not problems, '; '.join(problems) or + '%d events: |M|^2 = 0, helicity -1, colour %s (default card: %.3e)' + % (npoints, sorted(set(col)), me_default[0])) + + +# --------------------------------------------------------------------------- +# 3. a non-blocking stream +# --------------------------------------------------------------------------- + +def check_nonblocking_stream(repo, backend, npoints=16384, ncalls=4): + section('g g > t t~ g on a non-blocking stream (%s)' % backend) + name = 'ggttxg' + (sub,) = output(repo, name, 'sm', 'g g > t t~ g', (backend,)) + rng = random.Random(3) + gpu = Umami(backend, library(name, sub, backend), card(name)) + points = [rambo(gpu.masses[2:], 1000., rng) for _ in range(npoints)] + rnd = [rng.random() for _ in range(npoints)] + outputs = (OUT_MATRIX_ELEMENT, OUT_HELICITY_INDEX) + stream = gpu.memory.nonblocking_stream() + # the first call on the non-blocking stream: it also runs the helicity filtering + calls = [gpu(points, random_helicity=rnd, outputs=outputs, stream=stream)] + reference = gpu(points, random_helicity=rnd, outputs=outputs) # default stream + calls += [gpu(points, random_helicity=rnd, outputs=outputs, stream=stream) + for _ in range(ncalls - 1)] + gpu.memory.stream_destroy(stream) + gpu.free() + me_ref, hel_ref = reference + problems = [] + for icall, (me, hel) in enumerate(calls): + bad_me = sum(1 for a, b in zip(me, me_ref) if not abs(a - b) <= 1e-12 * abs(b)) + bad_hel = sum(1 for a, b in zip(hel, hel_ref) if a != b) + if bad_me or bad_hel: + problems.append('call %d: %d |M|^2 and %d helicities differ' % (icall, bad_me, bad_hel)) + if not all(v > 0 for v in me_ref) or min(hel_ref) < 0: + problems.append('default stream: |M|^2 <= 0 or no helicity') + check('nonblocking_stream', not problems, '; '.join(problems[:3]) or + '%d calls of %d events identical to the default stream' % (ncalls, npoints)) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.split('\n')[0]) + parser.add_argument('--repo', required=True) + parser.add_argument('--backend', required=True) + parser.add_argument('--only', choices=CHECKS, + help='run this check here and write its result to NAME.json') + args = parser.parse_args() + if args.only: + step = globals()['check_' + args.only] + try: + step(args.repo, args.backend) + except Exception as error: + check(args.only, False, 'error: %s' % error) + json.dump(results, open(args.only + '.json', 'w')) + return 0 + if args.backend not in GPU_BACKENDS: + parser.error('the checks compare a GPU backend (%s) with the CPU one' + % ', '.join(GPU_BACKENDS)) + start = time.time() + summary, failures = {}, [] + for name in CHECKS: + if os.path.exists(name + '.json'): + os.remove(name + '.json') + status = subprocess.call([sys.executable, os.path.abspath(__file__), '--repo', args.repo, + '--backend', args.backend, '--only', name]) + try: + outcome = json.load(open(name + '.json')) + except (OSError, ValueError): + outcome = {name: (False, 'crashed (exit status %d)' % status)} + print('FAIL %s: %s' % (name, outcome[name][1]), flush=True) + for key, (ok, detail) in outcome.items(): + summary[key] = ('ok' if ok else 'FAILED') + ' (%s)' % detail + if not ok: + failures.append(key) + summary['node'] = os.uname().nodename + summary['gpu'] = os.environ.get('GPU_NAME', 'none') + summary['backend'] = args.backend + summary['walltime'] = '%ds' % (time.time() - start) + with open('summary.txt', 'w') as f: + for key, value in summary.items(): + f.write('%s=%s\n' % (key, value)) + print('\n%d check(s) failed: %s' % (len(failures), ', '.join(failures)) if failures + else '\nall checks passed') + return 1 if failures else 0 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/.github/gpu_ci/madmatrix_gpu.sh b/.github/gpu_ci/madmatrix_gpu.sh new file mode 100755 index 000000000..cc6d297d5 --- /dev/null +++ b/.github/gpu_ci/madmatrix_gpu.sh @@ -0,0 +1,46 @@ +#!/bin/bash -l +# Runs on the self-hosted runner (GPU node): job madmatrix of gpu_runner_ci.yml. +# The GPU backend of madmatrix against its CPU backend, through umami: see +# madmatrix_checks.py for what is checked. +# Environment: as pp_ttx_mg7.sh (BACKEND, MODULES, GPU_ARCH, VENV, WORKDIR); no madspace +# and no PDF set needed (standalone outputs only). +set -eo pipefail +HERE=$(cd "$(dirname "$0")" && pwd) +REPO=$(cd "$HERE/../.." && pwd) +section() { echo; echo "=================== $* ($(date +%T))"; } + +section "Environment" +if [ -n "$MODULES" ]; then source "$HERE/load_modules.sh"; fi +source "$VENV/bin/activate" +python3 --version +case $BACKEND in + cuda) nvidia-smi + GPU_NAME=$(nvidia-smi --query-gpu=name --format=csv,noheader | sed -n 1p) ;; + hip) rocm-smi || true + GPU_NAME=$(rocminfo 2> /dev/null | awk -F': *' '/Marketing Name/ {name=$2} /Device Type: *GPU/ {print name; exit}' || true) ;; + *) echo "::error::the madmatrix checks need a GPU backend, not '$BACKEND'"; exit 2 ;; +esac +export GPU_NAME +if [ -n "$GPU_ARCH" ]; then + export MADGRAPH_CUDA_ARCHITECTURE=$GPU_ARCH MADGRAPH_HIP_ARCHITECTURE=$GPU_ARCH +fi + +# default configuration, as .github/actions/checkout_mg5 does for the other CI jobs +cp "$REPO/input/.mg7_configuration_default.txt" "$REPO/input/mg7_configuration.txt" +cp "$REPO/Template/LO/Source/.make_opts" "$REPO/Template/LO/Source/make_opts" +cat >> "$REPO/input/mg7_configuration.txt" << EOF +auto_update = 0 +automatic_html_opening = False +notification_center = False +nb_core = ${SLURM_CPUS_PER_TASK:-4} +EOF + +rm -rf "$WORKDIR" +mkdir -p "$WORKDIR" +cd "$WORKDIR" +STATUS=0 +python3 "$HERE/madmatrix_checks.py" --repo "$REPO" --backend "$BACKEND" 2>&1 | tee checks.log || STATUS=1 + +section "Summary" +cat summary.txt 2> /dev/null || echo "no summary.txt" +exit $STATUS diff --git a/.github/gpu_ci/umami_harness.py b/.github/gpu_ci/umami_harness.py new file mode 100644 index 000000000..2b772029a --- /dev/null +++ b/.github/gpu_ci/umami_harness.py @@ -0,0 +1,183 @@ +"""Calling a madmatrix process library through umami from Python (ctypes), on the CPU or +on the GPU, for the checks of the GPU CI (interference_checks.py, madmatrix_checks.py). +""" + +import ctypes +import ctypes.util +import math +import subprocess + +GPU_BACKENDS = ('cuda', 'hip') + +# umami.h +META_PARTICLE_COUNT, META_MASSES = 1, 5 +IN_MOMENTA, IN_ALPHA_S, IN_FLAVOR_INDEX, IN_RANDOM_COLOR, IN_RANDOM_HELICITY = 0, 1, 2, 3, 4 +OUT_MATRIX_ELEMENT, OUT_COLOR_INDEX, OUT_HELICITY_INDEX, OUT_GPU_STREAM = 0, 2, 3, 5 + + +class Memory: + """Buffers for umami: host memory for a CPU library, device memory for a GPU one + (through the CUDA/HIP runtime the library itself is linked with).""" + + def __init__(self, backend, library): + self.gpu = backend in GPU_BACKENDS + self.hosts = [] # the host arrays umami reads on a CPU, alive until release + if not self.gpu: + return + name, prefix = ('cudart', 'cuda') if backend == 'cuda' else ('amdhip64', 'hip') + path = None + ldd = subprocess.run(['ldd', library], capture_output=True, text=True).stdout + for line in ldd.splitlines(): + if 'lib%s' % name in line and '=>' in line: + path = line.split('=>')[1].split()[0] + self.rt = ctypes.CDLL(path or ctypes.util.find_library(name) or 'lib%s.so' % name) + self.malloc = getattr(self.rt, prefix + 'Malloc') + self.malloc.argtypes = [ctypes.POINTER(ctypes.c_void_p), ctypes.c_size_t] + self.memcpy = getattr(self.rt, prefix + 'Memcpy') + self.memcpy.argtypes = [ctypes.c_void_p, ctypes.c_void_p, ctypes.c_size_t, ctypes.c_int] + self.free = getattr(self.rt, prefix + 'Free') + self.free.argtypes = [ctypes.c_void_p] + self.sync = getattr(self.rt, prefix + 'DeviceSynchronize') + self.stream_create = getattr(self.rt, prefix + 'StreamCreateWithFlags') + self.stream_create.argtypes = [ctypes.POINTER(ctypes.c_void_p), ctypes.c_uint] + self.stream_destroy = getattr(self.rt, prefix + 'StreamDestroy') + self.stream_destroy.argtypes = [ctypes.c_void_p] + self.pointers = [] + + def buffer(self, host): + """A buffer for umami holding (a copy of) the ctypes array host.""" + if not self.gpu: + self.hosts.append(host) + return ctypes.addressof(host) + pointer = ctypes.c_void_p() + assert self.malloc(ctypes.byref(pointer), ctypes.sizeof(host)) == 0, 'device malloc' + # cudaMemcpyHostToDevice == hipMemcpyHostToDevice == 1 + assert self.memcpy(pointer, ctypes.addressof(host), ctypes.sizeof(host), 1) == 0 + self.pointers.append(pointer) + return pointer.value + + def fetch(self, pointer, host): + """Copy a umami output back into the ctypes array host.""" + if not self.gpu: + return + assert self.sync() == 0, 'device synchronize' + # cudaMemcpyDeviceToHost == hipMemcpyDeviceToHost == 2 + assert self.memcpy(ctypes.addressof(host), pointer, ctypes.sizeof(host), 2) == 0 + + def release(self): + self.hosts = [] + if self.gpu: + for pointer in self.pointers: + self.free(pointer) + self.pointers = [] + + def nonblocking_stream(self): + """A stream that does not synchronise with the default stream (as a torch side + stream): cudaStreamNonBlocking == hipStreamNonBlocking == 1.""" + stream = ctypes.c_void_p() + assert self.stream_create(ctypes.byref(stream), 1) == 0, 'stream create' + return stream.value + + +class Umami: + """One umami instance of a process library.""" + + def __init__(self, backend, library, param_card): + self.lib = ctypes.CDLL(library) + self.memory = Memory(backend, library) + self.handle = ctypes.c_void_p() + status = self.lib.umami_initialize(ctypes.byref(self.handle), param_card.encode()) + assert status == 0, 'umami_initialize status %d' % status + npar = ctypes.c_int() + self.lib.umami_get_meta(META_PARTICLE_COUNT, ctypes.byref(npar)) + self.npar = npar.value + masses = (ctypes.c_double * self.npar)() + assert self.lib.umami_get_meta(META_MASSES, masses) == 0, 'UMAMI_META_MASSES' + self.masses = list(masses) + self.lib.umami_matrix_element.argtypes = [ + ctypes.c_void_p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_size_t, + ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p), + ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p)] + self.lib.umami_set_parameter.argtypes = [ctypes.c_void_p, ctypes.c_char_p, + ctypes.c_double, ctypes.c_double] + + def set_parameter(self, name, value): + return self.lib.umami_set_parameter(self.handle, name.encode(), value, 0.) + + def __call__(self, points, alpha_s=None, flavor=None, random_helicity=None, + random_color=None, outputs=(OUT_MATRIX_ELEMENT,), stream=None): + """umami_matrix_element on these phase-space points (lists of npar 4-momenta) + and per-event inputs; the outputs as lists, in the order asked for.""" + count = len(points) + mom = (ctypes.c_double * (4 * self.npar * count))() + for ievt, point in enumerate(points): # momenta[count * (npar * mu + ipart) + ievt] + for ipart in range(self.npar): + for mu in range(4): + mom[count * (self.npar * mu + ipart) + ievt] = point[ipart][mu] + in_keys, inputs = [IN_MOMENTA], [self.memory.buffer(mom)] + for key, values, ctype in ((IN_ALPHA_S, alpha_s, ctypes.c_double), + (IN_FLAVOR_INDEX, flavor, ctypes.c_int), + (IN_RANDOM_HELICITY, random_helicity, ctypes.c_double), + (IN_RANDOM_COLOR, random_color, ctypes.c_double)): + if values is not None: + in_keys.append(key) + inputs.append(self.memory.buffer((ctype * count)(*values))) + hosts = [((ctypes.c_double if key == OUT_MATRIX_ELEMENT else ctypes.c_int) * count)() + for key in outputs] + out_keys = list(outputs) + out_pointers = [self.memory.buffer(host) for host in hosts] + if stream is not None: + out_keys.append(OUT_GPU_STREAM) + out_pointers.append(stream) + status = self.lib.umami_matrix_element( + self.handle, count, count, 0, + len(in_keys), (ctypes.c_int * len(in_keys))(*in_keys), + (ctypes.c_void_p * len(inputs))(*inputs), + len(out_keys), (ctypes.c_int * len(out_keys))(*out_keys), + (ctypes.c_void_p * len(out_pointers))(*out_pointers)) + assert status == 0, 'umami_matrix_element status %d' % status + for pointer, host in zip(out_pointers, hosts): + self.memory.fetch(pointer, host) + self.memory.release() + return [list(host) for host in hosts] + + def free(self): + self.lib.umami_free(self.handle) + + +def rambo(masses_out, roots, rng, masses_in=(0., 0.)): + """One phase-space point: two beams along z (masses masses_in) in their centre-of-mass + frame, then RAMBO with the masses masses_out.""" + n = len(masses_out) + q = [] + for _ in range(n): + c, f = 2 * rng.random() - 1, 2 * math.pi * rng.random() + e = -math.log(rng.random() * rng.random()) + s = math.sqrt(1 - c * c) + q.append([e, e * s * math.cos(f), e * s * math.sin(f), e * c]) + big_q = [sum(qi[k] for qi in q) for k in range(4)] + mass = math.sqrt(big_q[0] ** 2 - sum(big_q[k] ** 2 for k in (1, 2, 3))) + b = [-big_q[k] / mass for k in (1, 2, 3)] + g, x = big_q[0] / mass, roots / mass + a = 1 / (1 + g) + p = [] + for qi in q: + bq = sum(b[k] * qi[k + 1] for k in range(3)) + p.append([x * (g * qi[0] + bq)] + + [x * (qi[k + 1] + b[k] * qi[0] + a * bq * b[k]) for k in range(3)]) + # rescale the three-momenta so that the energies add up to roots with the masses + xi = 1. + for _ in range(100): + energies = [math.sqrt((xi * pi[0]) ** 2 + m * m) for pi, m in zip(p, masses_out)] + f = sum(energies) - roots + df = sum(xi * pi[0] ** 2 / e for pi, e in zip(p, energies)) + xi -= f / df + if abs(f) < 1e-12 * roots: + break + out = [] + for pi, m in zip(p, masses_out): + vec = [xi * v for v in pi[1:]] + out.append([math.sqrt(sum(v * v for v in vec) + m * m)] + vec) + m1, m2 = masses_in + pz = math.sqrt((roots ** 2 - (m1 + m2) ** 2) * (roots ** 2 - (m1 - m2) ** 2)) / (2 * roots) + return [[math.hypot(pz, m1), 0., 0., pz], [math.hypot(pz, m2), 0., 0., -pz]] + out diff --git a/.github/workflows/gpu_runner_ci.yml b/.github/workflows/gpu_runner_ci.yml index fb2a5490c..c848450ae 100644 --- a/.github/workflows/gpu_runner_ci.yml +++ b/.github/workflows/gpu_runner_ci.yml @@ -9,6 +9,7 @@ # pp_jj cluster runner p p > j j with mg7 on the GPU (cuts on several jets) # madspace_tests cluster runner the madspace unit tests on the GPU, .github/gpu_ci/madspace_tests_gpu.sh # interference cluster runner squared split orders (interference), .github/gpu_ci/interference_gpu.sh +# madmatrix cluster runner GPU against CPU backend of madmatrix, .github/gpu_ci/madmatrix_gpu.sh # ... cluster runner more CI jobs: same needs, if and runs-on, and add them to # stop_runner's needs # stop_runner cluster runner clean up and end the allocation (no ssh key needed) @@ -249,6 +250,53 @@ jobs: path: results/ if-no-files-found: warn + madmatrix: + # the GPU backend of madmatrix against its CPU backend through umami: running flavour + # couplings (MSSM), a subprocess with no good helicity, a non-blocking stream + needs: [start_runner, build_madspace] + if: needs.start_runner.outputs.attempt == github.run_attempt + runs-on: [self-hosted, "${{ inputs.runner_label }}"] + timeout-minutes: 90 + defaults: + run: + shell: bash -l -eo pipefail {0} + env: + BACKEND: ${{ inputs.backend }} + MODULES: ${{ inputs.modules }} + GPU_ARCH: ${{ inputs.gpu_arch }} + VENV: ${{ needs.build_madspace.outputs.venv }} + steps: + - uses: actions/checkout@v5 + - name: run directory + run: echo "WORKDIR=${GLOBALSCRATCH:-$HOME}/mg7-gpu-ci/runs/$RUN_TAG/madmatrix" >> "$GITHUB_ENV" + - name: madmatrix on the GPU + run: .github/gpu_ci/madmatrix_gpu.sh + - name: collect logs + if: always() + run: | + mkdir -p results + (cd "$WORKDIR" && tar -cf - --ignore-failed-read summary.txt checks.log *.json *.mg5 *.log 2> /dev/null) | tar -xf - -C results || true + find results -type f | sort + - name: job summary + if: always() + run: | + { + echo "### madmatrix on ${{ inputs.cluster }} ($BACKEND)" + if [ -f results/summary.txt ]; then + echo "| | |" + echo "|---|---|" + sed 's/^\([^=]*\)=\(.*\)$/| \1 | \2 |/' results/summary.txt + else + echo "No summary.txt: the job stopped early, see the log of the 'madmatrix on the GPU' step." + fi + } >> "$GITHUB_STEP_SUMMARY" + - uses: actions/upload-artifact@v4 + if: always() + with: + name: gpu-madmatrix-${{ inputs.cluster }}-${{ inputs.backend }} + path: results/ + if-no-files-found: warn + pp_jj: # p p > j j: cuts on several objects (the jets), which p p > t t~ does not have needs: [start_runner, build_madspace] @@ -355,7 +403,7 @@ jobs: stop_runner: # must list every job that runs on the self-hosted runner - needs: [start_runner, build_madspace, pp_ttx, pp_jj, madspace_tests, interference] + needs: [start_runner, build_madspace, pp_ttx, pp_jj, madspace_tests, interference, madmatrix] if: always() && needs.start_runner.result == 'success' # In an attempt that did not start the runner (a re-run of failed jobs, or of one job), a job # for the runner label would wait 24 hours, holding the concurrency group of the cluster. diff --git a/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc index be4b73191..3e1ef631f 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc @@ -197,6 +197,7 @@ namespace madmatrix using A_ACCESS = HostAccessAmplitudes; // TRIVIAL ACCESS (no kernel splitting yet): buffer for one event using CD_ACCESS = HostAccessCouplings; // non-trivial access (dependent couplings): buffer includes all events using CI_ACCESS = HostAccessCouplingsFixed; // TRIVIAL access (independent couplings): buffer for one event + using CDPF_ACCESS = CD_ACCESS; // the running flavor couplings gathered into dpf_value (one event page, as CD_ACCESS) using F_ACCESS = HostAccessIflavorVec; // non-trivial access: buffer includes all events using NUM_ACCESS = HostAccessNumerators; // non-trivial access: buffer includes all events mgDebug( 0, __FUNCTION__ ); @@ -266,8 +267,8 @@ namespace madmatrix // an ordinary value-based view over it. The flavor index is constant across a SIMD lane // (guaranteed by the phase-space integrator), so each lane gets its own running value // while sharing the same flavor selection. This is the direct analogue of Fortran's - // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CD_ACCESS so - // get_coupling_def reads dpf_value with the right per-flavor stride (CD_ACCESS::flv_stride). + // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CDPF_ACCESS + // (= CD_ACCESS here) so get_coupling_def reads dpf_value with the right per-flavor stride. constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); // cppAlign is only defined for SIMD alignas( mgOnGpu::cppAlign ) fptype dpf_value[ndpfbuf]{}; diff --git a/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h b/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h index b17ea25c4..cdb6a0f11 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h +++ b/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h @@ -29,6 +29,7 @@ #define gpuMemcpy( dstData, srcData, srcBytes, func ) checkGpu( cudaMemcpy( dstData, srcData, srcBytes, func ) ) #define gpuMemset( data, value, bytes ) checkGpu( cudaMemset( data, value, bytes ) ) +#define gpuMemsetAsync( data, value, bytes, stream ) checkGpu( cudaMemsetAsync( data, value, bytes, stream ) ) #define gpuMemcpyHostToDevice cudaMemcpyHostToDevice #define gpuMemcpyDeviceToHost cudaMemcpyDeviceToHost #define gpuMemcpyDeviceToDevice cudaMemcpyDeviceToDevice @@ -98,6 +99,7 @@ #define gpuMemcpy( dstData, srcData, srcBytes, func ) checkGpu( hipMemcpy( dstData, srcData, srcBytes, func ) ) #define gpuMemset( data, value, bytes ) checkGpu( hipMemset( data, value, bytes ) ) +#define gpuMemsetAsync( data, value, bytes, stream ) checkGpu( hipMemsetAsync( data, value, bytes, stream ) ) #define gpuMemcpyHostToDevice hipMemcpyHostToDevice #define gpuMemcpyDeviceToHost hipMemcpyDeviceToHost #define gpuMemcpyDeviceToDevice hipMemcpyDeviceToDevice diff --git a/madgraph/iolibs/template_files/madmatrix/backend/gpu/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/gpu/SigmaKin.cc index 4c1ad38fa..e1c84c8ef 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/gpu/SigmaKin.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/gpu/SigmaKin.cc @@ -183,6 +183,7 @@ namespace madmatrix using A_ACCESS = DeviceAccessAmplitudes; // TRIVIAL ACCESS (no kernel splitting yet): buffer for one event using CD_ACCESS = DeviceAccessCouplings; // non-trivial access (dependent couplings): buffer includes all events using CI_ACCESS = DeviceAccessCouplingsFixed; // TRIVIAL access (independent couplings): buffer for one event + using CDPF_ACCESS = CI_ACCESS; // TRIVIAL access: the running flavor couplings gathered into dpf_value for this event using F_ACCESS = DeviceAccessIflavorVec; // non-trivial access: buffer includes all events using NUM_ACCESS = DeviceAccessNumerators; // non-trivial access: buffer includes all events mgDebug( 0, __FUNCTION__ ); @@ -257,24 +258,27 @@ namespace madmatrix // Dependent (event-by-event, running-alphas) flavor couplings (Step 3): the per-flavor // values are NOT baked in (they run per event). Gather the current values of the - // underlying dependent couplings for this event page into an AOSOA buffer dpf_value - // (one nx2*neppC SIMD record per (coupling,flavor) slot, matching CD_ACCESS), then build - // an ordinary value-based view over it. The flavor index is constant across a SIMD lane - // (guaranteed by the phase-space integrator), so each lane gets its own running value - // while sharing the same flavor selection. This is the direct analogue of Fortran's - // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CD_ACCESS so - // get_coupling_def reads dpf_value with the right per-flavor stride (CD_ACCESS::flv_stride). - constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); + // underlying dependent couplings for this event into dpf_value, one complex per + // (coupling,flavor) slot, then build an ordinary value-based view over it. This is the + // direct analogue of Fortran's FLV_xx%VAL(k)%P => GC_yyy(J). A GPU thread is one event, + // so dpf_value is a one-event buffer, read with the trivial CDPF_ACCESS (= CI_ACCESS) by + // the vertex routines; CD_ACCESS would add the offset of this event in the all-events + // coupling buffer, outside dpf_value for every event past the first page. + constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CDPF_ACCESS::flv_stride : 1 ); fptype dpf_value[ndpfbuf]{}; for( int idpf = 0; idpf < nDPF; idpf++ ) for( int imf = 0; imf < nMF; imf++ ) { const int idc = cDPF_idcoup[idpf * nMF + imf]; if( idc >= 0 ) - CD_ACCESS::kernelAccess( dpf_value + ( idpf * nMF + imf ) * CD_ACCESS::flv_stride ) = - CD_ACCESS::kernelAccessConst( COUPs[idc] ); + { + const cxtype coup = CD_ACCESS::kernelAccessConst( COUPs[idc] ); + fptype* slot = dpf_value + ( idpf * nMF + imf ) * CDPF_ACCESS::flv_stride; + slot[0] = cxreal( coup ); + slot[1] = cximag( coup ); + } } - FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; + FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; // Reset color flows (reset jamp_sv) at the beginning of a new event or event page for( int i = 0; i < njampso; i++ ) { jamp_sv[i] = cxzero_sv(); } @@ -489,6 +493,9 @@ namespace madmatrix const bool sampleSignedHelicity ) // input: see sigmaKin { const int ievt = blockDim.x * blockIdx.x + threadIdx.x; // index of event (thread) + // no helicity (0) unless one is chosen below: a vanishing or NaN |M|^2 chooses none, and + // allselhel is not cleared beforehand (the CPU buffers are zero-initialised) + allselhel[ievt] = 0; if constexpr( nampso > 1 ) { if( sampleSignedHelicity ) @@ -501,7 +508,7 @@ namespace madmatrix for( int ighel = 0; ighel < dcNGoodHel; ighel++ ) absTotal += fabs( ghelAllMEs[ighel * nevt + ievt] ); allMEs[ievt] = 0; - if( absTotal == 0 ) return; // selhel left alone, as in the positive case + if( absTotal == 0 ) return; // no helicity, as in the positive case fptype absSum = 0; int ighelLast = -1; // last non-vanishing helicity, in case rndhel is not below 1 for( int ighel = 0; ighel < dcNGoodHel; ighel++ ) @@ -612,6 +619,7 @@ namespace madmatrix // NB (see #877): in the array icolamp, the input index uses C indexing (iconfig -1) if( mgOnGpu::icolamp[iconfig - 1][icolC] ) targetamp[icolC] += jamp2_sv[icolC]; } + allselcol[ievt] = 0; // no color unless one is chosen below (vanishing or NaN jamp2 sum) for( int icolC = 0; icolC < ncolor_flow; icolC++ ) { if( allrndcol[ievt] < ( targetamp[icolC] / targetamp[ncolor_flow - 1] ) ) @@ -667,30 +675,39 @@ namespace madmatrix } // === PART 0 - INITIALISATION (before calculate_jamps) === + // In async mode (umami) all the work goes to the caller's stream ghelStreams[0], in stream + // order with what the caller does before and after, also when that stream does not + // synchronise with the default stream (a non-blocking stream, e.g. a torch side stream). + // Otherwise the default stream, which the helicity streams synchronise with. + gpuStream_t stream = async ? ghelStreams[0] : 0; // Reset the "matrix elements" - running sums of |M|^2 over helicities for the given event + // (the sizes are computed in size_t, from sizeof: the int product could overflow) const int nevt = gpublocks * gputhreads; - gpuMemset( allMEs, 0, nevt * sizeof( fptype ) ); - gpuMemset( ghelAllJamps, 0, cNGoodHel * njampso * mgOnGpu::nx2 * nevt * sizeof( fptype_amp ) ); - gpuMemset( colAllJamp2s, 0, ncolor_flow * nevt * sizeof( fptype_amp ) ); + gpuMemsetAsync( allMEs, 0, sizeof( fptype ) * nevt, stream ); + gpuMemsetAsync( ghelAllJamps, 0, sizeof( fptype_amp ) * cNGoodHel * njampso * mgOnGpu::nx2 * nevt, stream ); + gpuMemsetAsync( colAllJamp2s, 0, sizeof( fptype_amp ) * ncolor_flow * nevt, stream ); // The numerators buffer has NO helicity dimension: all good helicities accumulate in place via // atomicAdd, so it is zeroed once as [nevt][ndiagrams]. The denominators are derived from the // numerators in normalise_output, so the buffer is just [nevt]. - gpuMemset( ghelAllNumerators, 0, ndiagrams * nevt * sizeof( fptype_amp ) ); - gpuMemset( ghelAllDenominators, 0, nevt * sizeof( fptype_amp ) ); - gpuMemset( ghelAllMEs, 0, cNGoodHel * nevt * sizeof( fptype ) ); + gpuMemsetAsync( ghelAllNumerators, 0, sizeof( fptype_amp ) * ndiagrams * nevt, stream ); + gpuMemsetAsync( ghelAllDenominators, 0, sizeof( fptype_amp ) * nevt, stream ); + gpuMemsetAsync( ghelAllMEs, 0, sizeof( fptype ) * cNGoodHel * nevt, stream ); // === PART 1 - HELICITY LOOP: CALCULATE WAVEFUNCTIONS (one event per GPU thread) === // Use CUDA/HIP streams to process different helicities in parallel (one good helicity per stream) // (1) First, within each helicity stream, compute the QCD partial amplitudes jamp's for each helicity // In multichannel mode, also compute the running sums over helicities of numerators, denominators and squared jamp2s + // No good helicity (a subprocess whose |M|^2 vanishes identically, e.g. an FCNC channel with + // every Wilson coefficient at zero): nothing to compute, the |M|^2 stays 0 as on the CPU (a + // launch with gridDim.y = 0 is an invalid configuration) bool storeChannelWeights = allChannelIds != nullptr || allrnddiagram != nullptr; - if( async ) + if( cNGoodHel > 0 && async ) { gpuLaunchKernel2D( calculate_jamps, gpublocks, cNGoodHel, gputhreads, ghelStreams[0], 0, allmomenta, allcouplings, iflavorVec, ghelAllJamps, storeChannelWeights, ghelAllNumerators, ghelAllDenominators, colAllJamp2s, nevt, true ); color_sum_gpu( ghelAllMEs, ghelAllJamps, ghelAllBlasTmp, pBlasHandle, ghelStreams, cNGoodHel, gpublocks, gputhreads, true ); } - else + else if( cNGoodHel > 0 ) { for( int ighel = 0; ighel < cNGoodHel; ighel++ ) { @@ -706,12 +723,12 @@ namespace madmatrix // (3) Wait for all helicity streams to complete, then finally compute the ME sum over all helicities and choose one helicity and one color } // Event-by-event random choice of helicity #403 and ME sum over helicities (defer this after the helicity loop to avoid breaking streams parallelism) - gpuLaunchKernel( add_and_select_hel, gpublocks, gputhreads, allselhel, allrndhel, ghelAllMEs, allMEs, gpublocks * gputhreads, sampleSignedHelicity ); + gpuLaunchKernelStream( add_and_select_hel, gpublocks, gputhreads, stream, allselhel, allrndhel, ghelAllMEs, allMEs, gpublocks * gputhreads, sampleSignedHelicity ); - gpuLaunchKernel( normalise_output, gpublocks, gputhreads, allMEs, iflavorVec, ghelAllNumerators, ghelAllDenominators, allChannelIds, storeChannelWeights, mulChannelWeight, helcolDenominators[0] ); + gpuLaunchKernelStream( normalise_output, gpublocks, gputhreads, stream, allMEs, iflavorVec, ghelAllNumerators, ghelAllDenominators, allChannelIds, storeChannelWeights, mulChannelWeight, helcolDenominators[0] ); // Event-by-event random choice of color and diagram #402 - gpuLaunchKernel( select_col_and_diag, gpublocks, gputhreads, allselcol, allDiagramIdsOut, allrndcol, allrnddiagram, allChannelIds, colAllJamp2s, ghelAllNumerators, ghelAllDenominators, gpublocks * gputhreads ); + gpuLaunchKernelStream( select_col_and_diag, gpublocks, gputhreads, stream, allselcol, allDiagramIdsOut, allrndcol, allrnddiagram, allChannelIds, colAllJamp2s, ghelAllNumerators, ghelAllDenominators, gpublocks * gputhreads ); mgDebugFinalise(); } diff --git a/madgraph/iolibs/template_files/madmatrix/backend/gpu/umami.cc b/madgraph/iolibs/template_files/madmatrix/backend/gpu/umami.cc index d85651573..5c0402444 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/gpu/umami.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/gpu/umami.cc @@ -470,6 +470,10 @@ extern "C" InterfaceInstance* instance = static_cast( handle ); if( !instance->initialized ) { + // the helicity filtering runs on the default stream: wait for the inputs copied on + // gpu_stream, which need not synchronise with it (a non-blocking stream); its + // synchronous copies back to the host finish it before the work below + gpuStreamSynchronize( gpu_stream ); // the jamp scratch of the helicity filtering (one helicity, njampso jamps, a few // events): ghel_jamps, whose size does not depend on ncolor_flow and nampso initialize( diff --git a/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc index 2d78abac8..9f80409c8 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc @@ -198,6 +198,7 @@ namespace madmatrix using A_ACCESS = HostAccessAmplitudes; // TRIVIAL ACCESS (no kernel splitting yet): buffer for one event using CD_ACCESS = HostAccessCouplings; // non-trivial access (dependent couplings): buffer includes all events using CI_ACCESS = HostAccessCouplingsFixed; // TRIVIAL access (independent couplings): buffer for one event + using CDPF_ACCESS = CD_ACCESS; // the running flavor couplings gathered into dpf_value (one event page, as CD_ACCESS) using F_ACCESS = HostAccessIflavorVec; // non-trivial access: buffer includes all events using NUM_ACCESS = HostAccessNumerators; // non-trivial access: buffer includes all events mgDebug( 0, __FUNCTION__ ); @@ -267,8 +268,8 @@ namespace madmatrix // an ordinary value-based view over it. The flavor index is constant across a SIMD lane // (guaranteed by the phase-space integrator), so each lane gets its own running value // while sharing the same flavor selection. This is the direct analogue of Fortran's - // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CD_ACCESS so - // get_coupling_def reads dpf_value with the right per-flavor stride (CD_ACCESS::flv_stride). + // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CDPF_ACCESS + // (= CD_ACCESS here) so get_coupling_def reads dpf_value with the right per-flavor stride. constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); // cppAlign is only defined for SIMD alignas( mgOnGpu::cppAlign ) fptype dpf_value[ndpfbuf]{}; diff --git a/madmatrix/model_handling.py b/madmatrix/model_handling.py index dfaf00419..7282547d5 100644 --- a/madmatrix/model_handling.py +++ b/madmatrix/model_handling.py @@ -2854,10 +2854,11 @@ def format_coupling(self, call): 'COUPs[%s], %s' % (alias[coup], '1.0' if not sign else '-1.0')) elif name == 'flvCOUPs_dep': # dependent (running-alphas) flavored coupling: instantiate the vertex - # routine with CD_ACCESS (the default in the call is CI_ACCESS, as for the - # ordinary running couplings) so get_coupling_def reads the per-event AOSOA - # values gathered into flvCOUPs_dep with the right per-flavor stride. - call = call.replace('CI_ACCESS', 'CD_ACCESS') + # routine with CDPF_ACCESS, the access to the values gathered into + # flvCOUPs_dep (each backend's SigmaKin.cc defines it: CD_ACCESS on the + # CPU, a one-event buffer on the GPU, where a thread is one event), so + # get_coupling_def reads them with the right per-flavor stride. + call = call.replace('CI_ACCESS', 'CDPF_ACCESS') call = call.replace('m_pars->%s%s' % (sign, coup), '%s[%s], %s' % (name, alias[coup], '1.0' if not sign else '-1.0')) elif name == 'flvCOUPs': @@ -2871,6 +2872,10 @@ def format_coupling(self, call): if newcoup: self.couplings2order = self.couporderdep | self.couporderindep + if 'flvCOUPs_dep[' in call and ('CDPF_ACCESS' not in call or re.search(r'\bCOUPs\[|flvCOUPs\[', call)): + # one access type per vertex routine: CDPF_ACCESS cannot also read the + # all-events buffer of an ordinary running coupling (on the GPU) + raise Exception('a vertex mixing a running flavored coupling with another running coupling is not supported: %s' % call) model.cudacpp_wanted_ordered_couplings = self.wanted_ordered_dep_couplings + self.wanted_ordered_indep_couplings + self.wanted_ordered_flv_couplings return call From 41090ea3d5abbe900a322bf1f040019f8ad27fe0 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Tue, 13 Oct 2026 22:34:02 +0200 Subject: [PATCH 2/3] Review fixes: NaN-safe MSSM check, no good helicity in the device kernels, CDPF_ACCESS gather - madmatrix_checks.py mssm_gpu_vs_cpu: the running maximum of the relative differences dropped a nan (rel <= nan is False, so the next finite value replaced it): a nan GPU |M|^2 anywhere but in the last event passed. Every event outside the tolerance is now collected, a nan included, and the count and the first one are reported. - GPU MatrixElementKernelDevice (used by runTest.cc; check_sa.exe and mg7 go through umami): no assert on nGoodHel > 0 any more, sigmaKin handles a subprocess with no good helicity since the previous commit; the per-helicity buffers are sized for one helicity at least, never empty. - CPU/SIMD SigmaKin: the gather into dpf_value, its size and the FLV_COUPLING_ARRAY stride use CDPF_ACCESS, the access the vertex routines read it with, as on the GPU (= CD_ACCESS there: no change of behaviour). --- .github/gpu_ci/madmatrix_checks.py | 21 +++++++++++-------- .../madmatrix/backend/cpu/SigmaKin.cc | 6 +++--- .../backend/gpu/MatrixElementKernels.cc | 12 +++++++---- .../madmatrix/backend/simd/SigmaKin.cc | 6 +++--- 4 files changed, 26 insertions(+), 19 deletions(-) diff --git a/.github/gpu_ci/madmatrix_checks.py b/.github/gpu_ci/madmatrix_checks.py index 8eec47646..215611303 100644 --- a/.github/gpu_ci/madmatrix_checks.py +++ b/.github/gpu_ci/madmatrix_checks.py @@ -28,6 +28,7 @@ import argparse import glob import json +import math import os import random import re @@ -113,15 +114,17 @@ def check_mssm_gpu_vs_cpu(repo, backend, npoints=1024): flavor = [rng.randrange(nflavor) for _ in range(npoints)] (me_cpu,) = cpu(points, alpha_s=alpha_s, flavor=flavor) (me_gpu,) = gpu(points, alpha_s=alpha_s, flavor=flavor) - worst, iworst = 0., -1 - for ievt, (a, b) in enumerate(zip(me_cpu, me_gpu)): - rel = abs(a - b) / abs(a) if a != 0 else (0. if b == 0 else float('inf')) - if not rel <= worst: # also catches nan - worst, iworst = rel, ievt - if not worst < 1e-9: - problems.append('%s: event %d %s %.10e vs %s %.10e (rel %.1e)' - % (sub, iworst, backend, me_gpu[iworst], CPU_BACKEND, - me_cpu[iworst], worst)) + rels = [abs(a - b) / abs(a) if a != 0 else (0. if b == 0 else float('inf')) + for a, b in zip(me_cpu, me_gpu)] + # every event that is not within the tolerance, a nan included (it fails any + # comparison, so it must not go through a running maximum, which would drop it) + bad = [ievt for ievt, rel in enumerate(rels) if not rel < 1e-9] + worst = max((rel for rel in rels if not math.isnan(rel)), default=0.) + if bad: + first = bad[0] + problems.append('%s: %d of %d events differ, first event %d %s %.10e vs %s %.10e (rel %.1e)' + % (sub, len(bad), npoints, first, backend, me_gpu[first], + CPU_BACKEND, me_cpu[first], rels[first])) if not any(v != 0 for v in me_cpu): problems.append('%s: all |M|^2 vanish' % sub) scaling = me_gpu[1] / me_gpu[0] / (0.12 / 0.10) ** 2 - 1 if me_gpu[0] else float('inf') diff --git a/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc index 3e1ef631f..214075cf2 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc @@ -269,7 +269,7 @@ namespace madmatrix // while sharing the same flavor selection. This is the direct analogue of Fortran's // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CDPF_ACCESS // (= CD_ACCESS here) so get_coupling_def reads dpf_value with the right per-flavor stride. - constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); + constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CDPF_ACCESS::flv_stride : 1 ); // cppAlign is only defined for SIMD alignas( mgOnGpu::cppAlign ) fptype dpf_value[ndpfbuf]{}; for( int idpf = 0; idpf < nDPF; idpf++ ) @@ -277,10 +277,10 @@ namespace madmatrix { const int idc = cDPF_idcoup[idpf * nMF + imf]; if( idc >= 0 ) - CD_ACCESS::kernelAccess( dpf_value + ( idpf * nMF + imf ) * CD_ACCESS::flv_stride ) = + CDPF_ACCESS::kernelAccess( dpf_value + ( idpf * nMF + imf ) * CDPF_ACCESS::flv_stride ) = CD_ACCESS::kernelAccessConst( COUPs[idc] ); } - FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; + FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; // Reset color flows (reset jamp_sv) at the beginning of a new event or event page for( int i = 0; i < njampso; i++ ) { jamp_sv[i] = cxzero_sv(); } diff --git a/madgraph/iolibs/template_files/madmatrix/backend/gpu/MatrixElementKernels.cc b/madgraph/iolibs/template_files/madmatrix/backend/gpu/MatrixElementKernels.cc index 419241e5c..9600010c8 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/gpu/MatrixElementKernels.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/gpu/MatrixElementKernels.cc @@ -14,6 +14,7 @@ #include "MemoryBuffers.h" #include "color_sum.h" // for blasColorSumTmpSize +#include // for std::max #include // for fetestexcept #include #include @@ -309,7 +310,10 @@ namespace madmatrix sigmaKin_getGoodHel( m_momenta.data(), m_couplings.data(), m_iflavorVec.data(), m_matrixElements.data(), m_pHelNumerators->data(), m_pHelDenominators->data(), m_pHelJamps->data(), hstIsGoodHel.data(), nevt ); // ... 0d3. Set good helicity list in host static memory int nGoodHel = sigmaKin_setGoodHel( hstIsGoodHel.data() ); - assert( nGoodHel > 0 ); // SANITY CHECK: there should be at least one good helicity + // nGoodHel can be 0 (a subprocess whose |M|^2 vanishes identically): sigmaKin then skips the + // helicity loop and returns |M|^2 = 0, as on the CPU. The per-helicity buffers below are + // sized for one helicity at least, so that they are never empty (null) allocations. + const int nGoodHelAlloc = std::max( nGoodHel, 1 ); // Create one GPU stream for each good helicity for( int ighel = 0; ighel < nGoodHel; ighel++ ) gpuStreamCreate( &m_helStreams[ighel] ); @@ -325,10 +329,10 @@ namespace madmatrix } #endif // ... Create the "many-helicity" super-buffer of nGoodHel ME buffers (dynamically allocated because nGoodHel is determined at runtime) - m_pHelMEs.reset( new DeviceBufferSimple( nGoodHel * nevt ) ); + m_pHelMEs.reset( new DeviceBufferSimple( nGoodHelAlloc * nevt ) ); // ... Create the "many-helicity" super-buffer of nGoodHel ME buffers (dynamically allocated because nGoodHel is determined at runtime) // ... (calling reset here deletes the previously created "one-helicity" buffers used for helicity filtering) - m_pHelJamps.reset( new DeviceBufferAmp( nGoodHel * ProcessData::njampso * mgOnGpu::nx2 * nevt ) ); + m_pHelJamps.reset( new DeviceBufferAmp( nGoodHelAlloc * ProcessData::njampso * mgOnGpu::nx2 * nevt ) ); // ... Create the numerator and denominator buffers. These no longer carry a helicity dimension: // ... the numerators are accumulated in place over all good helicities via atomicAdd in calculate_jamps // ... ([nevt][ndiagrams]) and the denominators are derived from them ([nevt]). @@ -338,7 +342,7 @@ namespace madmatrix // Create the "many-helicity" super-buffer of temporary buffers for the cuBLAS/hipBLAS intermediate // results in color_sum_blas, and in mixed precision mode for the converted jamps too // (see blasColorSumTmpSize in color_sum.h, which is where the size is defined) - if( m_blasColorSum ) m_pHelBlasTmp.reset( new DeviceBufferSimple2( blasColorSumTmpSize( nGoodHel, nevt ) ) ); + if( m_blasColorSum ) m_pHelBlasTmp.reset( new DeviceBufferSimple2( blasColorSumTmpSize( nGoodHelAlloc, nevt ) ) ); #endif // Return the number of good helicities return nGoodHel; diff --git a/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc index 9f80409c8..969e0a8b2 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc +++ b/madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc @@ -270,7 +270,7 @@ namespace madmatrix // while sharing the same flavor selection. This is the direct analogue of Fortran's // FLV_xx%VAL(k)%P => GC_yyy(J). The vertex routines are instantiated with CDPF_ACCESS // (= CD_ACCESS here) so get_coupling_def reads dpf_value with the right per-flavor stride. - constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); + constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CDPF_ACCESS::flv_stride : 1 ); // cppAlign is only defined for SIMD alignas( mgOnGpu::cppAlign ) fptype dpf_value[ndpfbuf]{}; for( int idpf = 0; idpf < nDPF; idpf++ ) @@ -278,10 +278,10 @@ namespace madmatrix { const int idc = cDPF_idcoup[idpf * nMF + imf]; if( idc >= 0 ) - CD_ACCESS::kernelAccess( dpf_value + ( idpf * nMF + imf ) * CD_ACCESS::flv_stride ) = + CDPF_ACCESS::kernelAccess( dpf_value + ( idpf * nMF + imf ) * CDPF_ACCESS::flv_stride ) = CD_ACCESS::kernelAccessConst( COUPs[idc] ); } - FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; + FLV_COUPLING_ARRAY flvCOUPs_dep{ cDPF_partner1, cDPF_partner2, dpf_value }; // Reset color flows (reset jamp_sv) at the beginning of a new event or event page for( int i = 0; i < njampso; i++ ) { jamp_sv[i] = cxzero_sv(); } From 8974fe3beb23790515aae59ae8ae16511203c133 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Wed, 14 Oct 2026 14:11:16 +0200 Subject: [PATCH 3/3] GPU CI: run the madmatrix and interference checks with check_sa.exe, not ctypes Review of #255 (Daniele): check_sa.exe already evaluates madmatrix through umami on every backend, so use it instead of a second umami driver written with ctypes, which copied the key values of umami.h by hand. check_sa.exe perf gets per-event inputs and outputs as raw binary files, like --momenta and --dump-me: - --alphas , --flavors : alpha_s and flavour of each event; - --rnd-hel , --rnd-col : the random numbers of the helicity and colour choice; --dump-hel , --dump-col : the chosen indices; - --param =: umami_set_parameter (repeatable); - --nonblocking-stream (CUDA/HIP): umami runs on a non-blocking stream; the driver synchronises the device before the call and the stream after it. run_umami takes the key and buffer lists; --rambo-massless no longer eats the next argument. New gpuStreamCreateNonBlocking in GpuAbstraction.h. madmatrix_checks.py builds check_sa_.exe for the scalar and the GPU backend and compares their dumps on the same events: MSSM p p > go go with a random alpha_s and flavour per event, plus the alpha_s^2 scaling of every event (--alphas 0.10 / 0.12); b b~ > ta+ ta- / z a with ymtau = 0 on both backends; g g > t t~ g, 4 calls of 16384 events, non-blocking against default stream. A GPU abort now fails one check_sa.exe run, so the per-check subprocess re-exec is gone. interference_checks.py draws the helicity of u u~ > t t~ g QED^2==2 the same way (16 points x 4000 grid random numbers, interference_helicity_summed via --param). umami_harness.py is removed. Checked on CPU (simd_128 against scalar): mssm max rel diff 2.5e-15 with nDPF = 16, zero_me |M|^2 = 0 / helicity -1 on both, helicity_choice 15 of 16 points with mixed signs, helicity_summed exact. --- .github/gpu_ci/README.md | 4 +- .github/gpu_ci/interference_checks.py | 129 ++++---- .github/gpu_ci/madmatrix_checks.py | 301 +++++++++-------- .github/gpu_ci/madmatrix_gpu.sh | 2 +- .github/gpu_ci/umami_harness.py | 183 ---------- .github/workflows/gpu_runner_ci.yml | 4 +- .../madmatrix/backend/gpu/GpuAbstraction.h | 2 + .../template_files/madmatrix/check_sa.cc | 313 ++++++++++++++---- 8 files changed, 483 insertions(+), 455 deletions(-) delete mode 100644 .github/gpu_ci/umami_harness.py diff --git a/.github/gpu_ci/README.md b/.github/gpu_ci/README.md index fb0ab447d..d56dde6b9 100644 --- a/.github/gpu_ci/README.md +++ b/.github/gpu_ci/README.md @@ -23,8 +23,8 @@ allocation per workflow run, and every CI job of that run executes in it: | `pp_ttx` | cluster | `generate p p > t t~`, `output mg7`, `device = ["cuda"]` or `["hip"]` in `run_card.toml`, `bin/generate_events -f`, check the cross section in `info.json` (`pp_ttx_mg7.sh`). The logs and cards are uploaded as an artifact, and the cross section is shown in the run summary | | `pp_jj` | cluster | the same with `generate p p > j j` (`PROCESS`/`TAG` of `pp_ttx_mg7.sh`): its jet cuts act on several objects, which `p p > t t~` has none of | | `madspace_tests` | cluster | the madspace unit tests, `madspace/tests`, with `MADSPACE_TEST_DEVICE=cuda` or `hip` (`madspace_tests_gpu.sh`): the tests put their torch tensors on the GPU, so madspace runs them on its GPU runtime. Host-side tests (marker `cpu_only`) are skipped. pytest and scipy are installed in the cache when the modules lack them, and so is a CUDA torch wheel when the PyTorch module is a CPU build (manneback). The log and `junit.xml` are uploaded as an artifact | -| `interference` | cluster | squared split orders on the GPU backend of madmatrix (`interference_gpu.sh`, `interference_checks.py`): standalone `u u~ > u u~` constraints, the helicity choice of an interference matrix element through umami (exact and summed), and `p p > u u~ QCD^2==2` with mg7 against its CPU reference. Input `interference`, on by default and on push | -| `madmatrix` | cluster | the GPU backend of madmatrix against its scalar CPU backend through umami (`madmatrix_gpu.sh`, `madmatrix_checks.py`): MSSM `p p > go go` event by event (running flavour couplings), a subprocess with no good helicity (`b b~ > ta+ ta- / z a`, ymtau = 0), and a non-blocking stream handed to umami (`g g > t t~ g`). Each check runs in its own process | +| `interference` | cluster | squared split orders on the GPU backend of madmatrix (`interference_gpu.sh`, `interference_checks.py`): standalone `u u~ > u u~` constraints, the helicity choice of an interference matrix element through umami, with `check_sa.exe` (exact and summed), and `p p > u u~ QCD^2==2` with mg7 against its CPU reference. Input `interference`, on by default and on push | +| `madmatrix` | cluster | the GPU backend of madmatrix against its scalar CPU backend (`madmatrix_gpu.sh`, `madmatrix_checks.py`), with the `check_sa.exe` of each backend fed the same events from files: MSSM `p p > go go` event by event (running flavour couplings), a subprocess with no good helicity (`b b~ > ta+ ta- / z a`, ymtau = 0), and a non-blocking stream handed to umami (`g g > t t~ g`, `--nonblocking-stream`) | | `stop_runner` | cluster | clean up, then create the stop file: `runner_batch.sh` stops the runner and the allocation ends | | `watchdog` | vmadgraph | from the start of the runner to the end of the allocation: if the allocation ends before `stop_runner` completed, cancel the run (`ceci-gate watch-runner`, `remote_watch_runner.sh`) | diff --git a/.github/gpu_ci/interference_checks.py b/.github/gpu_ci/interference_checks.py index 28f4ce338..3f8d49173 100755 --- a/.github/gpu_ci/interference_checks.py +++ b/.github/gpu_ci/interference_checks.py @@ -12,12 +12,13 @@ pair loop of the colour sum (all ORDERED pairs of amplitude orders, not a triangle); and QED^2==2 again with FPTYPE=m (colour algebra in single precision). 2. the helicity drawn for an interference |M|^2 through umami (event generation), on - ``u u~ > t t~ g QED^2==2``, whose helicities have both signs at most points: every - draw has |M|^2 = sum|T|, a helicity always comes with the same sign, the average over - a uniform grid of random numbers is the helicity sum (what umami returns without the - random number) and the points with mixed signs give both signs; with the umami - parameter interference_helicity_summed (run card interference_helicity = "summed") - every draw returns the helicity sum. + ``u u~ > t t~ g QED^2==2``, whose helicities have both signs at most points, with + check_sa.exe perf on a grid of random numbers for each RAMBO point (--momenta, + --rnd-hel, --dump-me, --dump-hel): every draw has |M|^2 = sum|T|, a helicity always + comes with the same sign, the average over the uniform grid is the helicity sum (what + umami returns without the random number) and the points with mixed signs give both + signs; with the umami parameter interference_helicity_summed (run card + interference_helicity = "summed", --param) every draw returns the helicity sum. 3. mg7 ``p p > u u~ QCD^2==2`` on the backend: the negative cross section of the CPU reference, the LHE declaring the negative weights (IDWTUP = -4, XMAXUP the unit weight sigma_abs) and events of both signs. It runs with the PDF set of the @@ -29,19 +30,18 @@ """ import argparse -import ctypes +import array import glob import gzip import json import math import os -import random import re import subprocess import sys import time -from umami_harness import GPU_BACKENDS, Memory, rambo +GPU_BACKENDS = ('cuda', 'hip') # u u~ > u u~ QED^2==2 at the check_sa RAMBO point, from the Fortran split-order # standalone (madmatrix CPU: -5.5828746494657258e-02) @@ -138,60 +138,53 @@ def check_standalone(repo, backend): def check_helicity_choice(repo, backend, npoints=16, ngrid=4000): section('helicity choice of an interference |M|^2 through umami (%s)' % backend) proc_dir = standalone(repo, backend, 'uuxttxg_qed2_d', 'u u~ > t t~ g QED^2==2', 'd') - out_dir = os.path.dirname(os.path.dirname(proc_dir)) - # the process library finds its common library through its own rpath - library = glob.glob(os.path.join(out_dir, 'lib', 'libmadmatrix_P*.so'))[0] - lib = ctypes.CDLL(library) - memory = Memory(backend, library) - handle = ctypes.c_void_p() - card = os.path.join(out_dir, 'Cards', 'param_card.dat').encode() - assert lib.umami_initialize(ctypes.byref(handle), card) == 0, 'umami_initialize' - npar = ctypes.c_int() - lib.umami_get_meta(1, ctypes.byref(npar)) # UMAMI_META_PARTICLE_COUNT - npar = npar.value - masses = (ctypes.c_double * npar)() - assert lib.umami_get_meta(5, masses) == 0, 'UMAMI_META_MASSES' # after initialize - # umami.h: UMAMI_IN_MOMENTA = 0, UMAMI_IN_RANDOM_HELICITY = 4, - # UMAMI_OUT_MATRIX_ELEMENT = 0, UMAMI_OUT_HELICITY_INDEX = 3 - lib.umami_matrix_element.argtypes = [ - ctypes.c_void_p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_size_t, - ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p), - ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p)] - - def evaluate(point, draw): - count = ngrid if draw else 1 - mom = (ctypes.c_double * (4 * npar * count))() - for ipart in range(npar): # momenta[count * (npar * mu + ipart) + ievt] - for mu in range(4): - start = count * (npar * mu + ipart) - mom[start:start + count] = [point[ipart][mu]] * count - rnd = (ctypes.c_double * count)(*[(k + 0.5) / count for k in range(count)]) - me = (ctypes.c_double * count)() - hel = (ctypes.c_int * count)() - inputs = [memory.buffer(mom)] + ([memory.buffer(rnd)] if draw else []) - in_keys = [0] + ([4] if draw else []) - outputs = [memory.buffer(me)] + ([memory.buffer(hel)] if draw else []) - out_keys = [0] + ([3] if draw else []) - status = lib.umami_matrix_element( - handle, count, count, 0, - len(in_keys), (ctypes.c_int * len(in_keys))(*in_keys), - (ctypes.c_void_p * len(inputs))(*inputs), - len(out_keys), (ctypes.c_int * len(out_keys))(*out_keys), - (ctypes.c_void_p * len(outputs))(*outputs)) - assert status == 0, 'umami_matrix_element status %d' % status - memory.fetch(outputs[0], me) - if draw: - memory.fetch(outputs[1], hel) - memory.release() - return list(me), list(hel) - - rng = random.Random(7) + work = os.path.abspath('uuxttxg_qed2_d') + + def check_sa(args, log): + status = run(['./check_sa.exe', 'perf'] + [str(a) for a in args], cwd=proc_dir, + log=log) + if status != 0: + raise RuntimeError('check_sa.exe perf %s failed, see %s' + % (' '.join(str(a) for a in args), log)) + + def doubles(path): + values = array.array('d') + with open(path, 'rb') as f: + values.frombytes(f.read()) + return values + + # npoints RAMBO points at 1000 GeV, each repeated ngrid times, with a uniform grid of + # helicity random numbers (k + 0.5) / ngrid for each point: one batch of 64000 events + points = os.path.join(work, 'points.bin') + check_sa(['--energy', 1000, '--dump-momenta', points, 1, npoints, 1], 'uuxttxg_points.log') + with open(points, 'rb') as f: + raw = f.read() + size = len(raw) // npoints + grid = os.path.join(work, 'grid.bin') + with open(grid, 'wb') as f: + for ipoint in range(npoints): + f.write(raw[ipoint * size:(ipoint + 1) * size] * ngrid) + rnd = os.path.join(work, 'rnd.bin') + with open(rnd, 'wb') as f: + array.array('d', [(k + 0.5) / ngrid for k in range(ngrid)] * npoints).tofile(f) + batch = [npoints * ngrid // 256, 256, 1] + outputs = {} + for tag, extra in (('sum', []), ('draw', ['--rnd-hel', rnd]), + ('summed', ['--rnd-hel', rnd, '--param', 'interference_helicity_summed=1'])): + me, hel = os.path.join(work, 'me_%s.bin' % tag), os.path.join(work, 'hel_%s.bin' % tag) + check_sa(['--momenta', grid, '--dump-me', me, '--dump-hel', hel] + extra + batch, + 'uuxttxg_%s.log' % tag) + helicities = array.array('i') + with open(hel, 'rb') as f: + helicities.frombytes(f.read()) + outputs[tag] = (doubles(me), helicities) + mixed = 0 problems = [] for ipoint in range(npoints): - point = rambo(list(masses)[2:], 1000., rng) - (helicity_sum,), _ = evaluate(point, False) - me, hel = evaluate(point, True) + events = slice(ipoint * ngrid, (ipoint + 1) * ngrid) + helicity_sum = outputs['sum'][0][ipoint * ngrid] + me, hel = outputs['draw'][0][events], outputs['draw'][1][events] abs_sum = abs(me[0]) signs = {} for value, ihel in zip(me, hel): @@ -219,19 +212,9 @@ def evaluate(point, draw): # run card [generation] interference_helicity = "summed": the instance keeps the # helicity sum as the |M|^2 of every draw (mg7 then writes the helicities as 9) - lib.umami_set_parameter.argtypes = [ctypes.c_void_p, ctypes.c_char_p, - ctypes.c_double, ctypes.c_double] - status = lib.umami_set_parameter(handle, b'interference_helicity_summed', 1., 0.) - worst = 0. - for ipoint in range(4): - point = rambo(list(masses)[2:], 1000., rng) - (helicity_sum,), _ = evaluate(point, False) - me, _hel = evaluate(point, True) - worst = max([worst] + [abs(v - helicity_sum) / abs(helicity_sum) for v in me]) - check('helicity_summed', status == 0 and worst < 1e-12, - 'umami_set_parameter status %d, max |M|^2 / helicity sum - 1 = %.1e' - % (status, worst)) - lib.umami_free(handle) + worst = max(abs(v - s) / abs(s) for v, s in zip(outputs['summed'][0], outputs['sum'][0])) + check('helicity_summed', worst < 1e-12, + 'max |M|^2 / helicity sum - 1 = %.1e over %d draws' % (worst, npoints * ngrid)) # --------------------------------------------------------------------------- diff --git a/.github/gpu_ci/madmatrix_checks.py b/.github/gpu_ci/madmatrix_checks.py index 215611303..797d9f78f 100644 --- a/.github/gpu_ci/madmatrix_checks.py +++ b/.github/gpu_ci/madmatrix_checks.py @@ -1,44 +1,50 @@ #!/usr/bin/env python3 -"""The GPU backend of madmatrix against its scalar CPU backend, through umami. +"""The GPU backend of madmatrix against its scalar CPU backend, with check_sa.exe. Run by madmatrix_gpu.sh (job madmatrix of gpu_runner_ci.yml) in its run directory. +Every evaluation is a ``check_sa.exe perf`` run of the standalone output, one executable +per backend (check_sa_.exe), on the same events: the momenta, alpha_s, flavour +and random numbers of every event go in as raw binary files (--momenta, --alphas, +--flavors, --rnd-hel, --rnd-col) and the |M|^2, helicities and colours come out the +same way (--dump-me, --dump-hel, --dump-col). An abort of the GPU runtime fails the run +of check_sa.exe, so only the check it belongs to. 1. mssm_gpu_vs_cpu: MSSM ``p p > go go`` (both subprocesses), whose squark-quark-gluino couplings are running flavour couplings (nDPF > 0, gathered per event into - dpf_value), on random points with a random alpha_s and flavour per event: the GPU + dpf_value), on RAMBO points with a random alpha_s and flavour per event: the GPU |M|^2 equal the CPU ones event by event (the GPU gather used to address its one-event dpf_value with the offset of the event in the all-events coupling buffer, out of - bounds past the first event page), and scale as alpha_s^2. -2. zero_me: ``b b~ > ta+ ta- / z a`` (s-channel Higgs only) with the tau Yukawa coupling - set to 0 in the param card, so that no helicity survives the helicity filtering: the - GPU library returns |M|^2 = 0 and no helicity (-1), and for the colour either none - (-1) or the only colour flow (0, what a -ffast-math build makes of 0/0, as the CPU - one), never stale memory; it used to abort on the invalid launch configuration + bounds past the first event page), and scale as alpha_s^2 (--alphas 0.10 and 0.12). +2. zero_me: ``b b~ > ta+ ta- / z a`` (s-channel Higgs only): |M|^2 > 0 with the default + param card, then the tau Yukawa coupling set to 0, so that no helicity survives the + helicity filtering: both backends return |M|^2 = 0 and no helicity (-1), and for the + colour either none (-1) or the only colour flow (0, what a -ffast-math build makes of + 0/0), never stale memory; the GPU used to abort on the invalid launch configuration (gridDim.y = 0) of the helicity loop. -3. nonblocking_stream: ``g g > t t~ g`` with umami handed a non-blocking stream (as a - torch side stream, which does not synchronise with the default stream) gives the same - |M|^2 and helicities as on the default stream, call after call: every memset, kernel - and copy of the matrix element must go to the stream umami is given. - -Each check runs in a process of its own (--only NAME): an abort of the GPU runtime -fails that check, not the others. Writes summary.txt (key=value lines) in the current -directory; exit status 1 if a check fails. +3. nonblocking_stream: ``g g > t t~ g``, 4 calls of 16384 events, with umami handed a + non-blocking stream (--nonblocking-stream, as a torch side stream, which does not + synchronise with the default stream) gives the same |M|^2 and helicities as on the + default stream: every memset, kernel and copy of the matrix element must go to the + stream umami is given, the helicity filtering of the first call included. + +Any CPU backend but scalar (simd_128, ...) works too, to try the script without a GPU; +check 3 is then skipped. Writes summary.txt (key=value lines) in the current directory; +exit status 1 if a check fails. """ import argparse +import array import glob -import json import math import os import random import re +import shutil import subprocess import sys import time -from umami_harness import (GPU_BACKENDS, OUT_COLOR_INDEX, OUT_HELICITY_INDEX, - OUT_MATRIX_ELEMENT, Umami, rambo) - +GPU_BACKENDS = ('cuda', 'hip') CPU_BACKEND = 'scalar' CHECKS = ('mssm_gpu_vs_cpu', 'zero_me', 'nonblocking_stream') @@ -59,63 +65,96 @@ def run(cmd, cwd=None, log=None): return subprocess.call(cmd, cwd=cwd, stdout=out, stderr=subprocess.STDOUT) +def write(path, typecode, values): + """A raw binary file of check_sa.exe: 'd' double, 'I' 32-bit unsigned int.""" + with open(path, 'wb') as f: + array.array(typecode, values).tofile(f) + return os.path.abspath(path) + + +def read(path, typecode): + """A file written by check_sa.exe: 'd' (--dump-me), 'i' (--dump-hel, --dump-col).""" + values = array.array(typecode) + with open(path, 'rb') as f: + values.frombytes(f.read()) + return values + + def output(repo, name, model, process, backends): - """output standalone + make for these backends in every subprocess: the output - directory and its subprocess names.""" + """output standalone + make for these backends in every subprocess, keeping the + check_sa.exe of each as check_sa_.exe: the subprocess directories.""" with open(name + '.mg5', 'w') as f: f.write('import model %s\ngenerate %s\noutput standalone %s -f\n' % (model, process, name)) if run([sys.executable, os.path.join(repo, 'bin', 'madgraph'), name + '.mg5'], log=name + '.log') != 0: raise RuntimeError('madgraph failed, see %s.log' % name) - subprocesses = sorted(os.path.basename(d) - for d in glob.glob(os.path.join(name, 'SubProcesses', 'P*'))) - if not subprocesses: + proc_dirs = sorted(glob.glob(os.path.join(name, 'SubProcesses', 'P*'))) + if not proc_dirs: raise RuntimeError('no subprocess directory for %s' % process) jobs = os.environ.get('SLURM_CPUS_PER_TASK', '4') - for sub in subprocesses: + for proc_dir in proc_dirs: for backend in backends: - log = '%s_%s_%s_make.log' % (name, sub, backend) + log = '%s_%s_%s_make.log' % (name, os.path.basename(proc_dir), backend) + # check_sa.exe is relinked for the backend of each make (its libraries carry + # the backend in their name) if run(['make', '-j', jobs, 'BACKEND=' + backend, 'FPTYPE=d'], - cwd=os.path.join(name, 'SubProcesses', sub), log=log) != 0: + cwd=proc_dir, log=log) != 0: raise RuntimeError('make failed, see %s' % log) - return subprocesses + shutil.copy(os.path.join(proc_dir, 'check_sa.exe'), + os.path.join(proc_dir, 'check_sa_%s.exe' % backend)) + return proc_dirs -def library(name, sub, backend): - return os.path.join(name, 'lib', 'libmadmatrix_%s_%s.so' % (sub, backend)) +def check_sa(proc_dir, backend, args, log): + """check_sa_.exe perf in proc_dir (it reads ../../Cards/param_card.dat).""" + cmd = ['./check_sa_%s.exe' % backend, 'perf'] + [str(a) for a in args] + status = run(cmd, cwd=proc_dir, log=log) + if status != 0: + raise RuntimeError('%s failed (exit status %d), see %s' % (' '.join(cmd), status, log)) -def card(name): - return os.path.join(name, 'Cards', 'param_card.dat') +def dump(prefix, what): + return os.path.abspath('%s_%s.bin' % (prefix, what)) # --------------------------------------------------------------------------- # 1. running flavour couplings (MSSM p p > go go): GPU against CPU # --------------------------------------------------------------------------- -def check_mssm_gpu_vs_cpu(repo, backend, npoints=1024): +def check_mssm_gpu_vs_cpu(backend, repo, nevt=1024): section('MSSM p p > go go, running flavour couplings: %s against %s' % (backend, CPU_BACKEND)) name = 'mssm_gogo' - subprocesses = output(repo, name, 'MSSM_SLHA2', 'p p > go go', (CPU_BACKEND, backend)) + proc_dirs = output(repo, name, 'MSSM_SLHA2', 'p p > go go', (CPU_BACKEND, backend)) rng = random.Random(11) problems, details = [], [] - for sub in subprocesses: - data = open(os.path.join(name, 'SubProcesses', sub, 'ProcessData.h')).read() + for proc_dir in proc_dirs: + sub = os.path.basename(proc_dir) + data = open(os.path.join(proc_dir, 'ProcessData.h')).read() ndpf = int(re.search(r'constexpr int nDPF = (\d+);', data).group(1)) nflavor = int(re.search(r'constexpr int nmaxflavor = (\d+);', data).group(1)) if ndpf == 0: problems.append('%s: nDPF = 0, no running flavour coupling to check' % sub) - cpu = Umami(CPU_BACKEND, library(name, sub, CPU_BACKEND), card(name)) - gpu = Umami(backend, library(name, sub, backend), card(name)) - points = [rambo(gpu.masses[2:], 3000., rng, gpu.masses[:2]) for _ in range(npoints)] - # the first two events: one point at two values of alpha_s - points[1] = points[0] - alpha_s = [0.10, 0.12] + [rng.uniform(0.09, 0.14) for _ in range(npoints - 2)] - flavor = [rng.randrange(nflavor) for _ in range(npoints)] - (me_cpu,) = cpu(points, alpha_s=alpha_s, flavor=flavor) - (me_gpu,) = gpu(points, alpha_s=alpha_s, flavor=flavor) + prefix = '%s_%s' % (name, sub) + # one batch of nevt events, the same on both backends + batch = [1, nevt, 1] + momenta = dump(prefix, 'momenta') + check_sa(proc_dir, CPU_BACKEND, ['--energy', 3000, '--dump-momenta', momenta] + batch, + prefix + '_rambo.log') + alphas = write(prefix + '_alphas.bin', 'd', [rng.uniform(0.09, 0.14) for _ in range(nevt)]) + flavors = write(prefix + '_flavors.bin', 'I', [rng.randrange(nflavor) for _ in range(nevt)]) + me = {} + for tag, run_backend, alphas_arg in (('cpu', CPU_BACKEND, alphas), ('gpu', backend, alphas), + ('gpu010', backend, '0.10'), ('gpu012', backend, '0.12')): + out = dump(prefix, 'me_' + tag) + check_sa(proc_dir, run_backend, ['--momenta', momenta, '--alphas', alphas_arg, + '--flavors', flavors, '--dump-me', out] + batch, + '%s_%s.log' % (prefix, tag)) + me[tag] = read(out, 'd') + if not len(me['cpu']) == len(me['gpu']) == nevt: + problems.append('%s: %d CPU and %d GPU |M|^2 for %d events' + % (sub, len(me['cpu']), len(me['gpu']), nevt)) rels = [abs(a - b) / abs(a) if a != 0 else (0. if b == 0 else float('inf')) - for a, b in zip(me_cpu, me_gpu)] + for a, b in zip(me['cpu'], me['gpu'])] # every event that is not within the tolerance, a nan included (it fails any # comparison, so it must not go through a running maximum, which would drop it) bad = [ievt for ievt, rel in enumerate(rels) if not rel < 1e-9] @@ -123,17 +162,20 @@ def check_mssm_gpu_vs_cpu(repo, backend, npoints=1024): if bad: first = bad[0] problems.append('%s: %d of %d events differ, first event %d %s %.10e vs %s %.10e (rel %.1e)' - % (sub, len(bad), npoints, first, backend, me_gpu[first], - CPU_BACKEND, me_cpu[first], rels[first])) - if not any(v != 0 for v in me_cpu): + % (sub, len(bad), nevt, first, backend, me['gpu'][first], + CPU_BACKEND, me['cpu'][first], rels[first])) + if not any(v != 0 for v in me['cpu']): problems.append('%s: all |M|^2 vanish' % sub) - scaling = me_gpu[1] / me_gpu[0] / (0.12 / 0.10) ** 2 - 1 if me_gpu[0] else float('inf') - if not abs(scaling) < 1e-9: - problems.append('%s: |M|^2 not ~ alpha_s^2 (ratio - 1 = %.1e)' % (sub, scaling)) + # |M|^2 ~ alpha_s^2, event by event + ratio = (0.12 / 0.10) ** 2 + scaling = [b / a / ratio - 1 if a != 0 else float('inf') + for a, b in zip(me['gpu010'], me['gpu012'])] + bad = [ievt for ievt, s in enumerate(scaling) if not abs(s) < 1e-9] + if bad or len(scaling) != nevt: + problems.append('%s: |M|^2 not ~ alpha_s^2 in %d of %d events (first: ratio - 1 = %.1e)' + % (sub, len(bad), nevt, scaling[bad[0]] if bad else float('nan'))) details.append('%s nDPF %d, %d flavours, max rel diff %.1e' % (sub, ndpf, nflavor, worst)) print(' ' + details[-1]) - cpu.free() - gpu.free() check('mssm_gpu_vs_cpu', not problems, '; '.join(problems[:3]) or '; '.join(details)) @@ -141,109 +183,108 @@ def check_mssm_gpu_vs_cpu(repo, backend, npoints=1024): # 2. a subprocess with no good helicity # --------------------------------------------------------------------------- -def check_zero_me(repo, backend, npoints=64): - section('no good helicity: b b~ > ta+ ta- / z a with ymtau = 0 (%s)' % backend) +def check_zero_me(backend, repo, nevt=64): + section('no good helicity: b b~ > ta+ ta- / z a with ymtau = 0 (%s, %s)' % (backend, CPU_BACKEND)) name = 'zero_me' - (sub,) = output(repo, name, 'sm', 'b b~ > ta+ ta- / z a', (CPU_BACKEND, backend)) - text = open(card(name)).read() + (proc_dir,) = output(repo, name, 'sm', 'b b~ > ta+ ta- / z a', (CPU_BACKEND, backend)) + card = os.path.join(name, 'Cards', 'param_card.dat') + text = open(card).read() zero_card, count = re.subn(r'(?mi)^(\s*15\s+)\S+(\s*#\s*ymtau)', r'\g<1>0.000000e+00\2', text) if count != 1: - raise RuntimeError('no ymtau in %s' % card(name)) - zero_path = os.path.join(name, 'Cards', 'param_card_ymtau0.dat') - open(zero_path, 'w').write(zero_card) + raise RuntimeError('no ymtau in %s' % card) rng = random.Random(5) - # the helicity filtering runs once per library (and process): the GPU library only - # ever sees the card with ymtau = 0, the CPU library the default one (|M|^2 != 0) - gpu = Umami(backend, library(name, sub, backend), zero_path) - cpu = Umami(CPU_BACKEND, library(name, sub, CPU_BACKEND), card(name)) - points = [rambo(gpu.masses[2:], 500., rng, gpu.masses[:2]) for _ in range(npoints)] - rnd = [rng.random() for _ in range(npoints)] - outputs = (OUT_MATRIX_ELEMENT, OUT_HELICITY_INDEX, OUT_COLOR_INDEX) - me, hel, col = gpu(points, random_helicity=rnd, random_color=rnd, outputs=outputs) - (me_default,) = cpu(points) - problems = [] - if any(v != 0 for v in me): - problems.append('|M|^2 not 0: %s' % me[:4]) - if set(hel) != {-1} or not set(col) <= {-1, 0}: - problems.append('helicities %s, colours %s (expected -1, and -1 or 0)' - % (sorted(set(hel))[:4], sorted(set(col))[:4])) - if not all(v > 0 for v in me_default): - problems.append('|M|^2 with the default card not > 0: %s' % me_default[:4]) - gpu.free() - cpu.free() + batch = [1, nevt, 1] + momenta = dump(name, 'momenta') + check_sa(proc_dir, CPU_BACKEND, ['--energy', 500, '--dump-momenta', momenta, + '--dump-me', dump(name, 'me_default')] + batch, + name + '_default.log') + me_default = read(dump(name, 'me_default'), 'd') + rnd = write(name + '_rnd.bin', 'd', [rng.random() for _ in range(nevt)]) + # check_sa.exe reads the param card at every run: from here on ymtau = 0 (the + # helicity filtering runs at the first call of each run) + open(card, 'w').write(zero_card) + problems, seen = [], [] + for run_backend in (backend, CPU_BACKEND): + prefix = '%s_%s' % (name, run_backend) + check_sa(proc_dir, run_backend, ['--momenta', momenta, '--rnd-hel', rnd, '--rnd-col', rnd, + '--dump-me', dump(prefix, 'me'), + '--dump-hel', dump(prefix, 'hel'), + '--dump-col', dump(prefix, 'col')] + batch, + prefix + '.log') + me, hel, col = (read(dump(prefix, 'me'), 'd'), read(dump(prefix, 'hel'), 'i'), + read(dump(prefix, 'col'), 'i')) + if len(me) != nevt or any(v != 0 for v in me): + problems.append('%s: |M|^2 not 0: %s' % (run_backend, list(me[:4]))) + if set(hel) != {-1} or not set(col) <= {-1, 0}: + problems.append('%s: helicities %s, colours %s (expected -1, and -1 or 0)' + % (run_backend, sorted(set(hel))[:4], sorted(set(col))[:4])) + seen.append('%s colour %s' % (run_backend, sorted(set(col)))) + if len(me_default) != nevt or not all(v > 0 for v in me_default): + problems.append('|M|^2 with the default card not > 0: %s' % list(me_default[:4])) check('zero_me', not problems, '; '.join(problems) or - '%d events: |M|^2 = 0, helicity -1, colour %s (default card: %.3e)' - % (npoints, sorted(set(col)), me_default[0])) + '%d events: |M|^2 = 0, helicity -1, %s (default card: %.3e)' + % (nevt, ', '.join(seen), me_default[0])) # --------------------------------------------------------------------------- # 3. a non-blocking stream # --------------------------------------------------------------------------- -def check_nonblocking_stream(repo, backend, npoints=16384, ncalls=4): +def check_nonblocking_stream(backend, repo, ncalls=4): section('g g > t t~ g on a non-blocking stream (%s)' % backend) + if backend not in GPU_BACKENDS: + check('nonblocking_stream', True, 'skipped: %s is not a GPU backend' % backend) + return name = 'ggttxg' - (sub,) = output(repo, name, 'sm', 'g g > t t~ g', (backend,)) + (proc_dir,) = output(repo, name, 'sm', 'g g > t t~ g', (backend,)) + # ncalls calls (iterations) of 64 x 256 = 16384 events, the RAMBO points of check_sa.exe + # (the same in both runs: one seed per iteration) + batch = [64, 256, ncalls] + nevt = 64 * 256 * ncalls rng = random.Random(3) - gpu = Umami(backend, library(name, sub, backend), card(name)) - points = [rambo(gpu.masses[2:], 1000., rng) for _ in range(npoints)] - rnd = [rng.random() for _ in range(npoints)] - outputs = (OUT_MATRIX_ELEMENT, OUT_HELICITY_INDEX) - stream = gpu.memory.nonblocking_stream() - # the first call on the non-blocking stream: it also runs the helicity filtering - calls = [gpu(points, random_helicity=rnd, outputs=outputs, stream=stream)] - reference = gpu(points, random_helicity=rnd, outputs=outputs) # default stream - calls += [gpu(points, random_helicity=rnd, outputs=outputs, stream=stream) - for _ in range(ncalls - 1)] - gpu.memory.stream_destroy(stream) - gpu.free() - me_ref, hel_ref = reference + rnd = write(name + '_rnd.bin', 'd', [rng.random() for _ in range(nevt)]) + me, hel = {}, {} + for tag, extra in (('default', []), ('nonblocking', ['--nonblocking-stream'])): + prefix = '%s_%s' % (name, tag) + check_sa(proc_dir, backend, extra + ['--energy', 1000, '--rnd-hel', rnd, + '--dump-me', dump(prefix, 'me'), + '--dump-hel', dump(prefix, 'hel')] + batch, + prefix + '.log') + me[tag], hel[tag] = read(dump(prefix, 'me'), 'd'), read(dump(prefix, 'hel'), 'i') problems = [] - for icall, (me, hel) in enumerate(calls): - bad_me = sum(1 for a, b in zip(me, me_ref) if not abs(a - b) <= 1e-12 * abs(b)) - bad_hel = sum(1 for a, b in zip(hel, hel_ref) if a != b) + if not len(me['default']) == len(me['nonblocking']) == nevt: + problems.append('%d and %d |M|^2 for %d events' % (len(me['default']), len(me['nonblocking']), nevt)) + for icall in range(ncalls): + events = slice(icall * nevt // ncalls, (icall + 1) * nevt // ncalls) + pairs = list(zip(me['nonblocking'][events], me['default'][events])) + bad_me = sum(1 for a, b in pairs if not abs(a - b) <= 1e-12 * abs(b)) + bad_hel = sum(1 for a, b in zip(hel['nonblocking'][events], hel['default'][events]) if a != b) if bad_me or bad_hel: problems.append('call %d: %d |M|^2 and %d helicities differ' % (icall, bad_me, bad_hel)) - if not all(v > 0 for v in me_ref) or min(hel_ref) < 0: + if not all(v > 0 for v in me['default']) or min(hel['default'], default=-1) < 0: problems.append('default stream: |M|^2 <= 0 or no helicity') check('nonblocking_stream', not problems, '; '.join(problems[:3]) or - '%d calls of %d events identical to the default stream' % (ncalls, npoints)) + '%d calls of %d events identical to the default stream' % (ncalls, nevt // ncalls)) def main(): parser = argparse.ArgumentParser(description=__doc__.split('\n')[0]) parser.add_argument('--repo', required=True) parser.add_argument('--backend', required=True) - parser.add_argument('--only', choices=CHECKS, - help='run this check here and write its result to NAME.json') + parser.add_argument('--only', choices=CHECKS, action='append', + help='run this check only (repeatable)') args = parser.parse_args() - if args.only: - step = globals()['check_' + args.only] - try: - step(args.repo, args.backend) - except Exception as error: - check(args.only, False, 'error: %s' % error) - json.dump(results, open(args.only + '.json', 'w')) - return 0 - if args.backend not in GPU_BACKENDS: - parser.error('the checks compare a GPU backend (%s) with the CPU one' - % ', '.join(GPU_BACKENDS)) + if args.backend == CPU_BACKEND: + parser.error('the checks compare a backend with %s: give another one' % CPU_BACKEND) start = time.time() - summary, failures = {}, [] - for name in CHECKS: - if os.path.exists(name + '.json'): - os.remove(name + '.json') - status = subprocess.call([sys.executable, os.path.abspath(__file__), '--repo', args.repo, - '--backend', args.backend, '--only', name]) + for name in args.only or CHECKS: try: - outcome = json.load(open(name + '.json')) - except (OSError, ValueError): - outcome = {name: (False, 'crashed (exit status %d)' % status)} - print('FAIL %s: %s' % (name, outcome[name][1]), flush=True) - for key, (ok, detail) in outcome.items(): - summary[key] = ('ok' if ok else 'FAILED') + ' (%s)' % detail - if not ok: - failures.append(key) + globals()['check_' + name](args.backend, args.repo) + except Exception as error: + check(name, False, 'error: %s' % error) + summary = {key: ('ok' if ok else 'FAILED') + ' (%s)' % detail + for key, (ok, detail) in results.items()} + failures = [key for key, (ok, _detail) in results.items() if not ok] summary['node'] = os.uname().nodename summary['gpu'] = os.environ.get('GPU_NAME', 'none') summary['backend'] = args.backend diff --git a/.github/gpu_ci/madmatrix_gpu.sh b/.github/gpu_ci/madmatrix_gpu.sh index cc6d297d5..cf0ed62ef 100755 --- a/.github/gpu_ci/madmatrix_gpu.sh +++ b/.github/gpu_ci/madmatrix_gpu.sh @@ -1,6 +1,6 @@ #!/bin/bash -l # Runs on the self-hosted runner (GPU node): job madmatrix of gpu_runner_ci.yml. -# The GPU backend of madmatrix against its CPU backend, through umami: see +# The GPU backend of madmatrix against its CPU backend, with check_sa.exe: see # madmatrix_checks.py for what is checked. # Environment: as pp_ttx_mg7.sh (BACKEND, MODULES, GPU_ARCH, VENV, WORKDIR); no madspace # and no PDF set needed (standalone outputs only). diff --git a/.github/gpu_ci/umami_harness.py b/.github/gpu_ci/umami_harness.py deleted file mode 100644 index 2b772029a..000000000 --- a/.github/gpu_ci/umami_harness.py +++ /dev/null @@ -1,183 +0,0 @@ -"""Calling a madmatrix process library through umami from Python (ctypes), on the CPU or -on the GPU, for the checks of the GPU CI (interference_checks.py, madmatrix_checks.py). -""" - -import ctypes -import ctypes.util -import math -import subprocess - -GPU_BACKENDS = ('cuda', 'hip') - -# umami.h -META_PARTICLE_COUNT, META_MASSES = 1, 5 -IN_MOMENTA, IN_ALPHA_S, IN_FLAVOR_INDEX, IN_RANDOM_COLOR, IN_RANDOM_HELICITY = 0, 1, 2, 3, 4 -OUT_MATRIX_ELEMENT, OUT_COLOR_INDEX, OUT_HELICITY_INDEX, OUT_GPU_STREAM = 0, 2, 3, 5 - - -class Memory: - """Buffers for umami: host memory for a CPU library, device memory for a GPU one - (through the CUDA/HIP runtime the library itself is linked with).""" - - def __init__(self, backend, library): - self.gpu = backend in GPU_BACKENDS - self.hosts = [] # the host arrays umami reads on a CPU, alive until release - if not self.gpu: - return - name, prefix = ('cudart', 'cuda') if backend == 'cuda' else ('amdhip64', 'hip') - path = None - ldd = subprocess.run(['ldd', library], capture_output=True, text=True).stdout - for line in ldd.splitlines(): - if 'lib%s' % name in line and '=>' in line: - path = line.split('=>')[1].split()[0] - self.rt = ctypes.CDLL(path or ctypes.util.find_library(name) or 'lib%s.so' % name) - self.malloc = getattr(self.rt, prefix + 'Malloc') - self.malloc.argtypes = [ctypes.POINTER(ctypes.c_void_p), ctypes.c_size_t] - self.memcpy = getattr(self.rt, prefix + 'Memcpy') - self.memcpy.argtypes = [ctypes.c_void_p, ctypes.c_void_p, ctypes.c_size_t, ctypes.c_int] - self.free = getattr(self.rt, prefix + 'Free') - self.free.argtypes = [ctypes.c_void_p] - self.sync = getattr(self.rt, prefix + 'DeviceSynchronize') - self.stream_create = getattr(self.rt, prefix + 'StreamCreateWithFlags') - self.stream_create.argtypes = [ctypes.POINTER(ctypes.c_void_p), ctypes.c_uint] - self.stream_destroy = getattr(self.rt, prefix + 'StreamDestroy') - self.stream_destroy.argtypes = [ctypes.c_void_p] - self.pointers = [] - - def buffer(self, host): - """A buffer for umami holding (a copy of) the ctypes array host.""" - if not self.gpu: - self.hosts.append(host) - return ctypes.addressof(host) - pointer = ctypes.c_void_p() - assert self.malloc(ctypes.byref(pointer), ctypes.sizeof(host)) == 0, 'device malloc' - # cudaMemcpyHostToDevice == hipMemcpyHostToDevice == 1 - assert self.memcpy(pointer, ctypes.addressof(host), ctypes.sizeof(host), 1) == 0 - self.pointers.append(pointer) - return pointer.value - - def fetch(self, pointer, host): - """Copy a umami output back into the ctypes array host.""" - if not self.gpu: - return - assert self.sync() == 0, 'device synchronize' - # cudaMemcpyDeviceToHost == hipMemcpyDeviceToHost == 2 - assert self.memcpy(ctypes.addressof(host), pointer, ctypes.sizeof(host), 2) == 0 - - def release(self): - self.hosts = [] - if self.gpu: - for pointer in self.pointers: - self.free(pointer) - self.pointers = [] - - def nonblocking_stream(self): - """A stream that does not synchronise with the default stream (as a torch side - stream): cudaStreamNonBlocking == hipStreamNonBlocking == 1.""" - stream = ctypes.c_void_p() - assert self.stream_create(ctypes.byref(stream), 1) == 0, 'stream create' - return stream.value - - -class Umami: - """One umami instance of a process library.""" - - def __init__(self, backend, library, param_card): - self.lib = ctypes.CDLL(library) - self.memory = Memory(backend, library) - self.handle = ctypes.c_void_p() - status = self.lib.umami_initialize(ctypes.byref(self.handle), param_card.encode()) - assert status == 0, 'umami_initialize status %d' % status - npar = ctypes.c_int() - self.lib.umami_get_meta(META_PARTICLE_COUNT, ctypes.byref(npar)) - self.npar = npar.value - masses = (ctypes.c_double * self.npar)() - assert self.lib.umami_get_meta(META_MASSES, masses) == 0, 'UMAMI_META_MASSES' - self.masses = list(masses) - self.lib.umami_matrix_element.argtypes = [ - ctypes.c_void_p, ctypes.c_size_t, ctypes.c_size_t, ctypes.c_size_t, - ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p), - ctypes.c_size_t, ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_void_p)] - self.lib.umami_set_parameter.argtypes = [ctypes.c_void_p, ctypes.c_char_p, - ctypes.c_double, ctypes.c_double] - - def set_parameter(self, name, value): - return self.lib.umami_set_parameter(self.handle, name.encode(), value, 0.) - - def __call__(self, points, alpha_s=None, flavor=None, random_helicity=None, - random_color=None, outputs=(OUT_MATRIX_ELEMENT,), stream=None): - """umami_matrix_element on these phase-space points (lists of npar 4-momenta) - and per-event inputs; the outputs as lists, in the order asked for.""" - count = len(points) - mom = (ctypes.c_double * (4 * self.npar * count))() - for ievt, point in enumerate(points): # momenta[count * (npar * mu + ipart) + ievt] - for ipart in range(self.npar): - for mu in range(4): - mom[count * (self.npar * mu + ipart) + ievt] = point[ipart][mu] - in_keys, inputs = [IN_MOMENTA], [self.memory.buffer(mom)] - for key, values, ctype in ((IN_ALPHA_S, alpha_s, ctypes.c_double), - (IN_FLAVOR_INDEX, flavor, ctypes.c_int), - (IN_RANDOM_HELICITY, random_helicity, ctypes.c_double), - (IN_RANDOM_COLOR, random_color, ctypes.c_double)): - if values is not None: - in_keys.append(key) - inputs.append(self.memory.buffer((ctype * count)(*values))) - hosts = [((ctypes.c_double if key == OUT_MATRIX_ELEMENT else ctypes.c_int) * count)() - for key in outputs] - out_keys = list(outputs) - out_pointers = [self.memory.buffer(host) for host in hosts] - if stream is not None: - out_keys.append(OUT_GPU_STREAM) - out_pointers.append(stream) - status = self.lib.umami_matrix_element( - self.handle, count, count, 0, - len(in_keys), (ctypes.c_int * len(in_keys))(*in_keys), - (ctypes.c_void_p * len(inputs))(*inputs), - len(out_keys), (ctypes.c_int * len(out_keys))(*out_keys), - (ctypes.c_void_p * len(out_pointers))(*out_pointers)) - assert status == 0, 'umami_matrix_element status %d' % status - for pointer, host in zip(out_pointers, hosts): - self.memory.fetch(pointer, host) - self.memory.release() - return [list(host) for host in hosts] - - def free(self): - self.lib.umami_free(self.handle) - - -def rambo(masses_out, roots, rng, masses_in=(0., 0.)): - """One phase-space point: two beams along z (masses masses_in) in their centre-of-mass - frame, then RAMBO with the masses masses_out.""" - n = len(masses_out) - q = [] - for _ in range(n): - c, f = 2 * rng.random() - 1, 2 * math.pi * rng.random() - e = -math.log(rng.random() * rng.random()) - s = math.sqrt(1 - c * c) - q.append([e, e * s * math.cos(f), e * s * math.sin(f), e * c]) - big_q = [sum(qi[k] for qi in q) for k in range(4)] - mass = math.sqrt(big_q[0] ** 2 - sum(big_q[k] ** 2 for k in (1, 2, 3))) - b = [-big_q[k] / mass for k in (1, 2, 3)] - g, x = big_q[0] / mass, roots / mass - a = 1 / (1 + g) - p = [] - for qi in q: - bq = sum(b[k] * qi[k + 1] for k in range(3)) - p.append([x * (g * qi[0] + bq)] + - [x * (qi[k + 1] + b[k] * qi[0] + a * bq * b[k]) for k in range(3)]) - # rescale the three-momenta so that the energies add up to roots with the masses - xi = 1. - for _ in range(100): - energies = [math.sqrt((xi * pi[0]) ** 2 + m * m) for pi, m in zip(p, masses_out)] - f = sum(energies) - roots - df = sum(xi * pi[0] ** 2 / e for pi, e in zip(p, energies)) - xi -= f / df - if abs(f) < 1e-12 * roots: - break - out = [] - for pi, m in zip(p, masses_out): - vec = [xi * v for v in pi[1:]] - out.append([math.sqrt(sum(v * v for v in vec) + m * m)] + vec) - m1, m2 = masses_in - pz = math.sqrt((roots ** 2 - (m1 + m2) ** 2) * (roots ** 2 - (m1 - m2) ** 2)) / (2 * roots) - return [[math.hypot(pz, m1), 0., 0., pz], [math.hypot(pz, m2), 0., 0., -pz]] + out diff --git a/.github/workflows/gpu_runner_ci.yml b/.github/workflows/gpu_runner_ci.yml index c848450ae..b588feead 100644 --- a/.github/workflows/gpu_runner_ci.yml +++ b/.github/workflows/gpu_runner_ci.yml @@ -251,7 +251,7 @@ jobs: if-no-files-found: warn madmatrix: - # the GPU backend of madmatrix against its CPU backend through umami: running flavour + # the GPU backend of madmatrix against its CPU backend with check_sa.exe: running flavour # couplings (MSSM), a subprocess with no good helicity, a non-blocking stream needs: [start_runner, build_madspace] if: needs.start_runner.outputs.attempt == github.run_attempt @@ -275,7 +275,7 @@ jobs: if: always() run: | mkdir -p results - (cd "$WORKDIR" && tar -cf - --ignore-failed-read summary.txt checks.log *.json *.mg5 *.log 2> /dev/null) | tar -xf - -C results || true + (cd "$WORKDIR" && tar -cf - --ignore-failed-read summary.txt checks.log *.mg5 *.log 2> /dev/null) | tar -xf - -C results || true find results -type f | sort - name: job summary if: always() diff --git a/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h b/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h index cdb6a0f11..9579602f7 100644 --- a/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h +++ b/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h @@ -52,6 +52,7 @@ #define gpuStream_t cudaStream_t #define gpuStreamCreate( pStream ) checkGpu( cudaStreamCreate( pStream ) ) +#define gpuStreamCreateNonBlocking( pStream ) checkGpu( cudaStreamCreateWithFlags( pStream, cudaStreamNonBlocking ) ) #define gpuStreamDestroy( stream ) checkGpu( cudaStreamDestroy( stream ) ) #define gpuMallocAsync( ptr, size, stream ) checkGpu( cudaMallocAsync( ptr, size, stream ) ) #define gpuFreeAsync( ptr, stream ) checkGpu( cudaFreeAsync( ptr, stream ) ) @@ -122,6 +123,7 @@ #define gpuStream_t hipStream_t #define gpuStreamCreate( pStream ) checkGpu( hipStreamCreate( pStream ) ) +#define gpuStreamCreateNonBlocking( pStream ) checkGpu( hipStreamCreateWithFlags( pStream, hipStreamNonBlocking ) ) #define gpuStreamDestroy( stream ) checkGpu( hipStreamDestroy( stream ) ) #define gpuMallocAsync( ptr, size, stream ) checkGpu( hipMallocAsync( ptr, size, stream ) ) #define gpuFreeAsync( ptr, stream ) checkGpu( hipFreeAsync( ptr, stream ) ) diff --git a/madgraph/iolibs/template_files/madmatrix/check_sa.cc b/madgraph/iolibs/template_files/madmatrix/check_sa.cc index aa0cf77c7..eea86a67e 100644 --- a/madgraph/iolibs/template_files/madmatrix/check_sa.cc +++ b/madgraph/iolibs/template_files/madmatrix/check_sa.cc @@ -50,6 +50,8 @@ #include #include #include +#include +#include #include #define STRINGIFY( s ) #s @@ -173,6 +175,42 @@ namespace return true; } + // Raw binary per-event values (--alphas , --flavors, --rnd-hel, --rnd-col): one + // value of type T per event, with no header, as array.tofile writes them in Python. + template + bool read_binary_values( const std::string& path, const char* what, std::vector& values ) + { + std::ifstream in( path, std::ios::binary ); + if( !in ) + { + std::cerr << "ERROR! cannot open " << what << " file '" << path << "'" << std::endl; + return false; + } + T value; + while( in.read( reinterpret_cast( &value ), sizeof( T ) ) ) + values.push_back( value ); + if( values.empty() || in.gcount() != 0 ) + { + std::cerr << "ERROR! " << what << " file '" << path << "' is empty or not a whole number of " + << sizeof( T ) << "-byte values" << std::endl; + return false; + } + return true; + } + + // Perf-mode options that add inputs or outputs to the umami call, or change how it runs + struct PerfIO + { + std::string alphas; // --alphas: a value, or a file of doubles + std::string flavorsFile; // --flavors: a file of 32-bit unsigned ints + std::string rndHelFile; // --rnd-hel: a file of doubles + std::string rndColFile; // --rnd-col: a file of doubles + std::string dumpHelFile; // --dump-hel: write the helicity indices + std::string dumpColFile; // --dump-col: write the color indices + std::vector> parameters; // --param = + bool nonblockingStream = false; // --nonblocking-stream (CUDA/HIP) + }; + int usage( const char* argv0, int ret = 1 ) { std::cout @@ -180,7 +218,9 @@ namespace << " " << argv0 << " [matrix] [-v|--verbose] []\n" << " " << argv0 << " perf [-v|--verbose] [-f|--flavor ] [--rambo-massless]" << " [-e|--events ] [--momenta ] [--dump-momenta ] [--dump-me ]" - << " [--energy ] [<#blocksPerGrid> <#threadsPerBlock>] <#iterations>\n" + << " [--energy ] [--alphas ] [--flavors ] [--rnd-hel ]" + << " [--rnd-col ] [--dump-hel ] [--dump-col ] [--param =]" + << " [--nonblocking-stream] [<#blocksPerGrid> <#threadsPerBlock>] <#iterations>\n" << " " << argv0 << " -p [opts] (legacy alias for `perf`)\n" << "\n" << "Subcommands:\n" @@ -205,6 +245,20 @@ namespace << " --dump-me (perf only) Write the matrix element of every event, raw doubles.\n" << " --energy (perf only) Ecms for RAMBO (default 1500 GeV).\n" << "\n" + << "Per-event inputs and outputs (perf only). A holds one raw value per event\n" + << "(no header): doubles, or 32-bit ints for the flavor and the indices. An input file\n" + << "needs a value for every event of the run, in the order of the momenta.\n" + << " --alphas alpha_s of every event (default: the param card).\n" + << " --flavors Flavor index of every event (instead of -f).\n" + << " --rnd-hel Random numbers of the helicity choice (umami otherwise uses 0.5\n" + << " and returns the helicity sum for an interference |M|^2).\n" + << " --rnd-col Random numbers of the color choice (otherwise 0.5).\n" + << " --dump-hel Write the chosen helicity index of every event (-1: none).\n" + << " --dump-col Write the chosen color index of every event (-1: none).\n" + << " --param = umami_set_parameter after umami_initialize (repeatable).\n" + << " --nonblocking-stream (CUDA/HIP) Give umami a non-blocking stream, which does not\n" + << " synchronise with the default stream.\n" + << "\n" << "Both modes read the parameters, alpha_s included, from ../../Cards/param_card.dat.\n" << "\n" << "perf-mode defaults if positional args are omitted:\n" @@ -312,55 +366,28 @@ namespace return slha.get_block_entry( "sminputs", 3, 1.180000e-01 ); } - // Run sigmaKin via UMAMI for `nevt` events and copy back the MEs. - // The momenta (UMAMI SoA layout), the per-event flavor and the per-event - // alpha_s buffers must be set by the caller. On GPU the buffers are device - // pointers and `hstMEs` receives the host-side copy; on CPU `umamiMEs` is the - // output buffer. + // Run sigmaKin via UMAMI for `nevt` events. The input and output buffers (the + // momenta in the UMAMI SoA layout, the per-event flavor and alpha_s, ...) are set + // by the caller, who also copies the outputs back: on GPU they are device pointers. bool run_umami( UmamiHandle handle, unsigned int nevt, mgOnGpu::TimerMap& timermap, double& wavetime, -#ifdef MGONGPUCPP_GPUIMPL - const DeviceBufferBase& devUmamiMomenta, - const DeviceBufferBase& devFlv, - const DeviceBufferBase& devAlphaS, - DeviceBufferBase& devUmamiMEs, - std::vector& hstMEs -#else - const std::vector& umamiMomenta, - const std::vector& flvVec, - const std::vector& alphasVec, - std::vector& umamiMEs -#endif - ) + const std::vector& in_keys, + const std::vector& inputs, + const std::vector& out_keys, + const std::vector& outputs ) { - constexpr unsigned int UmamiInKeyNum = 3; timermap.start( "3a SigmaKin" ); - UmamiInputKey in_keys[UmamiInKeyNum] = { UMAMI_IN_MOMENTA, UMAMI_IN_FLAVOR_INDEX, UMAMI_IN_ALPHA_S }; - UmamiOutputKey out_keys[1] = { UMAMI_OUT_MATRIX_ELEMENT }; -#ifdef MGONGPUCPP_GPUIMPL - const void* inputs[UmamiInKeyNum] = { devUmamiMomenta.data(), devFlv.data(), devAlphaS.data() }; - void* outputs[1] = { devUmamiMEs.data() }; -#else - const void* inputs[UmamiInKeyNum] = { umamiMomenta.data(), flvVec.data(), alphasVec.data() }; - void* outputs[1] = { umamiMEs.data() }; -#endif UmamiStatus st = umami_matrix_element( - handle, nevt, nevt, 0, UmamiInKeyNum, in_keys, inputs, 1, out_keys, outputs ); + handle, nevt, nevt, 0, in_keys.size(), in_keys.data(), inputs.data(), out_keys.size(), out_keys.data(), outputs.data() ); wavetime += timermap.stop(); if( st != UMAMI_SUCCESS ) { std::cerr << "ERROR! umami_matrix_element failed (status=" << st << ")" << std::endl; return false; } - -#ifdef MGONGPUCPP_GPUIMPL - timermap.start( "3b CpDTHmes" ); - gpuMemcpy( hstMEs.data(), devUmamiMEs.data(), nevt * sizeof( double ), gpuMemcpyDeviceToHost ); - wavetime += timermap.stop(); -#endif return true; } @@ -803,7 +830,8 @@ namespace const std::string& lheFile = "", const std::string& momentaFile = "", const std::string& dumpMomentaFile = "", - const std::string& dumpMEFile = "" ) + const std::string& dumpMEFile = "", + const PerfIO& io = PerfIO() ) { const unsigned int nevt = gpublocks * gputhreads; @@ -828,28 +856,62 @@ namespace } const bool eventsFromFile = !lheEvents.empty(); - std::ofstream dumpMomenta, dumpMEs; - if( !dumpMomentaFile.empty() ) + std::ofstream dumpMomenta, dumpMEs, dumpHel, dumpCol; + for( auto [file, path, what] : { std::make_tuple( &dumpMomenta, &dumpMomentaFile, "momenta" ), + std::make_tuple( &dumpMEs, &dumpMEFile, "matrix element" ), + std::make_tuple( &dumpHel, &io.dumpHelFile, "helicity" ), + std::make_tuple( &dumpCol, &io.dumpColFile, "color" ) } ) { - dumpMomenta.open( dumpMomentaFile, std::ios::binary ); - if( !dumpMomenta ) + if( path->empty() ) continue; + file->open( *path, std::ios::binary ); + if( !*file ) { - std::cerr << "ERROR! cannot write momenta file '" << dumpMomentaFile << "'" << std::endl; + std::cerr << "ERROR! cannot write " << what << " file '" << *path << "'" << std::endl; return 2; } } - if( !dumpMEFile.empty() ) + + // Per-event inputs from files: one value for each event of the run, in the order + // of the momenta (the padding events of the last batch repeat its last real event) + double alphasValue = param_card_alpha_s(); + std::vector alphasIn, rndHelIn, rndColIn; + std::vector flavorsIn; + if( is_float( io.alphas.c_str() ) ) + alphasValue = atof( io.alphas.c_str() ); + else if( !io.alphas.empty() && !read_binary_values( io.alphas, "alpha_s", alphasIn ) ) + return 2; + if( !io.flavorsFile.empty() && !read_binary_values( io.flavorsFile, "flavor", flavorsIn ) ) return 2; + if( !io.rndHelFile.empty() && !read_binary_values( io.rndHelFile, "helicity random number", rndHelIn ) ) return 2; + if( !io.rndColFile.empty() && !read_binary_values( io.rndColFile, "color random number", rndColIn ) ) return 2; + const std::size_t nevtRun = eventsFromFile ? lheEvents.size() : (std::size_t)niter * nevt; + for( auto [size, what] : { std::make_pair( alphasIn.size(), "alpha_s" ), + std::make_pair( flavorsIn.size(), "flavor" ), + std::make_pair( rndHelIn.size(), "helicity random number" ), + std::make_pair( rndColIn.size(), "color random number" ) } ) { - dumpMEs.open( dumpMEFile, std::ios::binary ); - if( !dumpMEs ) + if( size != 0 && size < nevtRun ) { - std::cerr << "ERROR! cannot write matrix element file '" << dumpMEFile << "'" << std::endl; + std::cerr << "ERROR! the " << what << " file has " << size << " values for " << nevtRun << " events" << std::endl; + return 2; + } + } + for( unsigned int iflav : flavorsIn ) + { + if( iflav >= CPPProcess::nmaxflavor ) + { + std::cerr << "ERROR! flavor index " << iflav << " in the flavor file is out of range [0, " + << CPPProcess::nmaxflavor << ")" << std::endl; return 2; } } mgOnGpu::TimerMap timermap; + // The per-event inputs and outputs of one batch (host side) + std::vector flvVec( nevt, flavorID ); + std::vector alphasVec( nevt, alphasValue ); + std::vector rndHelVec( nevt, 0.5 ), rndColVec( nevt, 0.5 ); + std::vector helVec( nevt ), colVec( nevt ); #ifdef MGONGPUCPP_GPUIMPL timermap.start( "00 GpuInit" ); GpuRuntime gpuRuntime( false ); @@ -863,21 +925,70 @@ namespace DeviceBufferBase devUmamiMomenta( (std::size_t)4 * CPPProcess::npar * nevt ); DeviceBufferBase devUmamiMEs( nevt ); DeviceBufferBase devFlv( nevt ); - std::vector flvVec( nevt, flavorID ); + DeviceBufferBase devAlphaS( nevt ); + DeviceBufferBase devRndHel( nevt ), devRndCol( nevt ); + DeviceBufferBase devHel( nevt ), devCol( nevt ); std::vector hstUmamiMEs( nevt ); - // perf-mode runs a single flavor, so the device-side flavor buffer is filled once. + // the inputs that are the same for every event are copied to the device once gpuMemcpy( devFlv.data(), flvVec.data(), nevt * sizeof( unsigned int ), gpuMemcpyHostToDevice ); - std::vector alphasVec( nevt, param_card_alpha_s() ); - DeviceBufferBase devAlphaS( nevt ); gpuMemcpy( devAlphaS.data(), alphasVec.data(), nevt * sizeof( double ), gpuMemcpyHostToDevice ); + void* umamiMomentaPtr = devUmamiMomenta.data(); + void* umamiMEsPtr = devUmamiMEs.data(); + unsigned int* flvPtr = devFlv.data(); + double* alphasPtr = devAlphaS.data(); + double* rndHelPtr = devRndHel.data(); + double* rndColPtr = devRndCol.data(); + int* helPtr = devHel.data(); + int* colPtr = devCol.data(); + gpuStream_t stream = nullptr; // nullptr: the default stream + if( io.nonblockingStream ) gpuStreamCreateNonBlocking( &stream ); #else HostBufferRndNumMomenta hstRndmom( nevt ); HostBufferMomenta hstMomenta( nevt ); HostBufferWeights hstWeights( nevt ); std::vector umamiMomenta( (std::size_t)4 * CPPProcess::npar * nevt ); std::vector umamiMEs( nevt ); - std::vector flvVec( nevt, flavorID ); - std::vector alphasVec( nevt, param_card_alpha_s() ); + void* umamiMomentaPtr = umamiMomenta.data(); + void* umamiMEsPtr = umamiMEs.data(); + unsigned int* flvPtr = flvVec.data(); + double* alphasPtr = alphasVec.data(); + double* rndHelPtr = rndHelVec.data(); + double* rndColPtr = rndColVec.data(); + int* helPtr = helVec.data(); + int* colPtr = colVec.data(); +#endif + + // The keys and buffers of the umami call, the same for every batch + std::vector in_keys = { UMAMI_IN_MOMENTA, UMAMI_IN_FLAVOR_INDEX, UMAMI_IN_ALPHA_S }; + std::vector inputs = { umamiMomentaPtr, flvPtr, alphasPtr }; + std::vector out_keys = { UMAMI_OUT_MATRIX_ELEMENT }; + std::vector outputs = { umamiMEsPtr }; + if( !rndHelIn.empty() ) + { + in_keys.push_back( UMAMI_IN_RANDOM_HELICITY ); + inputs.push_back( rndHelPtr ); + } + if( !rndColIn.empty() ) + { + in_keys.push_back( UMAMI_IN_RANDOM_COLOR ); + inputs.push_back( rndColPtr ); + } + if( dumpHel.is_open() ) + { + out_keys.push_back( UMAMI_OUT_HELICITY_INDEX ); + outputs.push_back( helPtr ); + } + if( dumpCol.is_open() ) + { + out_keys.push_back( UMAMI_OUT_COLOR_INDEX ); + outputs.push_back( colPtr ); + } +#ifdef MGONGPUCPP_GPUIMPL + if( stream != nullptr ) + { + out_keys.push_back( UMAMI_OUT_GPU_STREAM ); + outputs.push_back( stream ); + } #endif std::unique_ptr prnk( @@ -889,6 +1000,15 @@ namespace std::cerr << "ERROR! umami_initialize failed" << std::endl; return 2; } + for( const auto& [name, value] : io.parameters ) + { + if( umami_set_parameter( umami_handle, name.c_str(), value, 0. ) != UMAMI_SUCCESS ) + { + std::cerr << "ERROR! umami_set_parameter( " << name << ", " << value << " ) failed" << std::endl; + umami_free( umami_handle ); + return 2; + } + } // Retrieve masses std::vector masses; @@ -1010,20 +1130,46 @@ namespace #endif rambtime += timermap.stop(); - double wavetime = 0; - if( !run_umami( umami_handle, nevt, timermap, wavetime, + // The per-event inputs of this batch read from files + const std::size_t base = (std::size_t)iiter * nevt; + for( auto [batch, values] : { std::make_pair( &alphasVec, &alphasIn ), + std::make_pair( &rndHelVec, &rndHelIn ), + std::make_pair( &rndColVec, &rndColIn ) } ) + { + if( values->empty() ) continue; + for( unsigned int ievt = 0; ievt < nevt; ++ievt ) + ( *batch )[ievt] = ( *values )[base + std::min( ievt, nreal - 1 )]; + } + if( !flavorsIn.empty() ) + { + for( unsigned int ievt = 0; ievt < nevt; ++ievt ) + flvVec[ievt] = flavorsIn[base + std::min( ievt, nreal - 1 )]; + } #ifdef MGONGPUCPP_GPUIMPL - devUmamiMomenta, devFlv, devAlphaS, devUmamiMEs, hstUmamiMEs -#else - umamiMomenta, flvVec, alphasVec, umamiMEs + if( !alphasIn.empty() ) gpuMemcpy( alphasPtr, alphasVec.data(), nevt * sizeof( double ), gpuMemcpyHostToDevice ); + if( !rndHelIn.empty() ) gpuMemcpy( rndHelPtr, rndHelVec.data(), nevt * sizeof( double ), gpuMemcpyHostToDevice ); + if( !rndColIn.empty() ) gpuMemcpy( rndColPtr, rndColVec.data(), nevt * sizeof( double ), gpuMemcpyHostToDevice ); + if( !flavorsIn.empty() ) gpuMemcpy( flvPtr, flvVec.data(), nevt * sizeof( unsigned int ), gpuMemcpyHostToDevice ); + // A non-blocking stream does not wait for the inputs copied and transposed on the + // default stream: finish them first, as any caller on such a stream must + if( stream != nullptr ) checkGpu( gpuDeviceSynchronize() ); #endif - ) ) + + double wavetime = 0; + if( !run_umami( umami_handle, nevt, timermap, wavetime, in_keys, inputs, out_keys, outputs ) ) { umami_free( umami_handle ); return 3; } #ifdef MGONGPUCPP_GPUIMPL + timermap.start( "3b CpDTHmes" ); + // nor do the copies below wait for the work on a non-blocking stream + if( stream != nullptr ) gpuStreamSynchronize( stream ); + gpuMemcpy( hstUmamiMEs.data(), umamiMEsPtr, nevt * sizeof( double ), gpuMemcpyDeviceToHost ); + if( dumpHel.is_open() ) gpuMemcpy( helVec.data(), helPtr, nevt * sizeof( int ), gpuMemcpyDeviceToHost ); + if( dumpCol.is_open() ) gpuMemcpy( colVec.data(), colPtr, nevt * sizeof( int ), gpuMemcpyDeviceToHost ); + wavetime += timermap.stop(); if( verbose || dumpMomenta.is_open() ) { timermap.start( "3c CpDTHmom" ); @@ -1047,6 +1193,10 @@ namespace } if( dumpMEs.is_open() ) dumpMEs.write( reinterpret_cast( mes ), (std::streamsize)nreal * sizeof( double ) ); + if( dumpHel.is_open() ) + dumpHel.write( reinterpret_cast( helVec.data() ), (std::streamsize)nreal * sizeof( int ) ); + if( dumpCol.is_open() ) + dumpCol.write( reinterpret_cast( colVec.data() ), (std::streamsize)nreal * sizeof( int ) ); timermap.start( "4@ UpdtStat" ); for( unsigned int ievt = 0; ievt < nreal; ++ievt ) @@ -1105,7 +1255,8 @@ namespace std::cout << "NumBlocksPerGrid = " << gpublocks << std::endl << "NumThreadsPerBlock = " << gputhreads << std::endl << "NumIterations = " << niter << std::endl - << "FlavorIndex = " << flavorID << " / " << CPPProcess::nmaxflavor << std::endl + << "FlavorIndex = " << ( flavorsIn.empty() ? std::to_string( flavorID ) : "per event" ) + << " / " << CPPProcess::nmaxflavor << std::endl << std::string( SEP79, '-' ) << std::endl << "NaN/abnormal MEs = " << nevtABN << std::endl << "Zero MEs = " << nevtZERO << std::endl @@ -1134,10 +1285,16 @@ namespace std::cout << std::string( SEP79, '*' ) << std::endl; umami_free( umami_handle ); - if( ( dumpMomenta.is_open() && !dumpMomenta.flush() ) || ( dumpMEs.is_open() && !dumpMEs.flush() ) ) +#ifdef MGONGPUCPP_GPUIMPL + if( stream != nullptr ) gpuStreamDestroy( stream ); +#endif + for( std::ofstream* file : { &dumpMomenta, &dumpMEs, &dumpHel, &dumpCol } ) { - std::cerr << "ERROR! failed to write the dump files" << std::endl; - return 2; + if( file->is_open() && !file->flush() ) + { + std::cerr << "ERROR! failed to write the dump files" << std::endl; + return 2; + } } return 0; } @@ -1161,6 +1318,7 @@ int main( int argc, char** argv ) std::string dumpMomentaFile; // --dump-momenta: write the momenta of every event (perf mode only) std::string dumpMEFile; // --dump-me: write the matrix element of every event (perf mode only) double perfEnergy = -1.; // --energy: Ecms for the perf-mode RAMBO + PerfIO io; // the per-event inputs and outputs of perf mode (--alphas, ...) // Optional leading subcommand (no leading dash). int firstArg = 1; @@ -1184,7 +1342,6 @@ int main( int argc, char** argv ) flavorID = strtoul( argv[++argn], nullptr, 0 ); else if( arg == "--rambo-massless" ) { - std::string r = argv[++argn]; ramboType = RAMBO_MASSLESS; ramboTypeSet = true; } @@ -1210,6 +1367,27 @@ int main( int argc, char** argv ) } else if( arg == "--energy" && argn + 1 < argc && is_float( argv[argn + 1] ) ) perfEnergy = atof( argv[++argn] ); + else if( ( arg == "--alphas" || arg == "--flavors" || arg == "--rnd-hel" || arg == "--rnd-col" || + arg == "--dump-hel" || arg == "--dump-col" ) && argn + 1 < argc ) + { + std::string& value = ( arg == "--alphas" ? io.alphas : arg == "--flavors" ? io.flavorsFile : + arg == "--rnd-hel" ? io.rndHelFile : arg == "--rnd-col" ? io.rndColFile : + arg == "--dump-hel" ? io.dumpHelFile : io.dumpColFile ); + value = argv[++argn]; + mode = MODE_PERF; + } + else if( arg == "--param" && argn + 1 < argc && strchr( argv[argn + 1], '=' ) != nullptr ) + { + const std::string param = argv[++argn]; + const std::size_t eq = param.find( '=' ); + io.parameters.emplace_back( param.substr( 0, eq ), atof( param.c_str() + eq + 1 ) ); + mode = MODE_PERF; + } + else if( arg == "--nonblocking-stream" ) + { + io.nonblockingStream = true; + mode = MODE_PERF; + } else if( is_number( argv[argn] ) && nnum < 3 ) { numvec[nnum++] = strtoul( argv[argn], nullptr, 0 ); @@ -1266,6 +1444,13 @@ int main( int argc, char** argv ) return usage( argv[0] ); } if( perfEnergy > 0 ) kEnergy = perfEnergy; +#ifndef MGONGPUCPP_GPUIMPL + if( io.nonblockingStream ) + { + std::cerr << "ERROR: --nonblocking-stream needs a CUDA or HIP build." << std::endl; + return 1; + } +#endif if( flavorID >= CPPProcess::nmaxflavor ) { @@ -1275,5 +1460,5 @@ int main( int argc, char** argv ) } return run_perf_mode( verbose, gpublocks, gputhreads, niter, flavorID, ramboType, lheFile, - momentaFile, dumpMomentaFile, dumpMEFile ); + momentaFile, dumpMomentaFile, dumpMEFile, io ); }