Cheaper colour sum: DDM basis, Kleiss-Kuijf flows, reflection folding, BLAS, and the same for madmatrix - #60
Open
oliviermattelaer wants to merge 50 commits into
Open
Conversation
A color matrix entry is the full contraction of two color structures, so it only depends on the relative labelling of the indices. Any permutation of the external indices mapping the color basis onto itself therefore satisfies C[P(i)][P(j)] = C[i][j], and only one row per orbit has to be computed. ColorBasisSymmetry finds those permutations; for g g > n g the action is transitive, so the whole matrix follows from a single row. The matrix is no longer stored entry by entry either: the distinct values are kept once and the (i,j) grid only holds a compact index into them, so the N^2 dictionaries of ColorFactor objects are gone. The never-read inverted_col_matrix dictionary becomes a property built on demand. In ColorBasis, a recycled color factor only has its indices relabelled, and a relabelled simplified expression is still simplified, so it only needs putting back in canonical form rather than running the full simplification over every term. The equivalence is checked once per canonical representation and the old path is kept where it does not hold. ColorFactor.simplify also looks up similar strings through a dictionary instead of scanning what it has accumulated, which was quadratic in the number of terms. Measured on g g > 6g (5040 color structures, 34300 diagrams): color matrix 194s/2126MB -> 7s/92MB, color basis 167s -> 38s, get_color_amplitudes 4599MB -> 1023MB, peak 8.9GB -> 3.4GB. Generated output is unchanged: 1002 source files identical across standalone, madevent, matchbox and standalone_cpp, and matrix.f for g g > 5g is byte for byte the same. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
optimise_jamp looks for pairs of columns whose ratio is the same on many lines. It found the candidates by looking every column of the 0..nb_col range up in the matrix, for every non zero entry. That matrix is very sparse - a color flow only gets a small share of the amplitudes, 322 out of 7245 for g g > 5g - so nearly all of that work was spent discovering zeros, and the cost was the number of non zero entries times the number of columns. Index the non zero entries by line instead, and walk only those. The pairs are visited in the same order as before, so the sub-expressions found and the order they are defined in are unchanged. The substitution step gets the same treatment through an index by column, and the count of lines sharing a ratio now lives in a single dictionary keyed by the two columns and the ratio at once rather than in nested dictionaries. g g > 5g: get_JAMP_lines 33.5s -> 5.2s, and the whole standalone output 29.3s -> 14.3s. g g > 6g (5040 color flows, 126630 amplitudes, 8.1M non zero entries) now goes through in 270s where the scan was out of reach. Generated output is unchanged: the same 1002 source files across standalone, madevent, matchbox and standalone_cpp, and matrix.f for g g > 5g is byte for byte the same as before. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every amplitude coefficient carries a power of the number of colors, and Fraction(3)**power was built again for each of them: 231840 times for g g > 5g, 8.1M times for g g > 6g, for a handful of distinct powers. This is a small win only. The Fraction time left in get_JAMP_lines is not in building those powers but in the multiplications mixing complex, int and Fraction, which fall back on the reflected operators. Reordering the product to keep it rational until the end would avoid that, but it would also move where the rounding to floating point happens, and optimise_jamp compares the resulting ratios for equality, so the generated code would change. g g > 5g get_JAMP_lines: 3.50s -> 3.46s, within the run to run spread. Generated output unchanged, the same 1002 source files as before. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ion sum Two independent changes to the standalone output, both in the color matrix. Every line of the color matrix is the line of its orbit representative with the columns permuted, in every process looked at: pure gluon, one quark line, several quark lines, quark initiated. So instead of writing the N*(N+1)/2 entries out, write one line per orbit, the permutations of a spanning subset of the generators, and for each line which generator reaches it from which other line. INIT_CF follows those back to the representative on the first call and rebuilds the matrix. This is only taken when it is at least four times smaller, which leaves everything below about a hundred color structures on the path it was on. g g > 5g: 259560 -> 5760 numbers, matrix.f 3.5MB -> 2.3MB. g g > 6g: 12703320 -> 45360 numbers, the color block of matrix.f 73MB -> 0.2MB and the file itself 114MB -> 41MB. The contraction of the matrix with the JAMPs now sums over four accumulators rather than one. A single accumulator makes every term wait for the one before it to leave the adder, and that latency, not the arithmetic, is what the loop was spending its time on: 1.55x faster on the generated g g > 5g code, 1.6x to 2x in a standalone benchmark up to 5040 colors. Narrowing the coefficients to one byte or widening them to a double changes nothing, and no compiler flag does this by itself, -ffast-math included. Reordering the sum moves the rounding, so |M|^2 changes in the last bit: checked against the previous output on g g > 4g, g g > tt~ gg, uu~ > uu~ gg, gg > uu~ dd~ g, g g > 3g and g g > 5g, all agreeing to one ulp or exactly. The color basis and color matrix themselves are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A permutation of the external color indices which maps the color basis onto itself also permutes the columns of the JAMP matrix up to a sign, so the whole matrix is invariant and the sub-expressions the optimisation looks for come in orbits: every one of them is reused exactly as often as the others. optimise_jamp took them one at a time, in whatever order the argmax set happened to be in, so its result was not closed under the symmetry: only about a fifth of the definitions had their image among the definitions. Taking a whole orbit at a time instead keeps the matrix invariant at every step. Two sub-expressions of the same orbit never want the same entry of the matrix, and the contention between orbits is settled by applying one orbit at a time. That also compresses better, not worse, everywhere it was measured, since the orbit which is applied is no longer eaten into by the ones applied before it: g g > 3g 72 -> 66 operations, 4g 951 -> 795, 5g 22221 -> 15750, 6g 441943 -> 237510, t t~ g g 281 -> 264, u u~ > u u~ g g 153 -> 148, g g > u u~ d d~ g 1753 -> 1680. Only one line per orbit is then written out, plus the amplitude permutations, and INIT_JAMP walks each orbit once on the first call to work out the operands of every definition. For g g > 6g that is 157 recipes and 6 permutations instead of 441943 lines: matrix.f goes from 43.4 MB to 19.4 MB and the GET_JAMP block from 33.6 MB to 3.69 MB. The operands are read from one array holding the amplitudes first and the definitions after them, which is why AMP is declared longer. Below five thousand definitions the lines are still both smaller and faster written out, which is where the threshold comes from. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The walk which rebuilds the definitions only knew how to carry a sign, so any process whose color coefficients bring an i fell back to writing every definition out. That is every process with a heavy quark line: g g > t t~ g g, u u~ > u u~ g g, and g g > t t~ g g g g among them. All the factors which turn up are powers of i, so they form a group of four elements and can be carried as an exponent modulo four. The walk stays integer arithmetic: the image of A + i**e*B is A' + i**(e+eb-ea)*B', and the column it defines is i**ea times the image of the column A + i**e*B defines. Only the array holding the factor in front of the second operand has to be complex, and only when one of the exponents is odd, so the processes which were already handled keep a real one and are untouched. Looking the two operands up the other way round now uses the inverse factor. With a sign the two were the same thing, which is why it went unnoticed. g g > t t~ g g g g goes from 115 858 lines to 72 580, its GET_JAMP block from 4.35 MB to 1.51 MB, its -O2 compile from 835 s to 30 s and its -O2 run time from 0.0667 ms to 0.0403 ms per call. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…imisation The routine which rebuilds the definitions from one recipe per orbit and the plain list of their operands drive the very same loop; they only differ in how that list reaches memory. Rebuilding it needs the amplitude permutations in the source, and those are namp x nperm numbers whatever the number of definitions, so on a large process they cost more than the list they replace: on g g > 6g the tables of INIT_JAMP are 5.7 MB against the 6.1 MB of the list itself. Write the list out instead, and keep the orbit equivariant optimisation which is what makes it short in the first place. The recipes are still there behind jamp_emit for comparison. Dropping the walk also makes the definitions free to be reordered, which the recipes could not be without the generated code losing track of them. The ones introduced together use none of each other, so inside each of those groups they are sorted by the factor in front of the second operand and the loop runs over the ones adding it, then the ones subtracting it, then the rest. Only that last group multiplies, and for a process whose factors are all signs it is empty and the factor array disappears. g g > 6g: matrix.f 19.4 MB -> 16.2 MB, GET_JAMP block 9.20 MB -> 6.10 MB. g g > 5g: 1.30 MB -> 1.20 MB, and 0.0081 -> 0.0078 ms per call at -O2, then 0.0072 with the loop split. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An orbit can only be taken as a whole, so the orbit rounds stop while the JAMP lines still hold a good many terms: 7215 amplitude references on g g > 5g and 75 615 on g g > 6g, against 1456 and 10 096 for the plain scan. Long straight line expressions are what -O2 spends its time on, which is why the table emission was slower to compile than the plain scan on exactly those two processes and not on g g > t t~ g g g g, where the orbit rounds run to the end. Let the plain scan finish the job once they stall. Its sub-expressions are not orbits of anything, which keeps them out of the recipes, but the table emission does not care: there a definition costs three numbers of DATA and one indirect add, while a term left in a line costs a term of source and a direct add. So it is shorter, quicker to compile and fewer operations all at once. g g > 5g: 9990 definitions and 6480 terms become 12657 and 720, matrix.f 1.20 MB -> 0.82 MB, the GET_JAMP block 0.60 -> 0.23 MB, -O2 compile 16.5 -> 4.9 s and -O2 run time 0.00723 -> 0.00666 ms per call. The levels the loop runs over are now read off the operands rather than off the rounds, so that whatever the scan adds at the end lands where it belongs. myjamp_count was only ever set by get_JAMP_lines, so the scan crashed when optimise_jamp was called on its own, as the tests do. It is a class attribute now. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…esses Taking whole orbits only pays once there is enough of them to share. On a small matrix it can come out longer than the plain scan, which is free to take whatever it likes: on g g > t t~ g the orbit rounds ask for 46 additions where the scan asks for 43, and that showed up as the generated GET_JAMP going from 5.5 to 8.4 microseconds a call. Small matrices are cheap to optimise, so rather than guess where the turn is, run both up to twenty thousand entries and keep the shorter. Above that only the orbit version runs: it wins by a wide margin on everything that big, and the plain scan is the slow one there. g g > t t~ g is back to what the scan gives. The processes which were already ahead are untouched: g g > t t~ g g 317 -> 304 additions, g g > 3g 108 -> 106, g g > 4g 1083 -> 931, and g g > 5g still 12657 definitions. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Same optimisation, same gain: on g g > 5g the definitions go from 22221 to 12657 and matrix1_orig.f from 3.61 MB to 1.64 MB. The generated code has the same shape as before, only fewer lines of it, so the stored comparison files are untouched. The definitions are written out here rather than read from a table. Madevent rewrites the matrix element for helicity recycling, which is on by default, and that rewriting builds the routine from a template of its own where AMP is indexed by helicity as AMP(NCOMB,NGRAPHS), carrying over the JAMP lines and nothing else. The table wants the definitions at the end of a one dimensional AMP so that the loop takes both operands from one array, which that shape cannot give, and the DATA statements would be dropped on the way. So jamp_tables_allowed says no for madevent and only the optimisation is shared. get_JAMP_lines_split_order now passes the matrix element down as the place to read the color basis from: it hands over one list of color amplitudes per order, which does not carry one. It only asks for the orbit version when there is a single order to compute, since one set of definitions is all the template holds. Checked with a full run of g g > g g g: every number of every results.dat is identical to what main gives, only the timings differ. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The definitions could not be read from a table under helicity recycling: that
rewriting indexes AMP by helicity as AMP(NCOMB,NGRAPHS), so the definitions
cannot sit at the end of it and share one index with the amplitudes.
Read the amplitudes of the current helicity into a buffer first, and run the
definitions over that. The buffer is one dimensional whatever AMP looks like,
so the table is the same one the standalone code uses, and the definitions stay
one dimensional instead of gaining a helicity index they would never use.
The recycler needs no change. It rewrites every AMP( it finds in the JAMP block
to AMP( K,, and the only one left there is the gather itself, which is exactly
what wants the helicity index:
AMPBUF(ITMP) = AMP(ITMP) -> AMPBUF(ITMP) = AMP( K,ITMP)
The declarations reach the rewritten file because the template it is built from
comes out of the same replace_dict as the matrix element.
The copy is bought back. On g g > 5g, with the JAMP block run once per helicity
for 80 helicities, the definitions written out and read from AMP(K,i) strided
take 1.5377 ms a call against 1.2312 for the gather, and gfortran -O2 takes
436 s against 0.17 s.
The DATA statements have their own loop variable now, the templates already use
I for something else.
Checked with a full run of g g > g g g, tables forced on: every number of every
results.dat is what main gives, only the timings differ.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ix on it Reversing every color basis element maps the basis onto itself, and for a pure gluon process the color coefficients of a line and of its reverse differ by one overall sign: JAMP[reverse(i)] = (-1)^n JAMP[i]. Half the color flows therefore carry nothing of their own, and |M|^2 can be summed over one line per pair with a folded color matrix instead of over every line. That is half the color flows and a quarter of the color sum, which is 78% of the run time of g g > 6g. This is the reading half of it, not yet used by the generated code: - reverse_immutable, the reversal of a basis key - get_jamp_reflection, which reads the sign off the color coefficients rather than assuming it, so a process where the relation does not hold gets None. It holds for pure gluons (+1 for an even number, -1 for an odd one) and not for a quark line, where reversing does not commute with the fermion flow. - jamp_folded_color_matrix, C'[a][b] summed over the two lines of each pair with their signs. Checked against the unfolded contraction for g g > g g, g g > 3g and g g > 4g: same |M|^2 to all digits, same denominator, exactly half the lines and no self-paired line in any of them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Reversing a color flow gives another flow of the basis whose amplitude is the same up to one overall sign, so half of them carry nothing of their own. Sum over one per pair against a color matrix folded onto those, and the color sum does a quarter of the work. It is 78% of the run time of g g > 6g, so this is the largest single lever there is short of changing the color basis itself. GET_MATRIX gathers the flows it keeps out of JAMP through COLREP and runs the same loop over them. JAMP stays the full length, so the color flow decomposition, GET_INTER and the density matrix see exactly what they saw before. Two ways of getting the folded matrix to the generated code: - sign +1: reversing commutes with permuting the color indices, so a permutation of the flows is also a permutation of the pairs and the folded matrix still has one line per orbit. It is rebuilt at run time as before. - sign -1: a permutation may send a flow onto its own partner, which flips its weight, and the rebuilt form has nowhere to put that sign. The folded matrix is written out instead, which is affordable while it stays small (g g > 5g folds to 360 lines, 65k entries) and folding is declined above color_fold_max_written. The sign is read off the color coefficients rather than assumed, so a quark line, where reversing does not commute with the fermion flow, simply does not fold: g g > t t~ g g and u u~ > u u~ g g are untouched and their |M|^2 is unchanged to the last digit. g g > 5g: GET_MATRIX 0.3459 -> 0.0844 ms a call, 4.10x, the 4x expected. |M|^2 6.6739867626784571E-007 -> ...550E-007, 15 digits: the folded sum regroups 720^2 terms into 360^2 with larger coefficients, so it rounds differently. g g > 4g and g g > 3g are unchanged to the last digit. test_generate_helas_diagrams_gg_gg checks the emitted color matrix against a hand written 6x6, which is the unfolded one; it now asks for that explicitly. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The optimisation is given only the flows the color sum keeps, and the rest are one assignment each. Sharing means the definitions do not halve with the lines -- they are reused across flows -- but there are fewer of them: g g > 4g goes from 795 to 645, 19% off, with |M|^2 unchanged to the last digit. Only with sign +1. There both lines of a pair weigh the same, so a permutation of the lines is a permutation of the pairs with nothing attached and what is left of the matrix is still symmetric enough for the orbit version to work on it. With sign -1 a permutation may send a kept line onto its partner, which flips the weight for that pair alone; the matrix is then not symmetric in the form optimise_jamp wants, the symmetry is lost and with it the orbit optimisation and the table emission, which costs far more than the halving saves -- g g > 5g went to 41175 lines from 21575 when I tried. get_jamp_halving is where that is decided, and g g > 3g and g g > 5g keep every flow. The copies are still written: only the color sum skips those flows, while the color flow decomposition, GET_INTER and the density matrix read the whole array. g g > 3g, g g > t t~ g g and u u~ > u u~ g g are untouched, all 122 tests pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This reverts commit 48ab222.
get_color_data_lines sits on ProcessExporterFortran and is shared by madevent, madweight and the fks exporters, but only the standalone template was taught to sum over NCOLORFOLD. So those three were handed a folded matrix and kept looping over every line: g g > g g g wrote a 78 entry triangle for the 12 pairs into CF(300) and then read all 300, the last 222 never set. jamp_fold moves to the standalone exporter, where the template agrees. The others generate exactly what they did before the folding went in. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The three madevent templates can now sum |M|^2 over one color flow per reversal pair: NCOLORFOLD, a JFOLD gathered per split order, and the sum reading %(color_fold_array)s, which is JAMP when there is no folding so the generated code is unchanged in that case -- g g > g g and g g > g g g both come back at 4.454e+08 and 3.677e+07 pb, the numbers from before the patch. It is off because the folded numbers are wrong: 2.672e+09 against 4.454e+08 on g g > g g and 3.971e+09 against 3.677e+07 on g g > g g g. Both compile and run, so this is not a missing declaration, and matrix1_optim.f -- what actually runs under helicity recycling -- is rewritten by hel_recycle out of the template, which is where I would look first: the gather sits inside the helicity loop and reads JAMP(COLREP(ICF),ICFSO), and the recycler rewrites JAMP references. get_color_data_lines folds for whichever exporter asks, so template_matrix1.f had to be taught the same sum: it comes from the group hel template and feeds the recycler. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
hel_recycle scrapes the lines between the jamp block and the color sum into amp2_lines, starting at the first ENDDO after the jamp lines and stopping at a six space indented "DO I = 1, NCOLOR". The gather loop put an ENDDO in front of that start, and the folded sum sits two levels in, so the whole sum was scraped and emitted a second time inside the sde_strat branch -- without its DENOM line, which that scrape drops on purpose. The second copy overwrote TS(K), so |M|^2 came out DENOM times too large: 6x on g g > g g, 108x on g g > g g g, each exactly the DENOM of that process. A vector subscript gathers in one statement and leaves the scrape alone. The recycler is unchanged: it was reading what it was written to read. g g > g g and g g > g g g both come back at 4.454e+08 and 3.677e+07 pb with the folding on and the recycler doing its work, the same to every digit quoted, uncertainties included, as the unfolded runs. TS(K) is zeroed once and keeps its DENOM. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The color matrix is real and symmetric and MATRIX is real, so the imaginary
part of the sum is thrown away as it stands. Splitting JAMP into its two parts
computes what is kept and nothing else:
Re( sum_j C_ij JAMP_j conj(JAMP_i) ) = (sum_j C_ij ReJ_j) ReJ_i
+ (sum_j C_ij ImJ_j) ImJ_i
term by term, so it needs no cancellation argument and holds for every process,
not only the ones whose coefficients are real. The gather fills two REAL*8
arrays instead of one COMPLEX*16, and the four accumulators become eight.
g g > g g g, g g > g g g g and g g > t t~ g g are unchanged to the last digit.
u u~ > u u~ g g moves in the last digit, the same 1 ulp it already carried
before this.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both are the case where the color basis has no reversal symmetry, so the folding declines and the only change is the name of the loop bound: NCOLORFOLD equals NCOLOR, the gather is empty and the sum still reads JAMP. testIO_export_matrix_element_v4_standalone still fails, from 9fe3443 rather than from here: that commit moved the standalone sum to real arithmetic and left the reference on the complex form. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This reverts commit 9fe3443.
A pure gluon process picks up one factor of i per f^abc, the same for every term, so the coefficient matrix is either real throughout or imaginary throughout. Dividing that factor out leaves whole numbers for optimise_jamp to compare and hash, with nothing widened to complex, and it goes back onto the JAMP coefficients afterwards -- the definitions hold ratios, which it cancels out of. g g > 6g generates a byte identical matrix.f either way. It is worth less than it looks: the color flow step goes from 136s to 131s, under 4%, on the process where it costs the most. A quark line mixes the two phases and gets the old path. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
DSYMM takes the whole helicity set as one right hand side instead of walking the color matrix once per helicity. On a synthetic g g > 6g sized problem the kernel goes from 2.485 to 0.069 ms per helicity, and that is with DSYMM doing about twice the arithmetic: it reads the full symmetric matrix where the scalar loop reads a packed triangle. DSYMM is real, so the two parts of JAMP go through separately -- the color matrix is real and symmetric, so nothing else is needed. The full matrix is built once from the written triangle, whose off diagonal is doubled because the scalar sum walks it once; halving it there is exact, the entries are integers. Only in the steady state: USERHEL unset, the good helicities settled and no polarization selection. Discovery and filtering keep the scalar path, which also stays as the fallback whenever BLAS is not taken. blas defaults to None, meaning take it when DSYMM links and the folded matrix has at least blas_min_ncolor rows -- below that the call is not worth setting up. g g > 5g agrees with the scalar build to 1.4 ulp, which is what reassociat- ing the sum costs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The guard needed a bracket around the condition it extends, and writing that bracket in the template changed the file even with BLAZ off. Both halves are placeholders now, so a build without BLAS writes what it always wrote. g g > 6g, 100 calls a run over three runs: 232.5s scalar against 121.7s with BLAS, 1.91x, and that still counts the first twenty calls, which run scalar while the good helicities settle. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
It arrived here through a git add -A, in both halves: the placeholder in the three madevent templates and the value for it in the madevent exporter. It is someone else's fix and it comes back with their branch. Taking it out puts the madevent references back in agreement -- all 122 tests pass, and no stored reference was touched. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A compressed color matrix is handed over as a description and rebuilt at run time by INIT_CF. The standalone template reads CF out of a common block that the routine fills; madevent had CF as a local array, so nothing reached it and the sum ran over whatever was on the stack -- g g > g g g g came out at -9.1e+09 pb against 8.344e+06. It only showed up above the size where the matrix is compressed at all, which is why the small processes looked fine. The block is named per matrix subroutine. A grouped directory links several of them into one executable, and one /color_matrix/ between them would quietly hand every subprocess the matrix of whichever ran first. get_color_init_routine takes a suffix for that, empty for standalone, so its output is unchanged. g g > g g g g now gives 8.344e+06 +- 4.095e+04 pb folded and unfolded alike, the number main gives; u u~ > u u~ g, which is grouped, is unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
They came in with the merge and they are not mine to change. The color flow folding does change what the exporters write, so the three comparisons fail until someone who owns them regenerates them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
madevent read its color matrix out of a local array while INIT_CF filled a common block, so above the size where the matrix is compressed the sum ran over the stack. Needed here: this branch carries both the compressed color matrix and the madevent JAMP tables.
They were put back untouched when the color folding came in, since that session did not own them, and have been failing since. Three things have changed what the exporters write and all three are accounted for here: the folded color sum brings NCOLORFOLD and runs the sum over the folded flows; madevent now emits INIT_CF and calls it, since its color matrix is rebuilt at run time like the standalone one; and TMP_JAMP moved onto a declaration of its own so the table emission can leave it out. |M|^2 is unchanged: g g > g g g g and g g > g g g g g agree with main to one ulp standalone, and g g > g g g g through madevent gives 8.458e+06 +- 1.794e+04 pb either way, with every number of every results.dat identical.
The color matrix does not depend on the helicity, so the helicities are the columns of a single right hand side and the whole sum is two DSYMM calls instead of one triangular loop per helicity. DSYMM is real, so the two parts of JAMP go through separately; the matrix is real symmetric, so the two products add term by term with nothing left over. Where the matrix element is rewritten for helicity recycling it already walks every helicity inside one call, so there the batch is the loop it is already running. Everywhere else MATRIX is one helicity at a time and SMATRIX holds the loop, so one sweep fills the per-helicity values and the loop only reads them back -- nothing is computed twice and AMP2/JAMP2 still add up once. The file the recycling rewriter reads is left alone. Taken when a BLAS links and the color basis is big enough to pay for the call, and switchable at run time with the hidden blas_color_sum run card parameter. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ment
p p > l+ l- came out at ~538 pb with apply_flavor_grouping=False against
~1336 pb with it on. The cause is not broken_sym (l+ and l- are distinct,
so no identical-particle factor applies) but a missed source of flavor
multiplicity in the mg7 exporter.
A matrix element carries flavor multiplicity in two independent places:
* merged legs (pdg 81/82/...), enumerated by
get_external_flavors_with_iden -- the apply_flavor_grouping=True case;
* the several processes mapped onto one matrix element when their
matrix elements are identical -- the grouping-off case, where
u u~ > e+ e-, u u~ > mu+ mu-, c c~ > e+ e- and c c~ > mu+ mu- all
share a single matrix element.
madevent walks both in get_leshouche_lines, and its leshouche.inc lists
all 16 channels of p p > l+ l-. export_mg7 built its flavor list from
get_external_flavors_with_iden alone, so with grouping off -- where all
the multiplicity sits in the processes list -- subprocesses.json recorded
one channel per subprocess instead of four: c/s dropped from the initial
state and mu+ mu- dropped entirely. The exporter was already inconsistent
with itself, since pdg_color_types in the same function does loop over
all processes and duly listed +-4 and +-13.
Add HelasMatrixElement.get_flavor_pdg_combinations implementing madevent's
rule once, have get_leshouche_lines consume it (behaviour preserving, its
IDUP numbering is untouched), and expand the mg7 flavor list through it.
Verified:
- ungrouped p p > l+ l- now 1335.30 +- 1.65 against grouped
1336.14 +- 1.36, i.e. 0.28 sigma; subprocesses.json carries all 16
channels, structured as the grouped output already was;
- madevent output byte-identical across all four
apply_flavor_grouping/group_subprocesses combinations (full
SubProcesses tree diff, not just leshouche.inc);
- 46/46 IOTests pass;
- test_flavor_grouping_consistency_mg7, test_generation_from_file_1_mg7,
test_group_subprocess_mg7, test_e_e_collision_mg7 all pass.
Also add the acceptancetest_mg7_flavor_grouping CI job, which runs
test_flavor_grouping_consistency_mg7 -- previously not covered by any
workflow, which is why this stayed unnoticed.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
test_flavor_grouping_consistency intermittently reported "Incompatible cross-sections ... (1e-99 +- 0)", which reads like a physics disagreement but means the run produced nothing at all. Two things conspired to hide that. First, the guards could not catch it. Each configuration checked its own precision with err/(cross+1e-99); when a run fails, cross and error are both 0, so that is 0/1e-99 = 0 < 0.05 and it passes. The explicit zero-checks were dead code for the same reason: two of them tested val == 0 *after* adding 1e-99, and the MLM one asserted val > 0 on the padded value. So a zero always slipped through to the pairwise comparison, where it produced a misleading message. Assert on the raw cross instead, in all four affected tests (the three madevent ones and the mg7 counterpart), naming the run directory. Second, when such a run does fail the cause is usually a ZeroResult swallowed by nice_error_handling, which only logs a warning -- invisible at the CRITICAL level the tests run at -- leaving cross/error at 0. The new assertion turns that from a puzzle into a one-line diagnosis. Also fix multiple_try: after exhausting its retries it ended with a bare `raise` outside the except block, so instead of re-raising the real error it raised "RuntimeError: No active exception to reraise", discarding the only information explaining the failure (and doing so after several sleeps). Hit while debugging the above, where it masked an AttributeError from load_results_db. Raise my_error, which keeps the traceback; the non-__debug__ branch is unchanged. Verified: the new guard fires on cross=0 and stays silent on a good run; multiple_try now surfaces the original exception with its traceback under __debug__ and the wrapped message under -O. test_flavor_grouping_ consistency, _width, _mg7 and test_generation_from_file_1_mg7 all pass. Note this makes a future occurrence self-diagnosing; it does not fix the underlying transient, which I could not reproduce (~28 configurations, idle and under load). A likely contributor is gen_ximprove.get_helicity classifying a subprocess as having no phase space purely from `if stdout:` on ./gensym, without ever checking its return code. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
For a process whose colour structure is purely adjoint, the Jacobi identity lets any colour factor be written on the (n-2)! half-ladder structures instead of the (n-1)! traces. The reduction walks each diagram's colour string as a tree and expands every subtree hanging off the leg-1-to-leg-n spine as a nested commutator of adjoint generators, which costs 2^(n-2) per diagram. Selected by 'set color_basis auto|trace|ddm', with auto turning it on for the exporters which declare support_ddm_color_basis. Fortran standalone and the grouped madevent exporter do; everything else keeps the trace basis, and any process which is not fully adjoint falls back on its own. Two things this needed: - f.complex_conjugate now returns f untouched. The inherited implementation reverses the indices, which is right for Tr but flips the sign of the totally antisymmetric f. It was latent because full_simplify removes every f before the trace-basis colour matrix is built; with DDM the basis elements are products of (n-2) f's, so |M|^2 came out negative for odd n. - ColorMatrix.build_matrix_ddm expands each ladder once onto its 2^(n-2) traces and assembles the entry from cached trace-trace products. Contracting two ladders head on turns every one of the 2(n-2) f's into a pair of traces: gg>gggg spent 17.4s there against 0.089s in the trace basis, and gg>ggggg never finished. It is now 0.12s and 1.6s. madevent keeps the trace basis for the colour flows: the colour sum runs on the DDM structures while JAMP2 is filled from trace JAMPs rebuilt out of the DDM ones through the Kleiss-Kuijf relations, 3840 terms instead of 231840 at n=7. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Standalone is the mode used to time the matrix element, so it has to carry what madevent carries. It now builds the trace basis next to the DDM one and emits GET_JAMPF, the Kleiss-Kuijf rebuild of the trace JAMPs, together with the JAMP2 accumulation that madevent does for its color flow probabilities. The batched BLAS color sum takes its own branch in SMATRIX and never enters MATRIX, so the flow JAMPs are computed there as well; without that the BLAS timings were measuring a matrix element madevent could not use. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…oni-gluon-basis-69e1d7 # Conflicts: # madgraph/iolibs/export_v4.py # madgraph/iolibs/template_files/matrix_madevent_group_v4_hel.inc
The C++ backend now does for a multi-gluon process what the fortran ones do: the color sum runs on the (n-2)! Del Duca-Dixon-Maltoni structures, while the color flow is still picked among the (n-1)! trace ones, whose jamps are rebuilt from the DDM ones through the Kleiss-Kuijf relations rather than from the amplitudes. The two counts were the same number before, so ncolor did for both. They are now separate: ncolor stays the color sum, ncolor_flow is what a color flow is picked among, and every buffer holding jamp2 -- the accumulation in calculate_jamps, the jamp2_sv arrays, the colAllJamp2s super-buffer, the selection walking icolamp -- moves to ncolor_flow. Without a DDM basis ncolor_flow is ncolor and the generated code is what it was. g g > g g g g: ncolor 120 -> 24, CPPProcess.cc 512 -> 269 kB, |M|^2 unchanged to the last digit (1.5929925846563478e-04 against 1.5929925846563475e-04). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Reversing a color flow gives the same flow back up to one overall sign, so half of them carry nothing of their own and |M|^2 can be summed over one flow per pair against a color matrix folded onto them. The fortran output has done this for a while; the C++/GPU one summed the full ncolor x ncolor triangle. The reflection and the folded matrix are not reimplemented: the four methods the fortran exporters were carrying move to a ColorReflectionFolding mixin, which madmatrix inherits. The only thing the two backends disagree on is when a folding is worth taking, now a jamp_fold_worthwhile hook: fortran can rebuild a sign +1 matrix at run time and so refuses to write a large sign -1 one, while madmatrix always writes the matrix out and folded it is a quarter of the size. What is generated is the folded matrix alone, over ncolorfold flows, plus the colorFoldRep table saying which flow of each pair is kept. color_sum_cpu and color_sum_kernel gather through that table and sum over ncolorfold. jamp2 is untouched: color selection needs every flow separately, and it is accumulated from the full local jamp before the gather. The BLAS color sum is NOT folded - folding it means compacting allJamps, which the jamps are not written in - but it does not get a second matrix either: it multiplies the folded matrix spread back over the ncolor flows, with the dropped rows and columns left at zero, so it cannot drift away from what the kernel computes. Compacting the jamps is left to do. Measured on the color sum alone (FPTYPE=d, one call, arm64): g g > g g g 24 -> 12 flows, sign -1: 95 ns -> 28 ns (3.4x) g g > g g g g 120 -> 60 flows, sign +1: 3.6 us -> 0.71 us (5.1x) Validated against unfolded builds over 512 phase space points: max relative difference 6.5e-16 (5 ulp) for g g > g g g and 1.2e-15 (8 ulp) for g g > g g g g -- a re-grouped sum cannot be bit-identical, and a wrong sign is not a 1e-16 effect. Both agree with the fortran ./check to the same accuracy, and the folded matrices match the exact rational fold of the unfolded ones on every entry. u u~ > u u~ g, which does not fold, stays bit-identical. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The C++/cudacpp backend wrote one 'jamp_sv[i] += c*amp_sv[0]' per (color flow, amplitude) pair, which for g g > g g g g is 8160 statements and half a megabyte of CPPProcess.cc. The fortran side has had a common sub-expression pass over the same coefficient matrix for a long time (optimise_jamp), but the two shared no code. Everything about that pass which is not printing now lives in madgraph/iolibs/jamp_optimiser.py, mixed into both exporters. The fortran output is unchanged, byte for byte. The C++ emitter cannot print the definitions the way fortran does: fortran keeps AMP(NGRAPHS) to read the amplitudes back from, while here every amplitude passes through the single slot amp_sv[0] and is gone by the next diagram. So they are accumulated instead -- each amplitude is added into the definition that uses it while it is still there, and the definitions built from other definitions are emitted at the amplitude they become ready after. That keeps the temporaries down to one array next to jamp_sv, where storing every amplitude would have cost 510 more cxtype_sv for g g > g g g g. g g > g g g g: 8160 -> 1371 jamp statements, 512 -> 274 kB, CPPProcess.cc compiles in 0.96s instead of 2.04s and the matrix elements come out 35% faster. |M|^2 agrees with the expanded output to 7e-16 (the optimisation reassociates the sums, exactly as the fortran one does), and the mg7 cross-section acceptance tests are unchanged. --jamp_optim=False recovers the expanded output. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
run_matrix_mode passed the fptype copy of the masses to classic_rambo::get_momenta, whose signature takes std::vector<double>. At FPTYPE=d/m fptype is double so this compiled by accident; at FPTYPE=f it is float and there is no vector-to-vector conversion, so any generated process failed to build in single precision. Pass massesD, the double vector the masses are read into, instead. Both vectors are already in scope and both are still needed: the fptype copy keeps feeding RamboSamplingKernelHost. Widening the get_momenta signature would have been wrong here, since classic_rambo works in double throughout precisely so it reproduces the phase-space point of the Fortran/C++ 'check' drivers -- feeding it float-rounded masses would perturb that reference. Checked on g g > g g g with a clean rebuild at each precision: all three build, the phase-space point is identical, and the matrix elements agree to ~2 ulp of float (f vs d 2.7e-7 relative, m vs d 1.2e-8). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The cuBLAS path in color_sum.cc was GPU only: the C++ side had a hand written triangular SIMD loop and nothing else. The color matrix does not depend on the helicity, so on C++ too the jamps of every good helicity of one event page are the columns of one right hand side and the whole color sum becomes two SYMM calls, real and imaginary parts apart. Measured first. On the trace basis the C++ color sum is a real share of the matrix element: 31% for g g > g g g g and 14% for g g > t t~ g g g at ncolor=120 (sampling and a gutted-color-sum difference agree), against 5% to 11% at ncolor=24. Accelerate's SYMM does about twice the flops of the triangular loop and still wins by 9x on its own, so the color sum drops to 7% of the matrix element and the whole thing runs 1.30x faster on g g > g g g g, 1.13x on g g > t t~ g g g. Taken exactly where the Fortran color sum takes it - the same DSYMM probe, the same ncolor threshold - so with BLAS off, or below the threshold, color_sum.cc and CPPProcess.cc are character for character the files written before any of this existed. Above it both paths are written out and CPPBLAS=hasNoBlas still builds the old one. The batched call bypasses nothing: the jamp2 sums for the color selection, the numerators and the denominators all stay inside calculate_jamps and still happen once per helicity. What color_sum_cpu did on top of the sum was feed MEs_ighel, the running sum over helicities the event by event choice of helicity reads, so the batched call rebuilds those itself. Mixed precision is supported rather than switched off, through SSYMM. |M|^2 agrees with a non-BLAS build to one ulp (some points bit identical); the summation order differs, so it is not bit for bit. The C++ to Fortran gap at the same phase space point is 6.5e-15, forty times larger. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…-gluon-basis-69e1d7 The fold branch was cut before the batched BLAS color sum landed, so the two had never met. Folding redimensions colorMatrix and colorDenom onto the ncolorfold representatives, while color_sum_blas.inc still walked them over ncolor: reading past the end of a constexpr array is not a constant expression, so g g > g g g g -- the one process here where both fire -- did not compile. The batched sum now runs on the folded flows, gathering the representative of each reversal pair exactly as the scalar kernel does (the folded matrix already carries what the dropped partner contributes), with the scratch, the SYMM leading dimension and the final dot product all on ncolorfold. Where the basis does not fold ncolorfold is ncolor and the identity map leaves it as it was. g g > g g g g, ncolor=120 folded onto 60 with BLAS on: 1.5929925846563470e-04 against 1.5929925846563475e-04 unfolded, which is the re-grouped sum's few ulp and not a sign error. u u~ > u u~ g, which does not fold, is untouched.
The previous commit folded color_sum_cpu and color_sum_kernel but left the BLAS color sum doing the full ncolor x ncolor GEMM: it multiplied the folded matrix spread back over the ncolor flows, with the dropped rows and columns at zero. Correct, and it could not drift away from the kernel, but roughly four times more arithmetic than needed. Folding it means compacting the jamps, which the representative flows being non-contiguous makes a gather and not a copy. convertD2F_Jamps becomes gatherFold_Jamps: it reads through colorFoldRep and is compiled in every precision mode, not just the mixed one, since the conversion to fptype2 is a no-op when the two types are the same. The gathered jamps go in a second buffer carved out of ghelAllBlasTmp, and all four BLAS calls then run over ncolorfold against s_pNormalizedColorMatrixFold2. With nothing left reading it, UnfoldedNormalizedColorMatrix and s_pNormalizedColorMatrix2 are dropped, and with them hostColorFoldRep, whose only consumer they were. ncolorfold lived only in color_sum.cc and the buffer sizing needs it elsewhere, so it is exported as CPPProcess::ncolorfold from process_class.inc and color_sum.cc reads it back from there, as it already did for ncolor. A process which does not fold (ncolorfold == ncolor) would otherwise pay for this twice, in a scratch buffer of double the size and in an identity gather which copies the jamps onto themselves. So the gather and its buffer are taken only when the gather is not the identity or there is a conversion to do; those processes run exactly the code they ran before. The rule and the buffer size are stated once, in blasColorSumTmpSize, which both the allocation in MatrixElementKernels.cc and the reset in color_sum_gpu now call instead of carrying a copy of the formula each. Scratch buffer, nGoodHel * nevt fptype2 per unit: g g > g g g 24 -> 12 flows: 48 -> 48 (d/f), 97 -> 49 (mixed) g g > g g g g 120 -> 60 flows: 240 -> 240 (d/f), 481 -> 241 (mixed) u u~ > u u~ g no folding: 8 -> 8 (d/f), 17 -> 17 (mixed) NB: NOT validated on a GPU - no CUDA toolkit was available, so nvcc has never seen these sources and neither has cuBLAS. What was done instead is to compile the GPU branch of the generated color_sum.cc with a host shim (empty __device__ and __global__, emulated gridDim/threadIdx, a reference column-major GEMM honouring the op/ld/stride semantics) and run color_sum_gpu down both of its paths on the same jamps. Over 3 processes x 3 precision modes the BLAS color sum agrees with the kernel one (which is what HASBLAS=hasNoBlas runs) to 6.1e-16 in double and 6.0e-07 in float, agrees with a long double reference to the same accuracy, and reproduces the unfolded BLAS color sum it replaces exactly. A canary past the end of the scratch buffer is untouched in all nine. Four negative controls turn that test red: a gather ignoring colorFoldRep, a batched GEMM stride left at ncolor, an under-sized buffer, and the no-gather shortcut taken on a process which does fold. The CPU color sum is untouched and ./check_sa.exe is bit-identical for g g > g g g. colorFoldRep itself has not been through nvcc either - it came in with the previous commit. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…uon-basis-69e1d7 Both sides had pulled a mixin out of ProcessExporterFortran since the fork: this one moved the JAMP optimisation into jamp_optimiser.JampOptimiser so the C++ writer could reach it, the colour folding merge moved the reversal-pair work into ColorReflectionFolding. The class now carries both, and the hunk where the two overlapped keeps neither copy: jamp_global_phase comes from JampOptimiser and get_jamp_folding from ColorReflectionFolding. Rebuilt every combination that now exists, since the last merge of this kind compiled clean and then did not: g g > g g g trace ncol=24 fold=12 blas=0 1.8740711159594317e-02 g g > g g g ddm ncol=6 fold=6 blas=0 1.8740711159594321e-02 g g > g g g g trace ncol=120 fold=60 blas=1 1.5929925846563470e-04 g g > g g g g ddm ncol=24 fold=24 blas=0 1.5929925846563470e-04 u u~ > u u~ g (no fold, no blas) 1.5036487888197825e-03 all within a few ulp of the values before the merge, and the quark process bit-identical. CPPProcess.cc for g g > g g g g on the trace basis goes from 512 to 276 kB with 1373 shared sub-expressions; on the DDM basis the JAMPs are already few enough that it only goes 269 to 250 kB.
…i-gluon-basis-69e1d7 Folds the cuBLAS color sum onto one flow per reversal pair, the GPU half of the folding whose CPU half came in with 4bb8e20. Two additive conflicts: CPPProcess now carries both ncolor_flow (what a color flow is picked among) and ncolorfold (what |M|^2 is summed over), and color_sum.h keeps both new includes. The CPU paths were rebuilt and are unchanged: g g > g g g g on the trace basis folded 120 onto 60 with the host BLAS gives 1.5929925846563470e-04, on the DDM basis 1.5929925846563470e-04, and g g > g g g 1.8740711159594317e-02, all as before this merge. The GPU path itself is emitted but not compiled: there is no CUDA or HIP toolchain on this machine.
The fortran exporters look for the shared JAMP sub-expressions by whole orbits
of the permutations leaving the color basis invariant, not one at a time
(jamp_orbit). The C++/cudacpp writer had none of it: it ran the plain greedy
scan the previous commit gave it and stopped there.
Everything about that search which is not printing now lives in
jamp_optimiser.JampOptimiser next to the plain scan -- get_jamp_symmetry,
optimise_jamp_equivariant, optimise_jamp_best and their helpers, about 360
lines. Two hooks are left for the backends, because it is the emission which
decides what is usable: jamp_orbit_allowed (now also reading --jamp_orbit at
output time) and jamp_greedy_tail_enabled, which stays
'jamp_emit == tables' for fortran since INIT_JAMP cannot rebuild
sub-expressions that are not orbits of anything. What prints the result --
jamp_orbit_recipes, jamp_orbit_tables, get_jamp_decl_lines,
get_jamp_init_routine, jamp_gather, namp_dim -- stays in export_v4. The
fortran output is unchanged, byte for byte.
optimise_jamp_best earns its keep straight away: on g g > g g g g over the DDM
basis the plain scan is the shorter of the two (392 definitions against 399)
and is kept, so that build comes out identical.
Definitions and CPPProcess.cc, with the search off and on:
g g > g g g g ddm 392 -> 392 250 220 -> 250 220 B
g g > g g g g trace 951 -> 795 276 077 -> 266 886 B
u u~ > u u~ ggg trace 756 -> 732 487 327 -> 482 790 B
g g > t t~ ggg trace 3030 -> 2535 950 946 -> 899 241 B
g g > g g g g g ddm 8524 -> 7560 2 392 147 -> 2 289 786 B
The file moves less than the definitions do because the color flows are only
13 to 28% of CPPProcess.cc here, against 78% of the fortran matrix.f: the rest
is the HELAS calls. |M|^2 is unchanged to the last digit on every process
tested, and agrees with the fortran ./check at the same phase-space point to
6e-15 or better -- the pre-existing gap between the two backends, not something
this adds. Throughput is unchanged within the noise of the machine; the color
flows are not what these processes spend their time on.
No table emission. The stage which could be driven by a table is the one
building definitions out of other definitions, plus the color flows built out
of them: both sides are arrays that already exist. The amplitudes are not --
each one passes through the single slot amp_sv[0] and is gone by the next
diagram, so the capture stage is at least one written statement per amplitude
whatever the search does, and giving it an array to index would cost
nb_amp * cxtype_sv per SIMD page and per GPU thread (232 kB for
g g > g g g g g). The table-drivable part peaks at 4515 entries on the largest
process generated here, and a microbenchmark on the real cxtype_sv puts the
crossover at about ten thousand -- twice the fortran figure of five thousand:
200 defs lines 0.071 us table 0.128 us ratio 1.81
1000 defs lines 0.355 us table 0.645 us ratio 1.82
4515 defs lines 2.329 us table 3.549 us ratio 1.52
8000 defs lines 7.678 us table 8.213 us ratio 1.07
16000 defs lines 18.74 us table 18.34 us ratio 0.98
32000 defs lines 48.14 us table 37.08 us ratio 0.77
So tables would be 1.5x slower on the biggest thing the backend has, for 8% of
the source. jamp_emit stays fortran-only.
--jamp_orbit=False recovers the plain scan, on either backend.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Contributor
Author
|
Full unit suite on this branch: 917 tests, 3 red.
So one genuine reference to regenerate, and it is a reference rather than a
|
The mg7 exporter can run its color sum on the (n-2)! Del Duca-Dixon-Maltoni
structures while a color flow is still picked among the (n-1)! trace ones, but
set_channels_colors_map built active_colors -- the mask which becomes icolamp --
by walking self.color_basis, i.e. the basis of the color sum. The C++ selection
walks that mask over ncolor_flow entries, so a DDM build declared
icolamp[nconfig][ncolor] and read it as icolamp[nconfig][ncolor_flow]: for
g g > g g a [3][2] array read 6 wide, which reinterprets the neighbouring rows
and then runs off the end of the object. Reading icolamp the way the selection
does gave {1,2,4,5}, {0,2,3}, {0,1,5} for the three configs instead of
{0,1,3,5}, {1,2,3,4}, {0,2,4,5}.
|M|^2 does not depend on any of this, which is why the DDM port passed every
check run so far. Events do: over 3 x 100k g g > g g events the color flow
fractions came out (9.35, 9.38, 31.38, 9.15, 31.30, 9.44)% on the trace basis
against (10.54, 1.39, 32.32, 12.83, 32.46, 10.46)% on the DDM one -- 134 sigma
on the second flow -- at an unchanged cross section (1.2 sigma), since a wrong
color selection never moves the weight.
Give the exporter a color_flow_basis next to its color_basis and build
active_colors, the color_flows table and both nb_color values from it. The two
bases now write the same icolamp, and the same processes agree flow by flow
over 3 x 100k events each: g g > g g worst 2.0 sigma over 6 flows (xsec 0.1),
g g > g g g worst 1.9 sigma over 24 (xsec 0.5), and u u~ > g g, which never
leaves the trace basis, generates identical source and agrees at 0.1 sigma.
The GPU select_col_and_diag kernel had the same ncolor/ncolor_flow mismatch and
is fixed the same way; there is no CUDA toolchain here, so that path is
generated but not compiled.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
do_multiline breaks a long statement at the last space which fits in 72
characters. When the only such space is the statement's own indentation the
first chunk is all blanks, so the emitted physical line is empty and the
continuation which follows it is attached by the compiler to the *previous*
statement.
g g > g g g on the DDM color basis is the first process to hit it: the
Kleiss-Kuijf lines rebuilding the trace JAMPs are long and contain no space
of their own, and matrix1_optim.f came out as
JAMPF(1,1)=+2D0*(+IMAG1*JAMP(6,1))
$ JAMPF(2,1)=+2D0*(-IMAG1*JAMP(3,1)-...)
which gfortran rejects with "Unclassifiable statement", so the madevent run
dies in the compilation step. Break mid-token in that case, which is what the
no-space-at-all branch next to it already does.
The condition only fires where the old code emitted a blank chunk, so any
matrix element which built before is byte for byte what it was; of the arms
generated here only g g > g g g on the DDM basis contained one.
With this, g g > g g g builds on both color bases and the two agree exactly:
same cross section, and all 10000 color columns of all 2000 events identical
at equal seed.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Makes the colour sum of a matrix element cheaper, along five independent
lines that all land on the same code. Together they take
g g > 6gfromabout 2.2 s to 0.554 s per phase-space point, and none of them changes
a physics result.
What is in here
A (n-2)! colour basis for multi-gluon processes. The Jacobi identity lets
any fully adjoint colour factor be written on the (n-2)! Del Duca-Dixon-Maltoni
half-ladders instead of the (n-1)! traces. New
set color_basis auto|trace|ddm;autotakes it for the exporters that declaresupport_ddm_color_basis.Processes that are not fully adjoint fall back on their own, per amplitude.
Kleiss-Kuijf colour flows. A colour flow still has to be picked among the
(n-1)! trace structures, but those JAMPs need not be built from the amplitudes:
J_trace(rho) = sum_sigma c_{sigma rho} J_DDM(sigma), using the expansioncoefficients the colour matrix already computes. At
g g > 5gthat is 3840terms instead of 231840, 3.0 us against 22.4 us per helicity. LHE colour
flows are unchanged.
Colour-reflection folding. Reversing a basis element maps the basis onto
itself, so |M|^2 can be summed over one flow per reversal pair against a folded
matrix. 3.4x on the colour sum at
g g > 3g, 5.1x atg g > 4g.A batched BLAS colour sum. The colour matrix does not depend on helicity,
so every good helicity is a column of one right-hand side and the sum is two
SYMM calls instead of a triangular loop per helicity. Fortran and C++ host,
plus cuBLAS on the device.
The same treatment for the madmatrix (C++/cudacpp) backend, which had none
of it: the DDM basis and the KK flows, the folding on both the CPU and cuBLAS
paths, the host BLAS sum, and the JAMP common-subexpression and orbit-equivariant
searches, which are now shared with the Fortran exporters through the new
madgraph/iolibs/jamp_optimiser.pyrather than duplicated.Numbers
g g > 6g, standalone, this machine, per phase-space point:CPPProcess.ccforg g > g g g gon the trace basis: 512 -> 267 kB.Worth saying plainly: the DDM basis no longer wins in steady state once
folding is in (0.554 against 0.544, trace 2% ahead). It still buys 1.82x on the
first ten points and 1.87x in generation, and it is what makes the first point
cheap.
GET_AMPis 85% of the remaining time, which nothing here touches.Correctness
f.complex_conjugatewas silently wrong. The inherited implementationreverses the indices, which is right for
Trbut flips the sign of thetotally antisymmetric
f. It was latent becausefull_simplifyremoves everyfbefore the trace-basis colour matrix is built; with the DDM basis theelements are products of (n-2) f's, so |M|^2 came out negative for an odd
number of gluons. Fixed.
four phase-space points for
g g > 5g;g g > 6gis bit-identical.g g > g gthrough madevent with the same seed gives the identical crosssection and identical colour flows on all 8000 particles of 2000 events.
them at
g g > g g g g.a sign error; processes that do not fold stay bit-identical.
Notes for review
toolchain on the machine this was developed on. That covers the cuBLAS folding
and the device side of the colour sum.
tests/unit_tests/iolibsdrift, because thegenerated code changed (
NCOLOR_FLOWin the madevent template, and the JAMPdefinitions). They are left for the maintainer to regenerate.
does not yet know about: the batch sets
BLASDONEand switches off the scalarloop, so it would evaluate every good helicity where the loop did one per
mirror pair. That is a 1.86x regression on the merge of the two, not on
this branch, which has no de-duplication. The fix is written and sits on
claude/colour-merge-retimewith the merge of Crossing symmetry #49 and Sum the four gluon current into the cubic one carrying the same colour factor #57.🤖 Generated with Claude Code