From 77a333f122d9b517376280ae30ce43f61ea28a09 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Thu, 20 Aug 2026 08:11:43 +0200 Subject: [PATCH 1/3] madmatrix: refuse a squared-order component it cannot compute `output standalone_mg7` accepted `^2` syntax and wrote a matrix element with no squared-order machinery in it at all -- `grep split_orders madmatrix/` is still 0 hits, and so is `grep split_orders madgraph/iolibs/export_cpp.py`. That is not wrong everywhere, which is why it went unnoticed: this backend contracts one jamp vector per helicity with the color matrix once, so the |M|^2 it returns is the sum over every squared-order component, and that IS the requested number whenever the constraint keeps all of them. It is a wrong number when the constraint drops one. MG5 resolves what it can at generation -- `u u~ > t t~ QED^2==4` prunes to 2 diagrams and NSQAMPSO=1 -- but a component that needs two amplitudes cannot be reached that way. `u u~ > t t~ QED^2==2` keeps all three diagrams and asks for the QCD-EW cross term alone: fortran, masked -3.4453535703428253e-18 madmatrix 6.1741265892435804e-01 = the whole QCD^2 + int + EW^2 silently, from generated code differing from the unconstrained process by one comment line. `u u~ > u u~ QED^2==2` is the same failure with the sign visible: -5.5828747824035893e-02 asked for, +2.7756451 returned. WHAT IS SUPPORTED, and is byte-identical to before: NSQAMPSO==1, the constraint resolved at generation. `p p > j j QCD^2==4` (all 8 subprocesses build and run under FPTYPE=d), `u u~ > t t~ QED^2==4`, `QED^2==0`. NSQAMPSO>1 with every component kept, e.g. `u u~ > t t~ QED^2<=4`: the total is what was asked for. madmatrix 6.1741265892435804e-01 against the Fortran ANS(0) 0.61741265892435848, ~3 ulp. It now says out loud that the components behind that total are not available here, rather than leaving the user to assume they were applied. WHAT IS REFUSED: a constraint that drops a component. The message names the components, which are kept and which dropped, and points at `output standalone`. WHY REFUSE RATHER THAN MASK. Applying the mask needs the jamps split by amplitude order and paired in the color sum -- the Fortran GET_MATRIX contract, JAMP(NCOLOR,NAMPSO) -> RES(NSQAMPSO). Masking alone would keep the scalar signature, so it is not the return type that blocks it; it is that all three color sums (color_sum_cpu, the CUDA kernel, the BLAS path) take a single jamp vector per helicity, as do the color folding, the C-parity de-duplication and the good-helicity scan reading their output. The per-component vector on top of that would also need an ME buffer with a component index, which the SIMD lane layout and every consumer of allMEs assume away. madspace needs nothing and got nothing: MatrixElementOutput has no component kind and the umami metadata keys carry no squared-order count, so there is no ABI slot and no consumer for a component vector. The gate covers `output mg7` too, so madspace can never be handed a process whose mask it cannot honour. The predicate is factored out of set_chosen_SO_index as export_v4.chosen_squared_orders so the two cannot drift; the Fortran matrix.f, CHOSEN_SO_CONFIGS included, is byte-identical. NOT DONE, deliberately: the other C++ exporters share the gap. `output standalone_cpp` of `u u~ > t t~ QED^2==2` still returns 6.1741266e-01 (measured). Gating ProcessExporterCPP would cover standalone_cpp, mg7_v5, matchbox_cpp and pythia8 in one line, but matchbox_cpp has IOTest goldens and a golden regeneration is in flight, so that is left as its own change. Co-Authored-By: Claude Opus 5 --- madgraph/iolibs/export_v4.py | 111 +++++++++++++++++++++-------- madmatrix/output.py | 63 ++++++++++++++++ tests/acceptance_tests/test_cmd.py | 49 +++++++++++++ 3 files changed, 195 insertions(+), 28 deletions(-) diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index 49a59dff11..b72da22ce2 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -221,6 +221,86 @@ def export_helas(self, HELAS_PATH): raise Exception("V4 model not supported by this type of exporter. Please use UFO model") return +#=============================================================================== +# Squared split orders: which components the user's constraint keeps +#=============================================================================== +def chosen_squared_orders(process, squared_orders): + """Which entries of squared_orders the user's '^2' constraint keeps. + + squared_orders is the list of squared split-order tuples actually present + in the matrix element, as get_split_orders_mapping returns it; the process + carries the constraint. A False entry is a component that contributes to + the amplitude but that the user asked *not* to have summed into the total, + so a backend which cannot mask it does not compute what was asked for. + + This is the boolean form of what set_chosen_SO_index writes as the Fortran + CHOSEN_SO_CONFIGS DATA statement, kept here so that a backend can ask the + question without parsing that string. + """ + + user_squared_orders = process.get('squared_orders') + split_orders = process.get('split_orders') + + if len(user_squared_orders)==0: + return [True]*len(squared_orders) + + res = [] + for sqsos in squared_orders: + is_a_match = True + for user_sqso, value in user_squared_orders.items(): + if user_sqso == 'WEIGHTED' : + logger.debug('WEIGHTED^2%s%s encoutered. Please check behavior for' + \ + 'https://bazaar.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/revision/613', \ + (process.get_squared_order_type(user_sqso), sqsos[split_orders.index(user_sqso)])) + if user_sqso not in split_orders: + is_a_match = False + elif (process.get_squared_order_type(user_sqso) =='==' and \ + value!=sqsos[split_orders.index(user_sqso)]) or \ + (process.get_squared_order_type(user_sqso) in ['<=','='] and \ + value' and \ + value>=sqsos[split_orders.index(user_sqso)]): + is_a_match = False + break + res.append(is_a_match) + + return res + +def split_orders_not_supported_msg(format, process, split_orders, + squared_orders, chosen): + """Error text for a backend asked for a strict subset of squared orders. + + The backend sums every squared-order component into one |M|^2. That is the + number the user asked for as long as the constraint keeps all of them -- + which is why this is only an error when at least one is dropped. Naming the + components, and which of them survive, is the point: the difference between + what would be returned and what was asked for is otherwise invisible. + """ + + def name(orders): + return ' '.join('%s=%d' % (order, value) + for order, value in zip(split_orders, orders)) or 'ALL_ORDERS' + + kept = [name(o) for o, keep in zip(squared_orders, chosen) if keep] + dropped = [name(o) for o, keep in zip(squared_orders, chosen) if not keep] + + return """The '%(format)s' output format does not support squared split orders. +Process: %(proc)s +Its matrix element has %(n)d squared-order components: %(all)s. +The constraint keeps %(kept)s and drops %(dropped)s, but this backend has no +squared-order mask: it would sum all %(n)d components and return that total, +which is not the requested contribution. +Use 'output standalone' for this process -- the Fortran standalone carries the +split-order machinery (SMATRIX_SPLITORDERS) and returns the masked total as +well as each component. Alternatively, constrain the process so that all its +squared-order components are kept.""" % { + 'format': format, + 'proc': process.nice_string().replace('Process: ', ''), + 'n': len(squared_orders), + 'all': ', '.join(name(o) for o in squared_orders), + 'kept': ', '.join(kept) if kept else 'nothing', + 'dropped': ', '.join(dropped)} + #=============================================================================== # ProcessExporterFortran #=============================================================================== @@ -2017,34 +2097,9 @@ def set_chosen_SO_index(self, process, squared_orders): finds what indices of the squared_orders list the user intends to pick. It returns this as a string of comma-separated successive '.true.' or '.false.' for each index.""" - - user_squared_orders = process.get('squared_orders') - split_orders = process.get('split_orders') - - if len(user_squared_orders)==0: - return ','.join(['.true.']*len(squared_orders)) - - res = [] - for sqsos in squared_orders: - is_a_match = True - for user_sqso, value in user_squared_orders.items(): - if user_sqso == 'WEIGHTED' : - logger.debug('WEIGHTED^2%s%s encoutered. Please check behavior for' + \ - 'https://bazaar.launchpad.net/~maddevelopers/mg5amcnlo/3.0.1/revision/613', \ - (process.get_squared_order_type(user_sqso), sqsos[split_orders.index(user_sqso)])) - if user_sqso not in split_orders: - is_a_match = False - elif (process.get_squared_order_type(user_sqso) =='==' and \ - value!=sqsos[split_orders.index(user_sqso)]) or \ - (process.get_squared_order_type(user_sqso) in ['<=','='] and \ - value' and \ - value>=sqsos[split_orders.index(user_sqso)]): - is_a_match = False - break - res.append('.true.' if is_a_match else '.false.') - - return ','.join(res) + + return ','.join('.true.' if keep else '.false.' for keep in + chosen_squared_orders(process, squared_orders)) def get_split_orders_lines(self, orders, array_name, n=5): """ Return the split orders definition as defined in the list orders and diff --git a/madmatrix/output.py b/madmatrix/output.py index c72c7f8759..b8cfa48784 100644 --- a/madmatrix/output.py +++ b/madmatrix/output.py @@ -26,6 +26,7 @@ import madgraph.iolibs.export_v4 as export_v4 import madgraph.iolibs.export_cpp as export_cpp import madgraph.various.misc as misc +from madgraph import InvalidCmd from . import launch_plugin @@ -54,6 +55,10 @@ class ProcessExporterMadMatrix(export_cpp.ProcessExporterMG7): # AV - keep OM's default for this plugin (using grouped_mode=False, "can decide to merge uu~ and u~u anyway") sa_symmetry = True + # The name this exporter is reached by on the 'output' line, for the error + # messages that have to name it back to the user. + format_name = 'mg7' + # The color sum can run on the (n-2)! Del Duca-Dixon-Maltoni basis for a # multi-gluon process, but a color flow still has to be picked among the # (n-1)! trace structures, so the trace basis is built alongside and the @@ -206,8 +211,64 @@ def write_p_makefiles(self, model): rendered = self.read_template_file(pjoin(self.madmatrix_templates, name)) % replace_dict open(pjoin(self.dir_path, 'SubProcesses', name), 'w').write(rendered) + def check_split_orders(self, matrix_element): + """Refuse a process whose squared-order constraint drops a component. + + This backend has no squared-order machinery at all: it builds one jamp + vector per helicity and contracts it with the color matrix once, so the + |M|^2 it returns is the sum over *every* squared-order component the + amplitude has. That is exactly the requested number whenever the + constraint keeps all of them -- which covers the single-component case + (MG5 has already resolved the constraint at generation, e.g. + `p p > j j QCD^2==4`) and the multi-component case that keeps + everything (e.g. `u u~ > t t~ QED^2<=4`). Both stay supported, and + untouched. + + It is not the requested number when the constraint drops a component, + which happens when that component cannot be reached by dropping + diagrams -- an interference term, typically: `u u~ > t t~ QED^2==2` + keeps all three diagrams and asks for the QCD-EW cross term alone, + while this backend would return the full QCD^2 + interference + EW^2 + total. There is no mask to apply here, so it says so instead. + + Supporting it would mean splitting the jamps by amplitude order and + pairing them in the color sum, the Fortran GET_MATRIX contract + (JAMP(NCOLOR,NAMPSO) -> RES(NSQAMPSO)). All three color sums here -- + color_sum_cpu, the CUDA kernel and the BLAS path -- take a single jamp + vector per helicity, as do the color folding, the C-parity de-duplication + and the good-helicity scan that read their output. Masking alone would at + least keep the scalar signature; it is the jamp splitting underneath that + this backend has no room for, which is why this refuses rather than + computing the requested contribution. + """ + + process = matrix_element.get('processes')[0] + split_orders = process.get('split_orders') + if not split_orders: + return + squared_orders, _amp_orders = matrix_element.get_split_orders_mapping() + chosen = export_v4.chosen_squared_orders(process, squared_orders) + if all(chosen): + # The total is what was asked for. The components it is made of are + # not available here though, and the Fortran standalone would hand + # them back through SMATRIX_SPLITORDERS, so say which of the two + # this output is rather than let the user assume the other. + if len(squared_orders) > 1: + logger.warning( + "%s returns the summed |M|^2 over all %d squared-order " + "components of '%s', which is what its constraint asks " + "for. The individual components are not available from " + "this backend -- use 'output standalone' for those.", + self.__class__.format_name, len(squared_orders), + process.nice_string().replace('Process: ', '')) + return + raise InvalidCmd(export_v4.split_orders_not_supported_msg( + self.__class__.format_name, process, split_orders, + squared_orders, chosen)) + # AV - add debug printouts (in addition to the default one from OM's tutorial) def generate_subprocess_directory(self, matrix_element, cpp_helas_call_writer, proc_number=None): + self.check_split_orders(matrix_element) # Propagate the --mask toggle to the helas call writer that emits the # guarded wavefunction/amplitude calls, and the output command line as # a whole for the --jamp_optim toggle of the color-flow optimisation. @@ -240,6 +301,8 @@ def modify_grouping(self, matrix_element): # so that when running `make` in a P* folder, it builds check_sa.exe as well as the process library (predicatable behaviour) class ProcessExporterMadMatrixStandalone(ProcessExporterMadMatrix): + format_name = 'standalone_mg7' + # Each P* directory links madmatrix_standalone.mk (which itself includes # madmatrix.mk) as its 'makefile'; both have to be rendered in SubProcesses/ p_makefiles = ProcessExporterMadMatrix.p_makefiles + ['madmatrix_standalone.mk'] diff --git a/tests/acceptance_tests/test_cmd.py b/tests/acceptance_tests/test_cmd.py index 5100f7c763..380ad0e68c 100755 --- a/tests/acceptance_tests/test_cmd.py +++ b/tests/acceptance_tests/test_cmd.py @@ -1334,6 +1334,55 @@ def get_values(output_format, check_exe, build_source=False): 'all matrix elements vanished for u u~ > j j') self._assert_me_lists_close(mg7, standalone, atol=1e-7) + def test_standalone_mg7_refuses_a_dropped_squared_order(self): + """standalone_mg7 must refuse a squared-order constraint it cannot apply. + + The madmatrix backend contracts one jamp vector per helicity with the + color matrix once, so the |M|^2 it returns is the sum over every + squared-order component of the amplitude. Whenever the user's '^2' + constraint keeps all of them that is the requested number, and those + cases stay supported; when it drops one, the backend has no mask and + would hand back the full total instead. + + ``u u~ > t t~ QED^2==2`` is the case that matters: the interference + cannot be reached by dropping diagrams, so all three survive and the + Fortran masks the sum down to the QCD-EW cross term alone (-3.4e-18 at + the standard RAMBO point) while madmatrix returned the whole + 0.617412658924358 -- silently, with generated code differing from the + unconstrained process by one comment line. It has to say so instead. + + The supported cases are pinned alongside it, since the fix must not + turn them into errors: ``QED^2<=4`` keeps all three components and + ``QED^2==4`` is resolved by MG5 at generation down to a single one. + """ + self.do('import model sm') + + # Dropped component: refuse, and name what was dropped. + self.do('generate u u~ > t t~ QED^2==2') + try: + self.do('output standalone_mg7 %s -f' % self.out_dir) + except InvalidCmd as error: + message = str(error) + else: + self.fail('standalone_mg7 accepted a squared-order constraint it ' + 'cannot apply') + self.assertIn('does not support squared split orders', message) + self.assertIn('QED=0', message) # the components it would have summed + self.assertIn('QED=4', message) + + # All components kept: the plain sum is the requested number. + for process in ('u u~ > t t~ QED^2<=4', # three components, all kept + 'u u~ > t t~ QED^2==4'): # resolved at generation + if os.path.isdir(self.out_dir): + shutil.rmtree(self.out_dir) + self.do('generate %s' % process) + self.do('output standalone_mg7 %s -f' % self.out_dir) + proc_root = pjoin(self.out_dir, 'SubProcesses') + self.assertTrue([d for d in os.listdir(proc_root) + if d.startswith('P') and + os.path.isdir(pjoin(proc_root, d))], + 'standalone_mg7 wrote no subprocess for %s' % process) + def test_standalone_cpp(self): """test that the scalar C++ standalone exporter is working From 893beccf3d55da688c422f7f563f71e5f1e340ef Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Thu, 20 Aug 2026 09:34:16 +0200 Subject: [PATCH 2/3] madmatrix: squared split orders, with the jamps split by amplitude order The previous commit refused a '^2' constraint that drops a squared-order component, because this backend had no way to compute one. It has one now, and the refusal goes with it: `u u~ > u u~ QED^2==2` returns the QCD-EW cross term -5.5828746494657258e-02 where it used to return the whole +2.7756451. THE SHAPE, from the Fortran GET_MATRIX contract JAMP(NCOLOR,NAMPSO) -> RES(NSQAMPSO). jamp_sv becomes njampso = ncolor*nampso long, one vector per amplitude split order, and each amplitude is written into the vector of the order it belongs to. The color sum pairs them. THE PAIR LOOP RUNS OVER ALL nampso*nampso ORDERED PAIRS, and that is the one thing in here that is easy to get wrong. The color contraction keeps the triangular matrix with its off-diagonal doubled, which for a single jamp vector is just the statement that the (i,j) and (j,i) terms are conjugates. For two DIFFERENT jamp vectors they are not, so what the triangle computes for one ordered pair (m,n) is NOT that pair's contribution -- only (m,n) and (n,m) together are, and their sum is right. Halving the loop the way the color indices are halved gives roughly twice the answer. Masking stays safe for the same reason it is symmetric: a squared order is the SUM of the two amplitude orders (sqSoIndex[m][n] == sqSoIndex[n][m]), so a pair and its transpose are always kept or dropped together and the masked total stays real. A DEDICATED TEMPLATE, color_sum_splitorders.cc, because the sum changes shape rather than gaining a hole. Everything else is a hole whose empty fill is the text that was there before -- process_class.inc, process_matrix.inc, process_function_definitions.inc, process_sigmaKin_function.inc -- so a process without split orders is byte-for-byte unchanged, checked on `p p > j j QCD^2==4` (8 subprocesses), `u u~ > t t~ QED^2==4` and plain `u u~ > t t~`. SWITCHED OFF WHILE THE JAMPS ARE SPLIT The jamp CSE/orbit optimisation. Its shared sub-expressions are partial sums of amplitudes, and two amplitudes of different orders must not land in the same one; the flows are written out one amplitude at a time instead. The BLAS color sum, which multiplies a single jamp vector by the color matrix. Color choice is NOT split: jampso_sv[icol] = sum_M jamp_M[icol] is the color flow of the whole amplitude, and deliberately unmasked -- it picks which flow to write out, it is not an observable. With one amplitude order it is exactly today's jamp2. GPU REFUSES AT COMPILE TIME. The device jamp buffers are sized ncolor, not njampso, and the backend is a make-time flag rather than an output-time one, so the refusal cannot live in the exporter: color_sum_splitorders.cc #errors under MGONGPUCPP_GPUIMPL. The exporter says so at generation as well, so a GPU build is not the first the user hears of it. VERIFIED by building and running, at the standard RAMBO point with FPTYPE=d (mixed precision hides differences here), on `u u~ > u u~`: component madmatrix Fortran split-order driver rel QED=0 2.8276928588371750e+00 2.8276928588371737e+00 4.7e-16 QED=2 -5.5828746494657258e-02 -5.5828747824035893e-02 2.4e-08 QED=4 3.7810076327210495e-03 3.7810071013307570e-03 1.4e-07 The two EW rows do not agree to machine precision, and that is not new: the EW couplings differ between the Fortran and C++ builds at ~1e-7 relative. QED^2==4 prunes to a single amplitude order and never reaches any of this code, yet shows the same 1.4e-7 -- which is how one can tell. It is why the cross-backend tests in this tree compare at atol=1e-7 rather than in ulps. The check that does pin the interference algebra uses one backend, one set of momenta and one set of parameters: the three components sum to the unconstrained total to 1.6e-16 relative (2.7756451199752390 against 2.7756451199752394). And the pair loop with nothing masked reproduces the old single-jamp color sum bit for bit: `u u~ > t t~ QED^2<=4` gives 6.1741265892435804e-01 either way. Repeated on `u u~ > d d~ g`, where ncolor is 4 rather than 2 so the color triangle inside the pair loop is a real one: QED=0 2.0061172866596953e-04 (1.4e-15 from Fortran), QED=2 -1.5150831345723654e-05, QED=4 5.9366352531549790e-06. The two EW rows are 1.8e-11 and 1.2e-12 from Fortran in absolute terms -- 1.2e-6 and 2.0e-7 relative, the interference amplifying the coupling difference because it is a cancellation. The three components sum to the independently computed QED^2<=4 total BIT-IDENTICALLY. STILL NOT DONE, and unchanged by this: the other C++ exporters have no squared-order machinery either (`grep split_orders madgraph/iolibs/export_cpp.py` is 0 hits), so `output standalone_cpp` of `u u~ > t t~ QED^2==2` still builds and returns 6.1741266e-01. That wants either this treatment or a refusal, in its own change -- matchbox_cpp shares its base class and has IOTest goldens. Co-Authored-By: Claude Opus 5 --- madgraph/iolibs/export_v4.py | 83 +++++--- .../madmatrix/color_sum_splitorders.cc | 153 ++++++++++++++ .../madmatrix/process_class.inc | 2 +- .../process_function_definitions.inc | 6 +- .../madmatrix/process_matrix.inc | 10 +- .../madmatrix/process_sigmaKin_function.inc | 2 +- madmatrix/model_handling.py | 186 ++++++++++++++++-- madmatrix/output.py | 83 +++----- tests/acceptance_tests/test_cmd.py | 97 +++++---- 9 files changed, 476 insertions(+), 146 deletions(-) create mode 100644 madgraph/iolibs/template_files/madmatrix/color_sum_splitorders.cc diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index b72da22ce2..75f94213dd 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -266,40 +266,59 @@ def chosen_squared_orders(process, squared_orders): return res -def split_orders_not_supported_msg(format, process, split_orders, - squared_orders, chosen): - """Error text for a backend asked for a strict subset of squared orders. - - The backend sums every squared-order component into one |M|^2. That is the - number the user asked for as long as the constraint keeps all of them -- - which is why this is only an error when at least one is dropped. Naming the - components, and which of them survive, is the point: the difference between - what would be returned and what was asked for is otherwise invisible. +def split_order_tables(matrix_element): + """The squared split-order bookkeeping a backend needs, or None. + + Returns a dict with + + nampso how many amplitude split orders the amplitudes fall into + nsqampso how many squared orders their pairs produce + amp_so {amplitude number -> amplitude-order index}, 0-based + sqsoindex sqsoindex[m][n] -> squared-order index, 0-based. SYMMETRIC, + because a squared order is the SUM of the two amplitude + orders (Fortran SQSOINDEX), which is what lets a masked sum + stay real: (m,n) and (n,m) are kept or dropped together. + chosen [bool] per squared order, the user's constraint + names ['QED=0', ...] per squared order, for comments + + None when the process has no split orders, i.e. when there is one implicit + component and every backend already computes it. """ - def name(orders): - return ' '.join('%s=%d' % (order, value) - for order, value in zip(split_orders, orders)) or 'ALL_ORDERS' - - kept = [name(o) for o, keep in zip(squared_orders, chosen) if keep] - dropped = [name(o) for o, keep in zip(squared_orders, chosen) if not keep] - - return """The '%(format)s' output format does not support squared split orders. -Process: %(proc)s -Its matrix element has %(n)d squared-order components: %(all)s. -The constraint keeps %(kept)s and drops %(dropped)s, but this backend has no -squared-order mask: it would sum all %(n)d components and return that total, -which is not the requested contribution. -Use 'output standalone' for this process -- the Fortran standalone carries the -split-order machinery (SMATRIX_SPLITORDERS) and returns the masked total as -well as each component. Alternatively, constrain the process so that all its -squared-order components are kept.""" % { - 'format': format, - 'proc': process.nice_string().replace('Process: ', ''), - 'n': len(squared_orders), - 'all': ', '.join(name(o) for o in squared_orders), - 'kept': ', '.join(kept) if kept else 'nothing', - 'dropped': ', '.join(dropped)} + process = matrix_element.get('processes')[0] + split_orders = process.get('split_orders') + if not split_orders: + return None + squared_orders, amp_orders = matrix_element.get_split_orders_mapping() + if not squared_orders: + return None + + amp_so = {} + for iampso, (_orders, amp_numbers) in enumerate(amp_orders): + for namp in amp_numbers: + amp_so[namp] = iampso + + # The squared order a pair of amplitude orders lands in: add the two + # amplitude orders and look the sum up. A pair whose sum is not in the + # list cannot happen (the list is built from exactly these sums), but + # guard anyway rather than write a negative index into the generated code. + index_of = {tuple(sqso): i for i, sqso in enumerate(squared_orders)} + sqsoindex = [] + for m, (orders_m, _a) in enumerate(amp_orders): + row = [] + for n, (orders_n, _b) in enumerate(amp_orders): + key = tuple(om + on for om, on in zip(orders_m, orders_n)) + row.append(index_of.get(key, -1)) + sqsoindex.append(row) + + return {'nampso': len(amp_orders), + 'nsqampso': len(squared_orders), + 'amp_so': amp_so, + 'sqsoindex': sqsoindex, + 'chosen': chosen_squared_orders(process, squared_orders), + 'names': [' '.join('%s=%d' % (o, v) + for o, v in zip(split_orders, sqso)) + for sqso in squared_orders]} #=============================================================================== # ProcessExporterFortran diff --git a/madgraph/iolibs/template_files/madmatrix/color_sum_splitorders.cc b/madgraph/iolibs/template_files/madmatrix/color_sum_splitorders.cc new file mode 100644 index 0000000000..f9f2fcb7fc --- /dev/null +++ b/madgraph/iolibs/template_files/madmatrix/color_sum_splitorders.cc @@ -0,0 +1,153 @@ +// Copyright (C) 2020-2026 CERN and UCLouvain. +// Licensed under the GNU Lesser General Public License (version 3 or later). +// Created originally by: A. Valassi (Sep 2025) for the MG5aMC CUDACPP plugin. +// Further modified by: A. Valassi (2025). +// Integrated with the MadGraph7 project in Feb 2026. +// +// SQUARED SPLIT ORDERS. This is color_sum.cc for a process whose '^2' +// constraint leaves more than one amplitude split order, and it is a separate +// file because the sum itself changes shape rather than gaining a hole: the +// jamps arrive as nampso vectors end to end (njampso = nampso*ncolor) and the +// color sum pairs them, which is the Fortran GET_MATRIX contract +// JAMP(NCOLOR,NAMPSO) -> RES(NSQAMPSO). +// +// THE PAIR LOOP RUNS OVER ALL nampso*nampso ORDERED PAIRS, never a triangle. +// The color contraction below keeps the triangular matrix whose off-diagonal +// is doubled, which for a single jamp vector is just the statement that the +// (i,j) and (j,i) terms are conjugates. For two DIFFERENT jamp vectors they +// are not, so the value computed for one ordered pair (m,n) is NOT that pair's +// contribution -- only (m,n) and (n,m) together are, and they are equal to the +// true sum of the two. That is why both orderings must be summed, and why it +// is safe to mask on the squared order: sqSoIndex is symmetric (a squared +// order is the SUM of the two amplitude orders), so a pair and its transpose +// are always kept or dropped together and the masked total stays real. + +#include "color_sum.h" + +#include "mgOnGpuConfig.h" + +#include "MemoryAccessMatrixElements.h" + +#ifdef MGONGPUCPP_GPUIMPL +#error The squared split-order color sum is implemented for the CPU (SIMD) backend only. \ +The GPU jamp buffers are sized for a single jamp vector per helicity (ncolor, not njampso), \ +so a GPU build of this process would silently sum the wrong thing. Use a CPU backend, or \ +generate the process with a constraint that leaves a single amplitude split order. +#endif + +namespace mg5amcCpu +{ + constexpr int ncolor = CPPProcess::ncolor; // the number of leading colors + constexpr int nampso = CPPProcess::nampso; // the amplitude split orders + constexpr int njampso = CPPProcess::njampso; // ncolor * nampso: the jamps of every order, end to end + constexpr int nsqampso = CPPProcess::nsqampso; // the squared orders their pairs produce + + //-------------------------------------------------------------------------- + + // *** COLOR MATRIX BELOW *** +%(color_matrix_lines)s + + //-------------------------------------------------------------------------- + + // *** SQUARED SPLIT ORDERS BELOW *** +%(sqso_tables)s + + //-------------------------------------------------------------------------- + + void + color_sum_cpu( fptype* allMEs, // output: allMEs[nevt], add |M|^2 for one specific helicity + const cxtype_sv* allJamp_sv, // input: jamp_sv[njampso] (float/double) or jamp_sv[2*njampso] (mixed) for one specific helicity + const int ievt0 ) // input: first event number in current C++ event page (for CUDA, ievt depends on threadid) + { + // Pre-compute a constexpr triangular color matrix properly normalized #475 + struct TriangularNormalizedColorMatrix + { + // See https://stackoverflow.com/a/34465458 + __host__ __device__ constexpr TriangularNormalizedColorMatrix() + : value() + { + for( int icol = 0; icol < ncolor; icol++ ) + { + // Diagonal terms + value[icol][icol] = colorMatrix[icol][icol] / colorDenom[icol]; + // Off-diagonal terms + for( int jcol = icol + 1; jcol < ncolor; jcol++ ) + value[icol][jcol] = 2 * colorMatrix[icol][jcol] / colorDenom[icol]; + } + } + fptype2 value[ncolor][ncolor]; + }; + static constexpr auto cf2 = TriangularNormalizedColorMatrix(); + fptype_sv deltaMEs = { 0 }; +#if defined MGONGPU_CPPSIMD and defined MGONGPU_FPTYPE_DOUBLE and defined MGONGPU_FPTYPE2_FLOAT + fptype_sv deltaMEs_next = { 0 }; +#endif + // Gather the color flows the sum runs over, for each amplitude split order. + // The order index strides by ncolor, exactly as calculate_jamps wrote them. + // NB in mixed mode the two neppV vectors of allJamp_sv, at ijamp and at + // njampso+ijamp, are two halves of the event page and not two colors. + fptype2_sv jampR_sv[nampso][ncolor]; + fptype2_sv jampI_sv[nampso][ncolor]; + for( int iao = 0; iao < nampso; iao++ ) + { + for( int icol = 0; icol < ncolor; icol++ ) + { + const int ijamp = iao * ncolor + icol; +#if defined MGONGPU_CPPSIMD and defined MGONGPU_FPTYPE_DOUBLE and defined MGONGPU_FPTYPE2_FLOAT + jampR_sv[iao][icol] = fpvmerge( cxreal( allJamp_sv[ijamp] ), cxreal( allJamp_sv[njampso + ijamp] ) ); + jampI_sv[iao][icol] = fpvmerge( cximag( allJamp_sv[ijamp] ), cximag( allJamp_sv[njampso + ijamp] ) ); +#else + jampR_sv[iao][icol] = (fptype2_sv)( cxreal( allJamp_sv[ijamp] ) ); + jampI_sv[iao][icol] = (fptype2_sv)( cximag( allJamp_sv[ijamp] ) ); +#endif + } + } + // Loop over ORDERED pairs of amplitude split orders (see the note at the top: + // this must not be turned into a triangle over iao/jao) + for( int iao = 0; iao < nampso; iao++ ) + { + for( int jao = 0; jao < nampso; jao++ ) + { + // The squared order this pair contributes to, dropped if the process + // asked for a contribution that does not include it + if( !chosenSqso[sqSoIndex[iao][jao]] ) continue; + // Loop over icol + for( int icol = 0; icol < ncolor; icol++ ) + { + // Diagonal terms (the color contraction is still triangular) + fptype2_sv ztempR_sv = cf2.value[icol][icol] * jampR_sv[iao][icol]; + fptype2_sv ztempI_sv = cf2.value[icol][icol] * jampI_sv[iao][icol]; + // Loop over jcol + for( int jcol = icol + 1; jcol < ncolor; jcol++ ) + { + // Off-diagonal terms + ztempR_sv += cf2.value[icol][jcol] * jampR_sv[iao][jcol]; + ztempI_sv += cf2.value[icol][jcol] * jampI_sv[iao][jcol]; + } + // ztemp comes from the jamps of order iao, and is contracted with those + // of order jao (Fortran: ZTEMP from JAMP(:,M), times DCONJG(JAMP(I,N))) + fptype2_sv deltaMEs2 = ( jampR_sv[jao][icol] * ztempR_sv + jampI_sv[jao][icol] * ztempI_sv ); // may underflow #831 +#if defined MGONGPU_CPPSIMD and defined MGONGPU_FPTYPE_DOUBLE and defined MGONGPU_FPTYPE2_FLOAT + deltaMEs += fpvsplit0( deltaMEs2 ); + deltaMEs_next += fpvsplit1( deltaMEs2 ); +#else + deltaMEs += deltaMEs2; +#endif + } + } + } + // *** STORE THE RESULTS *** + using E_ACCESS = HostAccessMatrixElements; // non-trivial access: buffer includes all events + fptype* MEs = E_ACCESS::ieventAccessRecord( allMEs, ievt0 ); + // NB: color_sum ADDS |M|^2 for one helicity to the running sum of |M|^2 over helicities for the given event(s) + fptype_sv& MEs_sv = E_ACCESS::kernelAccess( MEs ); + MEs_sv += deltaMEs; // fix #435 +#if defined MGONGPU_CPPSIMD and defined MGONGPU_FPTYPE_DOUBLE and defined MGONGPU_FPTYPE2_FLOAT + fptype* MEs_next = E_ACCESS::ieventAccessRecord( allMEs, ievt0 + neppV ); + fptype_sv& MEs_sv_next = E_ACCESS::kernelAccess( MEs_next ); + MEs_sv_next += deltaMEs_next; +#endif + } + + //-------------------------------------------------------------------------- +} diff --git a/madgraph/iolibs/template_files/madmatrix/process_class.inc b/madgraph/iolibs/template_files/madmatrix/process_class.inc index 87c6242721..501a2c3793 100644 --- a/madgraph/iolibs/template_files/madmatrix/process_class.inc +++ b/madgraph/iolibs/template_files/madmatrix/process_class.inc @@ -58,7 +58,7 @@ static constexpr int ncolor = %(ncolor)s; // the number of leading colors: e.g. 1 for e+ e- -> mu+ mu- // The color structures the color flow is picked among. Same as ncolor, // unless the color sum runs on a smaller basis (see set_color_flow_lines_cpp). - static constexpr int ncolor_flow = %(ncolor_flow)s; + static constexpr int ncolor_flow = %(ncolor_flow)s;%(split_order_constants)s static constexpr int nmaxflavor = %(nmaxflavor)d; // the maximum number of flavor combinations // Hardcoded parameters for this process (constant class variables) diff --git a/madgraph/iolibs/template_files/madmatrix/process_function_definitions.inc b/madgraph/iolibs/template_files/madmatrix/process_function_definitions.inc index c115970337..c6db8e9ef1 100644 --- a/madgraph/iolibs/template_files/madmatrix/process_function_definitions.inc +++ b/madgraph/iolibs/template_files/madmatrix/process_function_definitions.inc @@ -84,7 +84,7 @@ namespace mg5amcCpu constexpr int ncomb = CPPProcess::ncomb; // #helicity combinations: e.g. 16 for e+ e- -> mu+ mu- (2**4 = fermion spin up/down ** npar) constexpr int ncolor = CPPProcess::ncolor; // the number of leading colors constexpr int ncolor_flow = CPPProcess::ncolor_flow; // the color structures the color flow is picked among - constexpr int nmaxflavor = CPPProcess::nmaxflavor; // the maximum number of flavor combinations + constexpr int nmaxflavor = CPPProcess::nmaxflavor; // the maximum number of flavor combinations%(jampso_aliases)s // [NB: I am currently unable to get the right value of nwf in CPPProcess.h - will hardcode it in CPPProcess.cc instead (#644)] //using CPPProcess::nwf; // #wavefunctions = #external (npar) + #internal: e.g. 5 for e+ e- -> mu+ mu- (1 internal is gamma or Z) @@ -576,9 +576,9 @@ namespace mg5amcCpu constexpr fptype_sv* jamp2_sv = nullptr; // no need for color selection during helicity filtering //std::cout << "sigmaKin_getGoodHel ihel=" << ihel << ( isGoodHel[ihel] ? " true" : " false" ) << std::endl; #if defined MGONGPU_CPPSIMD and defined MGONGPU_FPTYPE_DOUBLE and defined MGONGPU_FPTYPE2_FLOAT - cxtype_sv jamp_sv[2 * ncolor] = {}; // all zeros + cxtype_sv jamp_sv[2 * %(jamp_ncolor)s] = {}; // all zeros #else - cxtype_sv jamp_sv[ncolor] = {}; // all zeros + cxtype_sv jamp_sv[%(jamp_ncolor)s] = {}; // all zeros #endif calculate_jamps( ihel, allmomenta, allcouplings, hgFlavorVec, jamp_sv, false, allNumerators, allDenominators, jamp2_sv, ievt00 ); //maxtry? color_sum_cpu( allMEs, jamp_sv, ievt00 ); diff --git a/madgraph/iolibs/template_files/madmatrix/process_matrix.inc b/madgraph/iolibs/template_files/madmatrix/process_matrix.inc index f266f9b198..03aaa94d36 100644 --- a/madgraph/iolibs/template_files/madmatrix/process_matrix.inc +++ b/madgraph/iolibs/template_files/madmatrix/process_matrix.inc @@ -16,7 +16,7 @@ if( jamp2_sv ) // disable color choice if nullptr { for( int icol = 0; icol < ncolor_flow; icol++ ) - jamp2_sv[ncolor_flow * iParity + icol] += cxabs2( %(jamp_flow)s[icol] ); // may underflow #831 + jamp2_sv[ncolor_flow * iParity + icol] += cxabs2( %(jamp_flow_col)s ); // may underflow #831 } #else /* clang-format off */ assert( iParity == 0 ); // sanity check for J2_ACCESS @@ -25,7 +25,7 @@ { for( int icol = 0; icol < ncolor_flow; icol++ ) // NB: atomicAdd is needed after moving to cuda streams with one helicity per stream! - atomicAdd( &J2_ACCESS::kernelAccessIcol( colAllJamp2s, icol ), cxabs2( %(jamp_flow)s[icol] ) ); + atomicAdd( &J2_ACCESS::kernelAccessIcol( colAllJamp2s, icol ), cxabs2( %(jamp_flow_col)s ) ); } #endif /* clang-format on */ @@ -35,12 +35,12 @@ // In CUDA, copy the local jamp to the output global-memory jamp constexpr int ihel0 = 0; // the allJamps buffer already points to a specific helicity _within a super-buffer for dcNGoodHel helicities_ using J_ACCESS = DeviceAccessJamp; - for( int icol = 0; icol < ncolor; icol++ ) + for( int icol = 0; icol < %(jamp_ncolor)s; icol++ ) J_ACCESS::kernelAccessIcolIhelNhel( allJamps, icol, ihel0, dcNGoodHel ) = jamp_sv[icol]; #else // In C++, copy the local jamp to the output array passed as function argument - for( int icol = 0; icol < ncolor; icol++ ) - allJamp_sv[iParity * ncolor + icol] = jamp_sv[icol]; + for( int icol = 0; icol < %(jamp_ncolor)s; icol++ ) + allJamp_sv[iParity * %(jamp_ncolor)s + icol] = jamp_sv[icol]; #endif } // END LOOP ON IPARITY diff --git a/madgraph/iolibs/template_files/madmatrix/process_sigmaKin_function.inc b/madgraph/iolibs/template_files/madmatrix/process_sigmaKin_function.inc index ee74555d10..de51e2b556 100644 --- a/madgraph/iolibs/template_files/madmatrix/process_sigmaKin_function.inc +++ b/madgraph/iolibs/template_files/madmatrix/process_sigmaKin_function.inc @@ -122,7 +122,7 @@ %(cpp_blas_helicity_loop)s for( int ighel = 0; ighel < cNGoodHel; ighel++ ) { const int ihel = cGoodHel[ighel]; - cxtype_sv jamp_sv[nParity * ncolor] = {}; // fixed nasty bug (omitting 'nParity' caused memory corruptions after calling calculate_jamps) + cxtype_sv jamp_sv[nParity * %(jamp_ncolor)s] = {}; // fixed nasty bug (omitting 'nParity' caused memory corruptions after calling calculate_jamps) // **NB! in "mixed" precision, using SIMD, calculate_jamps computes MEs for TWO neppV pages with a single channelId! #924 bool storeChannelWeights = allChannelIds != nullptr || allrnddiagram != nullptr; calculate_jamps( ihel, allmomenta, allcouplings, iflavorVec, jamp_sv, storeChannelWeights, allNumerators, allDenominators, jamp2_sv, ievt00 ); diff --git a/madmatrix/model_handling.py b/madmatrix/model_handling.py index 5b5581fcb0..e41814347d 100644 --- a/madmatrix/model_handling.py +++ b/madmatrix/model_handling.py @@ -1677,6 +1677,17 @@ def get_process_class_definitions(self, write=True): replace_dict['nbhel'] = self.matrix_elements[0].get_helicity_combinations() # number of helicity combinations replace_dict['ndiagrams'] = len(self.matrix_elements[0].get('diagrams')) # AV FIXME #910: elsewhere matrix_element.get('diagrams') and max(config[0]... replace_dict['nmaxflavor'] = len(self.matrix_elements[0].get_external_flavors_with_iden()) # number of flavor combinations + # Only written when the jamps are actually split, so that a process + # without squared split orders keeps the header it always had + so = self.split_orders_info() + replace_dict['split_order_constants'] = '' if not self.split_orders_active() else ( + '\n // Squared split orders: the amplitudes fall into nampso amplitude' + '\n // orders, the jamps carry one vector per order (njampso long in total)' + '\n // and the color sum pairs them into nsqampso squared orders' + '\n // (see color_sum.cc, written from color_sum_splitorders.cc).' + '\n static constexpr int nampso = %d;' + '\n static constexpr int njampso = ncolor * nampso; // the jamps of every amplitude order, end to end' + '\n static constexpr int nsqampso = %d;' % (so['nampso'], so['nsqampso'])) replace_dict['nwave'] = 4 if (fd_gauge): replace_dict['nwave'] += 1 @@ -1699,6 +1710,12 @@ def get_process_function_definitions(self, write=True): """The complete class definition for the process""" replace_dict = super().get_process_function_definitions(write=False) # defines replace_dict['initProc_lines'] replace_dict['hardcoded_initProc_lines'] = replace_dict['initProc_lines'].replace( 'm_pars->', 'Parameters::') + replace_dict['jamp_ncolor'] = self.jamp_ncolor() + # Only pulled into scope when the jamps are split, so that a process + # without split orders keeps exactly the constants it always had + replace_dict['jampso_aliases'] = '' if not self.split_orders_active() else ( + '\n constexpr int nampso = CPPProcess::nampso; // the amplitude split orders' + '\n constexpr int njampso = CPPProcess::njampso; // ncolor * nampso: the jamps of every order, end to end') couplings2order_indep = [] ###replace_dict['ncouplings'] = len(self.couplings2order) ###replace_dict['ncouplingstimes2'] = 2 * replace_dict['ncouplings'] @@ -1924,6 +1941,7 @@ def get_sigmaKin_lines(self, color_amplitudes, write=True): replace_dict = super().get_sigmaKin_lines(color_amplitudes, write=False) replace_dict['proc_id'] = self.proc_id if self.proc_id>0 else 1 replace_dict['proc_id_source'] = 'MadMatrix exporter' + replace_dict['jamp_ncolor'] = self.jamp_ncolor() # Extract denominator (avoid to extend size for mirroring) den_factors = [str(me.get_denominator_factor()) for me in \ @@ -1961,6 +1979,9 @@ def get_sigmaKin_lines(self, color_amplitudes, write=True): def get_all_sigmaKin_lines(self, color_amplitudes, class_name): """Get sigmaKin_process for all subprocesses for CPPProcess.cc""" ret_lines = [] + # The jamps are one vector per amplitude split order, njampso long in + # total; 'ncolor' without them, so the default output is unchanged. + jamp_dim = 'njampso' if self.split_orders_active() else 'ncolor' if self.single_helicities: ###misc.sprint(type(self.helas_call_writer)) ###misc.sprint( 'before get_matrix_element_calls', self.matrix_elements[0].get_number_of_wavefunctions() ) # WRONG value of nwf, eg 7 for gg_tt @@ -2093,7 +2114,7 @@ def get_all_sigmaKin_lines(self, color_amplitudes, class_name): ret_lines.append(""" // Local variables for the given CUDA event (ievt) or C++ event page (ipagV) // [jamp: sum (for one event or event page) of the invariant amplitudes for all Feynman diagrams in a given color combination] - cxtype_sv jamp_sv[ncolor] = {}; // all zeros (NB: vector cxtype_v IS initialized to 0, but scalar cxtype is NOT, if "= {}" is missing!)""") + cxtype_sv jamp_sv[%s] = {}; // all zeros (NB: vector cxtype_v IS initialized to 0, but scalar cxtype is NOT, if "= {}" is missing!)""" % jamp_dim) # Shared sub-expressions of the color flows, filled in while the # amplitudes go by (see MadMatrixUFOHelasCallWriter.build_jamp_plan). # No "= {}": each one is assigned before it is ever read. @@ -2232,6 +2253,51 @@ def blas_available_flags(cls): return '' return cls._blas_flags + # AV - new method (add the split-order holes to process_matrix.inc) + def get_matrix_single_process(self, i, matrix_element, color_amplitudes, + class_name, write=True): + replace_dict = super().get_matrix_single_process( + i, matrix_element, color_amplitudes, class_name, write=False) + replace_dict['jamp_ncolor'] = self.jamp_ncolor() + # set_color_flow_lines_cpp fills jamp_flow / jamp_flow_col; it runs from + # get_process_class_definitions, before this, but be explicit rather + # than rely on the ordering of two independent methods. + if 'jamp_flow_col' not in replace_dict: + self.set_color_flow_lines_cpp(matrix_element, replace_dict) + if write: + return self.read_template_file(self.single_process_template) % replace_dict + return replace_dict + + # AV - new method + def jamp_ncolor(self): + """The length of a jamp array: 'ncolor', or 'njampso' (= ncolor*nampso) + once the jamps carry an amplitude-order index. Templates spell the size + through this hole so that a process without split orders gets exactly + the text it got before they existed.""" + return 'njampso' if self.split_orders_active() else 'ncolor' + + # AV - new method + def split_orders_info(self): + """The squared split-order tables for this process, or None. + + None means the process has no '^2' constraint, and then every hole this + fills reproduces the code that was written before split orders existed + -- which is how the default path stays byte-for-byte what it was. + """ + if not hasattr(self, '_split_orders_info'): + from madgraph.iolibs.export_v4 import split_order_tables + self._split_orders_info = split_order_tables(self.matrix_elements[0]) + return self._split_orders_info + + def split_orders_active(self): + """Whether the jamps carry an amplitude-order index at all. + + A single amplitude order is the same code as none: the pair loop of the + color sum has one term, so it is left switched off rather than writing + a 1x1 interference matrix.""" + so = self.split_orders_info() + return bool(so) and so['nampso'] > 1 + # AV - new method @classmethod def cpp_blas_wanted_for(cls, ncolor): @@ -2245,17 +2311,53 @@ def cpp_blas_wanted_for(cls, ncolor): return ncolor >= cls.blas_min_ncolor def cpp_blas_wanted(self): + # The BLAS color sum multiplies a single jamp vector per helicity by the + # color matrix. With the jamps split by amplitude order the sum is a + # pair loop instead, so the scalar color sum is used for those. + if self.split_orders_active(): + return False return self.cpp_blas_wanted_for( max(1, len(self.matrix_elements[0].get('color_basis')))) + # AV - new method + def get_sqso_table_lines(self): + """The interference matrix and the squared-order mask, as C++ tables. + + sqSoIndex[m][n] is the Fortran SQSOINDEX: the squared order a pair of + amplitude orders lands in. It is symmetric because a squared order is + the SUM of the two amplitude orders, which is what lets color_sum_cpu + mask on it and still get a real answer.""" + so = self.split_orders_info() + lines = [] + lines.append(' // The squared order each pair of amplitude orders' + ' contributes to (symmetric)') + lines.append(' static constexpr int sqSoIndex[nampso][nampso] = {') + lines.append(',\n'.join(' { %s }' % ', '.join(str(i) for i in row) + for row in so['sqsoindex'])) + lines.append(' };') + lines.append(' // The squared orders the process asked for:') + for k, (name, keep) in enumerate(zip(so['names'], so['chosen'])): + lines.append(' // %d) %s%s' % (k, name, + '' if keep else ' [dropped]')) + lines.append(' static constexpr bool chosenSqso[nsqampso] = { %s };' + % ', '.join('true' if k else 'false' + for k in so['chosen'])) + return '\n'.join(lines) + # AV - new method def edit_colorsum(self): """Generate color_sum.cc""" ###misc.sprint('Entering OneProcessExporterMadMatrix.edit_colorsum') - template = open(pjoin(self.template_path,'madmatrix','color_sum.cc'),'r').read() + # A process whose '^2' constraint leaves more than one amplitude split + # order gets the dedicated pair-loop color sum instead (see that file). + split = self.split_orders_active() + name = 'color_sum_splitorders.cc' if split else 'color_sum.cc' + template = open(pjoin(self.template_path,'madmatrix',name),'r').read() replace_dict = {} # Extract color matrix again (this was also in get_matrix_single_process called within get_all_sigmaKin_lines) replace_dict['color_matrix_lines'] = self.get_color_matrix_lines(self.matrix_elements[0]) + if split: + replace_dict['sqso_tables'] = self.get_sqso_table_lines() replace_dict['cpp_blas_color_sum'] = '' if self.cpp_blas_wanted(): replace_dict['cpp_blas_color_sum'] = strip_banner( @@ -2412,23 +2514,50 @@ def set_color_flow_lines_cpp(self, matrix_element, replace_dict): color_basis = matrix_element.get('color_basis') flow_basis = color_basis.get_flow_basis() if color_basis else None + # A color flow is picked from the WHOLE amplitude, so when the jamps are + # split by amplitude order the flow is built from their sum: the color + # flow amplitude of flow i is sum_M jamp_M[i], and squaring that gives + # back exactly today's jamp2 when there is a single amplitude order. + # (The squared-order mask is deliberately not applied here: this weight + # only chooses which color flow to write out, it is not an observable.) + so = self.split_orders_info() + split = bool(so) and so['nampso'] > 1 + jamp_src = 'jamp_sv' + so_lines = [] + if split: + jamp_src = 'jampso_sv' + so_lines = ['', + ' // The color flows of the whole amplitude: the' + ' jamps of every', + ' // amplitude split order summed back together' + ' (see color_sum.cc,', + ' // which instead keeps them apart and pairs' + ' them)', + ' cxtype_sv jampso_sv[ncolor] = {};', + ' for( int icol = 0; icol < ncolor; icol++ )', + ' for( int iao = 0; iao < nampso; iao++ )', + ' jampso_sv[icol] += jamp_sv[iao * ncolor +' + ' icol];'] + if flow_basis is None or flow_basis is color_basis: replace_dict['ncolor_flow'] = replace_dict['ncolor'] - replace_dict['jampflow_lines'] = '' - replace_dict['jamp_flow'] = 'jamp_sv' + replace_dict['jampflow_lines'] = '\n'.join(so_lines) + replace_dict['jamp_flow'] = jamp_src + replace_dict['jamp_flow_col'] = '%s[icol]' % jamp_src return projection = color_basis.get_flow_projection() - lines = ['', + lines = so_lines + ['', ' // The color flow jamps, rebuilt from the ones entering', ' // the color sum through the Kleiss-Kuijf relations', ' cxtype_sv jampf_sv[ncolor_flow] = {};'] for i, coeff_list in enumerate(projection): - terms = ''.join('%sjamp_sv[%d]' % (self.coeff(coefficient[0], - coefficient[1], - coefficient[2], - coefficient[3]), - number - 1) + terms = ''.join('%s%s[%d]' % (self.coeff(coefficient[0], + coefficient[1], + coefficient[2], + coefficient[3]), + jamp_src, + number - 1) for coefficient, number in coeff_list) lines.append(' jampf_sv[%d] = %s;' % (i, terms if terms else 'cxzero_sv()')) @@ -2436,6 +2565,7 @@ def set_color_flow_lines_cpp(self, matrix_element, replace_dict): replace_dict['ncolor_flow'] = max(1, len(flow_basis)) replace_dict['jampflow_lines'] = '\n'.join(lines) replace_dict['jamp_flow'] = 'jampf_sv' + replace_dict['jamp_flow_col'] = 'jampf_sv[icol]' logger.debug('Color sum on %d DDM structures, color flow on %d trace ' 'structures (%d Kleiss-Kuijf terms)', @@ -2814,6 +2944,20 @@ def jamp_statement(cls, target, terms, assign): pieces.append('%s %s%s' % ('-' if sign < 0 else '+', factor, name)) return '%s %s %s;' % (target, '=' if assign else '+=', ' '.join(pieces)) + def split_order_index(self, matrix_element): + """{amplitude number -> amplitude-order index}, or None. + + None whenever the jamps are a single vector, i.e. no '^2' constraint or + only one amplitude order, in which case every jamp index below is the + plain color index it always was.""" + key = id(matrix_element) + if getattr(self, '_so_key', None) != key: + from madgraph.iolibs.export_v4 import split_order_tables + so = split_order_tables(matrix_element) + self._so_key = key + self._so_index = so['amp_so'] if so and so['nampso'] > 1 else None + return self._so_index + def build_jamp_plan(self, matrix_element, color_amplitudes): """Work out how the color flows are built from shared sub-expressions, and return (ntmp, captures, combines, final): @@ -2828,6 +2972,13 @@ def build_jamp_plan(self, matrix_element, color_amplitudes): if not self.jamp_optim_enabled(): return None + # A shared sub-expression is a partial sum of amplitudes. With the jamps + # split by amplitude order, two amplitudes of different orders must not + # end up in the same partial sum -- it would be added to one jamp vector + # and the order information lost. Write the flows out one amplitude at a + # time instead. + if self.split_order_index(matrix_element) is not None: + return None all_element = self.jamp_matrix(color_amplitudes) if not all_element: return None @@ -2921,6 +3072,11 @@ def super_get_matrix_element_calls(self, matrix_element, color_amplitudes, multi # one (color flow, amplitude) pair at a time, as before) jamp_plan = self.build_jamp_plan(matrix_element, color_amplitudes) self.nb_tmp_jamp = jamp_plan[0] if jamp_plan else 0 + so_index = self.split_order_index(matrix_element) + ncolor_jamp = len(color_amplitudes) + # 'ncolor' unless the jamps carry an amplitude-order index, so that a + # process without split orders gets exactly the text it always got + jamp_dim = 'njampso' if so_index is not None else 'ncolor' if jamp_plan is not None: _ntmp, jamp_captures, jamp_combines, jamp_final = jamp_plan me = matrix_element.get('diagrams') @@ -2991,7 +3147,7 @@ def super_get_matrix_element_calls(self, matrix_element, color_amplitudes, multi 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 < ncolor; i++ ) { jamp_sv[i] = cxzero_sv(); } + for( int i = 0; i < """ + jamp_dim + """; i++ ) { jamp_sv[i] = cxzero_sv(); } // Numerators for the current event (CUDA) or SIMD event page (C++) // (denominators are no longer accumulated here: they are derived as the sum of numerators later) @@ -3110,6 +3266,13 @@ def _guard_open(group_mask): gmask = diag_group_mask.get(id(diagram)) before_guard = [] if jamp_plan is None: + # With split orders the jamps are nampso vectors end to end + # and this amplitude belongs to exactly one of them, so its + # color flows are offset onto that vector. so_index is None + # otherwise and the index is the plain color index. + jamp_offset = 0 + if so_index is not None: + jamp_offset = so_index.get(namp, 0) * ncolor_jamp for njamp, coeff in color[namp].items(): scoeff = OneProcessExporterMadMatrix.coeff(*coeff) # AV if scoeff[0] == '+' : scoeff = scoeff[1:] @@ -3118,6 +3281,7 @@ def _guard_open(group_mask): scoeff = scoeff.replace(',',', ') scoeff = scoeff.replace('*',' * ') scoeff = scoeff.replace('/',' / ') + njamp = njamp + jamp_offset if scoeff.startswith('-'): amp_block.append('jamp_sv[%s] -= %samp_sv[0];' % (njamp, scoeff[1:])) # AV else: amp_block.append('jamp_sv[%s] += %samp_sv[0];' % (njamp, scoeff)) # AV else: diff --git a/madmatrix/output.py b/madmatrix/output.py index b8cfa48784..d367b8c915 100644 --- a/madmatrix/output.py +++ b/madmatrix/output.py @@ -26,7 +26,6 @@ import madgraph.iolibs.export_v4 as export_v4 import madgraph.iolibs.export_cpp as export_cpp import madgraph.various.misc as misc -from madgraph import InvalidCmd from . import launch_plugin @@ -212,59 +211,41 @@ def write_p_makefiles(self, model): open(pjoin(self.dir_path, 'SubProcesses', name), 'w').write(rendered) def check_split_orders(self, matrix_element): - """Refuse a process whose squared-order constraint drops a component. - - This backend has no squared-order machinery at all: it builds one jamp - vector per helicity and contracts it with the color matrix once, so the - |M|^2 it returns is the sum over *every* squared-order component the - amplitude has. That is exactly the requested number whenever the - constraint keeps all of them -- which covers the single-component case - (MG5 has already resolved the constraint at generation, e.g. - `p p > j j QCD^2==4`) and the multi-component case that keeps - everything (e.g. `u u~ > t t~ QED^2<=4`). Both stay supported, and - untouched. - - It is not the requested number when the constraint drops a component, - which happens when that component cannot be reached by dropping - diagrams -- an interference term, typically: `u u~ > t t~ QED^2==2` - keeps all three diagrams and asks for the QCD-EW cross term alone, - while this backend would return the full QCD^2 + interference + EW^2 - total. There is no mask to apply here, so it says so instead. - - Supporting it would mean splitting the jamps by amplitude order and - pairing them in the color sum, the Fortran GET_MATRIX contract - (JAMP(NCOLOR,NAMPSO) -> RES(NSQAMPSO)). All three color sums here -- - color_sum_cpu, the CUDA kernel and the BLAS path -- take a single jamp - vector per helicity, as do the color folding, the C-parity de-duplication - and the good-helicity scan that read their output. Masking alone would at - least keep the scalar signature; it is the jamp splitting underneath that - this backend has no room for, which is why this refuses rather than - computing the requested contribution. + """Report what a squared-order constraint will produce here. + + Supported: the jamps carry an amplitude-order index and the color sum + pairs them (color_sum_splitorders.cc, the Fortran GET_MATRIX contract), + so a '^2' constraint that keeps only some squared orders gets the + contribution it asked for rather than the total. That is what makes the + interference case work -- `u u~ > t t~ QED^2==2` keeps all three + diagrams and wants the QCD-EW cross term alone, which no amount of + dropping diagrams at generation can produce. + + Not supported: a GPU build of such a process. The device jamp buffers + are sized for one jamp vector per helicity (ncolor, not njampso), and + the backend is a make-time choice rather than an output-time one, so + the refusal cannot live here: color_sum_splitorders.cc #errors under + MGONGPUCPP_GPUIMPL instead. Say so now rather than let a GPU build be + the first the user hears of it. """ - process = matrix_element.get('processes')[0] - split_orders = process.get('split_orders') - if not split_orders: - return - squared_orders, _amp_orders = matrix_element.get_split_orders_mapping() - chosen = export_v4.chosen_squared_orders(process, squared_orders) - if all(chosen): - # The total is what was asked for. The components it is made of are - # not available here though, and the Fortran standalone would hand - # them back through SMATRIX_SPLITORDERS, so say which of the two - # this output is rather than let the user assume the other. - if len(squared_orders) > 1: - logger.warning( - "%s returns the summed |M|^2 over all %d squared-order " - "components of '%s', which is what its constraint asks " - "for. The individual components are not available from " - "this backend -- use 'output standalone' for those.", - self.__class__.format_name, len(squared_orders), - process.nice_string().replace('Process: ', '')) + so = export_v4.split_order_tables(matrix_element) + if not so or so['nampso'] <= 1: return - raise InvalidCmd(export_v4.split_orders_not_supported_msg( - self.__class__.format_name, process, split_orders, - squared_orders, chosen)) + process = matrix_element.get('processes')[0] + kept = [n for n, k in zip(so['names'], so['chosen']) if k] + dropped = [n for n, k in zip(so['names'], so['chosen']) if not k] + logger.info( + "%s: '%s' has %d squared-order components (%s); keeping %s%s. " + "The jamps are split over %d amplitude orders and the color sum " + "pairs them; CPU backends only (a GPU build of this process will " + "not compile, by design).", + self.__class__.format_name, + process.nice_string().replace('Process: ', ''), + so['nsqampso'], ', '.join(so['names']), + ', '.join(kept) if kept else 'nothing', + '' if not dropped else ', dropping %s' % ', '.join(dropped), + so['nampso']) # AV - add debug printouts (in addition to the default one from OM's tutorial) def generate_subprocess_directory(self, matrix_element, cpp_helas_call_writer, proc_number=None): diff --git a/tests/acceptance_tests/test_cmd.py b/tests/acceptance_tests/test_cmd.py index 380ad0e68c..49caf89e07 100755 --- a/tests/acceptance_tests/test_cmd.py +++ b/tests/acceptance_tests/test_cmd.py @@ -1334,54 +1334,67 @@ def get_values(output_format, check_exe, build_source=False): 'all matrix elements vanished for u u~ > j j') self._assert_me_lists_close(mg7, standalone, atol=1e-7) - def test_standalone_mg7_refuses_a_dropped_squared_order(self): - """standalone_mg7 must refuse a squared-order constraint it cannot apply. - - The madmatrix backend contracts one jamp vector per helicity with the - color matrix once, so the |M|^2 it returns is the sum over every - squared-order component of the amplitude. Whenever the user's '^2' - constraint keeps all of them that is the requested number, and those - cases stay supported; when it drops one, the backend has no mask and - would hand back the full total instead. - - ``u u~ > t t~ QED^2==2`` is the case that matters: the interference - cannot be reached by dropping diagrams, so all three survive and the - Fortran masks the sum down to the QCD-EW cross term alone (-3.4e-18 at - the standard RAMBO point) while madmatrix returned the whole - 0.617412658924358 -- silently, with generated code differing from the - unconstrained process by one comment line. It has to say so instead. - - The supported cases are pinned alongside it, since the fix must not - turn them into errors: ``QED^2<=4`` keeps all three components and - ``QED^2==4`` is resolved by MG5 at generation down to a single one. + def test_standalone_mg7_split_orders_interference(self): + """standalone_mg7 must return the squared-order contribution asked for. + + The madmatrix jamps carry an amplitude-order index and the color sum + pairs them, so a '^2' constraint that keeps only some of the squared + orders gets that contribution and not the total. The case that matters + is an interference term, which cannot be reached by dropping diagrams + at generation: ``u u~ > u u~ QED^2==2`` keeps every diagram and wants + the QCD-EW cross term alone, -5.5828747824035893e-02 from the Fortran + split-order driver, where a backend with no mask returns the whole + +2.7756451181144683. + + The three components are checked to sum back to the unconstrained + total *as computed by this same backend*. That comparison is the one + that pins the pair loop: it uses one set of momenta and one set of + parameters, so it is sensitive to the interference algebra alone -- + in particular to the fact that the color contraction keeps its + doubled triangle while the loop over amplitude-order pairs must run + over all ordered pairs, since for two different jamp vectors a pair + and its transpose are not each other's conjugate. """ - self.do('import model sm') + devnull = open(os.devnull, 'w') + me_re = re.compile(r'Matrix element\s*=\s*([\d.eE+-]+)\s*GeV', + re.IGNORECASE) - # Dropped component: refuse, and name what was dropped. - self.do('generate u u~ > t t~ QED^2==2') - try: - self.do('output standalone_mg7 %s -f' % self.out_dir) - except InvalidCmd as error: - message = str(error) - else: - self.fail('standalone_mg7 accepted a squared-order constraint it ' - 'cannot apply') - self.assertIn('does not support squared split orders', message) - self.assertIn('QED=0', message) # the components it would have summed - self.assertIn('QED=4', message) - - # All components kept: the plain sum is the requested number. - for process in ('u u~ > t t~ QED^2<=4', # three components, all kept - 'u u~ > t t~ QED^2==4'): # resolved at generation + def value(constraint): if os.path.isdir(self.out_dir): shutil.rmtree(self.out_dir) - self.do('generate %s' % process) + self.do('generate u u~ > u u~ %s' % constraint) self.do('output standalone_mg7 %s -f' % self.out_dir) proc_root = pjoin(self.out_dir, 'SubProcesses') - self.assertTrue([d for d in os.listdir(proc_root) - if d.startswith('P') and - os.path.isdir(pjoin(proc_root, d))], - 'standalone_mg7 wrote no subprocess for %s' % process) + dirs = [d for d in os.listdir(proc_root) + if d.startswith('P') and os.path.isdir(pjoin(proc_root, d))] + self.assertTrue(dirs, 'no subprocess for %s' % constraint) + proc_dir = pjoin(proc_root, dirs[0]) + # FPTYPE=d: mixed precision hides and fakes differences here + self.assertEqual(0, subprocess.call(['make', 'FPTYPE=d'], + stdout=devnull, stderr=devnull, + cwd=proc_dir), + 'standalone_mg7 %s did not build' % constraint) + log = pjoin(proc_dir, 'check.log') + subprocess.call('./check_sa.exe 1000', shell=True, cwd=proc_dir, + stdout=open(log, 'w'), stderr=subprocess.STDOUT) + found = me_re.findall(open(log).read()) + self.assertTrue(found, 'no matrix element (see %s)' % log) + return float(found[0]) + + self.do('import model sm') + interference = value('QED^2==2') + # The Fortran split-order component, to the tolerance this backend is + # compared at elsewhere (the EW couplings differ in the last digits) + self.assertAlmostEqual(interference, -5.5828747824035893e-02, delta=1e-7) + # ... and emphatically not the unmasked total + self.assertLess(abs(interference), 1.0) + + components = [value('QED^2==0'), interference, value('QED^2==4')] + total = value('QED^2<=4') + self.assertAlmostEqual(sum(components), total, + delta=1e-12 * abs(total), + msg='the squared-order components do not add up ' + 'to the total this backend computes') def test_standalone_cpp(self): """test that the scalar C++ standalone exporter is working From a359fdec4f279e0b50a722aae8322860d6a6cd95 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Thu, 10 Sep 2026 14:18:28 +0200 Subject: [PATCH 3/3] madmatrix split orders: adapt the port to this branch The two commits above were written on a base that predates the big output format renaming and the removal of the color reflection folding, so replaying them here needs three adjustments that are not part of the feature: 'standalone_mg7' is now spelled 'standalone' (the Fortran one is 'standalone_fortran'). The exporter named itself by the old name, which only reached the user through the informational message, and the test drove the old name on the 'output' line -- where it is not an error but a directory name, so the test silently exercised the default format instead. color_sum_splitorders.cc summed over ncolorfold with a colorFoldRep[] map. There is no folding on this branch, so it sums over ncolor directly, which is the same code the folding degenerated to when it found nothing to fold. The expected interference is re-measured here: -5.5828746494657265e-02 from standalone_fortran on this base, which the C++ reproduces to 1.3e-16 relative rather than the ~1e-7 the couplings used to differ by. Verified on u u~ > u u~ at the standard RAMBO point, FPTYPE=d: QED=0 2.8276928588371750e+00 QED=2 -5.5828746494657258e-02 (Fortran: -5.5828746494657265e-02) QED=4 3.7810076327210495e-03 sum 2.7756451199752387e+00 total 2.7756451199752394e+00 (QED^2<=4, computed independently) and the no-split-orders path is byte-identical to the base commit for u u~ > u u~, g g > t t~ and p p > j j QCD^2==4. Co-Authored-By: Claude Opus 5 --- madmatrix/output.py | 2 +- tests/acceptance_tests/test_cmd.py | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/madmatrix/output.py b/madmatrix/output.py index d367b8c915..d9773eef08 100644 --- a/madmatrix/output.py +++ b/madmatrix/output.py @@ -282,7 +282,7 @@ def modify_grouping(self, matrix_element): # so that when running `make` in a P* folder, it builds check_sa.exe as well as the process library (predicatable behaviour) class ProcessExporterMadMatrixStandalone(ProcessExporterMadMatrix): - format_name = 'standalone_mg7' + format_name = 'standalone' # Each P* directory links madmatrix_standalone.mk (which itself includes # madmatrix.mk) as its 'makefile'; both have to be rendered in SubProcesses/ diff --git a/tests/acceptance_tests/test_cmd.py b/tests/acceptance_tests/test_cmd.py index 49caf89e07..bbe70b49d6 100755 --- a/tests/acceptance_tests/test_cmd.py +++ b/tests/acceptance_tests/test_cmd.py @@ -1334,17 +1334,17 @@ def get_values(output_format, check_exe, build_source=False): 'all matrix elements vanished for u u~ > j j') self._assert_me_lists_close(mg7, standalone, atol=1e-7) - def test_standalone_mg7_split_orders_interference(self): - """standalone_mg7 must return the squared-order contribution asked for. + def test_standalone_split_orders_interference(self): + """standalone (madmatrix) must return the squared-order contribution asked for. The madmatrix jamps carry an amplitude-order index and the color sum pairs them, so a '^2' constraint that keeps only some of the squared orders gets that contribution and not the total. The case that matters is an interference term, which cannot be reached by dropping diagrams at generation: ``u u~ > u u~ QED^2==2`` keeps every diagram and wants - the QCD-EW cross term alone, -5.5828747824035893e-02 from the Fortran + the QCD-EW cross term alone, -5.5828746494657265e-02 from the Fortran split-order driver, where a backend with no mask returns the whole - +2.7756451181144683. + +2.7756451199752394. The three components are checked to sum back to the unconstrained total *as computed by this same backend*. That comparison is the one @@ -1363,7 +1363,7 @@ def value(constraint): if os.path.isdir(self.out_dir): shutil.rmtree(self.out_dir) self.do('generate u u~ > u u~ %s' % constraint) - self.do('output standalone_mg7 %s -f' % self.out_dir) + self.do('output standalone %s -f' % self.out_dir) proc_root = pjoin(self.out_dir, 'SubProcesses') dirs = [d for d in os.listdir(proc_root) if d.startswith('P') and os.path.isdir(pjoin(proc_root, d))] @@ -1373,7 +1373,7 @@ def value(constraint): self.assertEqual(0, subprocess.call(['make', 'FPTYPE=d'], stdout=devnull, stderr=devnull, cwd=proc_dir), - 'standalone_mg7 %s did not build' % constraint) + 'standalone %s did not build' % constraint) log = pjoin(proc_dir, 'check.log') subprocess.call('./check_sa.exe 1000', shell=True, cwd=proc_dir, stdout=open(log, 'w'), stderr=subprocess.STDOUT) @@ -1385,7 +1385,7 @@ def value(constraint): interference = value('QED^2==2') # The Fortran split-order component, to the tolerance this backend is # compared at elsewhere (the EW couplings differ in the last digits) - self.assertAlmostEqual(interference, -5.5828747824035893e-02, delta=1e-7) + self.assertAlmostEqual(interference, -5.5828746494657265e-02, delta=1e-7) # ... and emphatically not the unmasked total self.assertLess(abs(interference), 1.0)