diff --git a/.github/gpu_ci/README.md b/.github/gpu_ci/README.md index 91857ded40..d56dde6b90 100644 --- a/.github/gpu_ci/README.md +++ b/.github/gpu_ci/README.md @@ -23,7 +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 | +| `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 3656ed907f..3f8d491734 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,14 +30,12 @@ """ import argparse -import ctypes -import ctypes.util +import array import glob import gzip import json import math import os -import random import re import subprocess import sys @@ -136,149 +135,56 @@ 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') - 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): @@ -306,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 new file mode 100644 index 0000000000..797d9f78fc --- /dev/null +++ b/.github/gpu_ci/madmatrix_checks.py @@ -0,0 +1,301 @@ +#!/usr/bin/env python3 +"""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 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 (--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``, 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 math +import os +import random +import re +import shutil +import subprocess +import sys +import time + +GPU_BACKENDS = ('cuda', 'hip') +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 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, 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) + 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 proc_dir in proc_dirs: + for backend in backends: + 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=proc_dir, log=log) != 0: + raise RuntimeError('make failed, see %s' % log) + shutil.copy(os.path.join(proc_dir, 'check_sa.exe'), + os.path.join(proc_dir, 'check_sa_%s.exe' % backend)) + return proc_dirs + + +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 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(backend, repo, nevt=1024): + section('MSSM p p > go go, running flavour couplings: %s against %s' % (backend, CPU_BACKEND)) + name = 'mssm_gogo' + proc_dirs = output(repo, name, 'MSSM_SLHA2', 'p p > go go', (CPU_BACKEND, backend)) + rng = random.Random(11) + problems, details = [], [] + 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) + 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'])] + # 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), 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) + # |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]) + check('mssm_gpu_vs_cpu', not problems, '; '.join(problems[:3]) or '; '.join(details)) + + +# --------------------------------------------------------------------------- +# 2. a subprocess with no good helicity +# --------------------------------------------------------------------------- + +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' + (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) + rng = random.Random(5) + 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, %s (default card: %.3e)' + % (nevt, ', '.join(seen), me_default[0])) + + +# --------------------------------------------------------------------------- +# 3. a non-blocking stream +# --------------------------------------------------------------------------- + +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' + (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) + 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 = [] + 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['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, 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, action='append', + help='run this check only (repeatable)') + args = parser.parse_args() + if args.backend == CPU_BACKEND: + parser.error('the checks compare a backend with %s: give another one' % CPU_BACKEND) + start = time.time() + for name in args.only or CHECKS: + try: + 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 + 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 0000000000..cf0ed62ef9 --- /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, 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). +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/workflows/gpu_runner_ci.yml b/.github/workflows/gpu_runner_ci.yml index fb2a5490c2..b588feeadd 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 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 + 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 *.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 be4b731919..214075cf24 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,9 +267,9 @@ 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). - constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); + // 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 * 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++ ) @@ -276,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/GpuAbstraction.h b/madgraph/iolibs/template_files/madmatrix/backend/gpu/GpuAbstraction.h index b17ea25c40..9579602f77 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 @@ -51,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 ) ) @@ -98,6 +100,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 @@ -120,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/backend/gpu/MatrixElementKernels.cc b/madgraph/iolibs/template_files/madmatrix/backend/gpu/MatrixElementKernels.cc index 419241e5c5..9600010c8c 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/gpu/SigmaKin.cc b/madgraph/iolibs/template_files/madmatrix/backend/gpu/SigmaKin.cc index 4c1ad38fa3..e1c84c8ef3 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 d85651573c..5c0402444b 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 2d78abac85..969e0a8b23 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,9 +268,9 @@ 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). - constexpr int ndpfbuf = ( nDPF > 0 ? nDPF * nMF * CD_ACCESS::flv_stride : 1 ); + // 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 * 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 +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(); } diff --git a/madgraph/iolibs/template_files/madmatrix/check_sa.cc b/madgraph/iolibs/template_files/madmatrix/check_sa.cc index aa0cf77c7c..eea86a67eb 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 ); } diff --git a/madmatrix/model_handling.py b/madmatrix/model_handling.py index dfaf00419b..7282547d5a 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