diff --git a/.github/workflows/IOtest.yml b/.github/workflows/IOtest.yml index 722a9495a..e6da9bb28 100644 --- a/.github/workflows/IOtest.yml +++ b/.github/workflows/IOtest.yml @@ -64,6 +64,17 @@ jobs: ./tests/test_manager.py testIO_test_pptt_fks_loonly testIO_test_pptt_fksreal testIO_test_pptt_fksrealew testIO_test_ppw_fksall testIO_test_tdecay_fksreal testIO_test_wprod_fksew -t0 + # FD gauge output reference files (fortran DHELAS + madmatrix HelAmps) + IOtest_fd_gauge: + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v5 + - name: test FD gauge output IO tests + run: | + cd $GITHUB_WORKSPACE + ./tests/test_manager.py testIO_FDgauge_standalone_fortran testIO_FDgauge_madmatrix -pA -t0 + + # C++ exporter IOTests (from unittest_20) IOtest_cpp_write: runs-on: ubuntu-24.04 diff --git a/.github/workflows/acceptancetest_mg7.yml b/.github/workflows/acceptancetest_mg7.yml index bc334e43b..a168a1be6 100644 --- a/.github/workflows/acceptancetest_mg7.yml +++ b/.github/workflows/acceptancetest_mg7.yml @@ -58,6 +58,23 @@ jobs: cd $GITHUB_WORKSPACE ./tests/test_manager.py test_mg7_ufo_aloha -pA -t0 -l INFO + # FD gauge: the madmatrix and the fortran standalone must compute the same + # value, and the scalar and vectorised madmatrix backends must agree with each + # other (the FD wavefunctions are written once per backend in helas_fd.h) + acceptancetest_mg7_fd_gauge: + needs: build_madspace + runs-on: ubuntu-24.04 + if: github.event_name == 'push' || github.event.pull_request.head.repo.fork == true + steps: + - uses: actions/checkout@v4 + - uses: ./.github/actions/checkout_mg5 + - uses: ./.github/actions/install_madspace + - uses: ./.github/actions/restore-pip-cache + - name: test FD gauge madmatrix vs fortran + run: | + cd $GITHUB_WORKSPACE + ./tests/test_manager.py test_standalone_mg7_fd_vs_fortran test_standalone_mg7_fd_simd_lanes -pA -t0 -l INFO + acceptancetest_mg7_output_directory: needs: build_madspace runs-on: ubuntu-24.04 diff --git a/aloha/aloha_writers.py b/aloha/aloha_writers.py index c9ec6da6c..4aaf70775 100755 --- a/aloha/aloha_writers.py +++ b/aloha/aloha_writers.py @@ -76,6 +76,13 @@ def __init__(self, abstract_routine, dirpath, options=None): #initialize global helper routine self.declaration = Declaration_list() self.options = options if options else {} + # a writer emitting one term of an assembled multi-structure routine + # (see write_combined_parts) works on the wavefunctions under the name + # the assembled routine gives them, adds its own coupling and adds up + # into the output instead of assigning it. + self.wf_rename = {} + self.coup_name = None + self.combined_part = False def pass_to_HELAS(self, indices, start=0): @@ -191,7 +198,24 @@ def get_P_sign(self, index): def get_foot_txt(self): """Prototype for language specific footer""" return '' - + + def has_fd_propagator(self): + """Does this routine have to apply the FD gauge propagator factor on + the wavefunction it builds? (language independent)""" + + if aloha.unitary_gauge != 3 or not self.offshell or 'P1N' in self.tag: + return False + return self.particles[self.outgoing-1].startswith(('V','S')) + + def rename_wf(self, name): + """Name under which a wavefunction is written. It is its own name + except in a term of an assembled multi-structure routine, where all the + structures share one name per leg: in FD gauge a leg can be a vector in + one structure and its Goldstone in the next one, and both are the same + argument.""" + + return self.wf_rename.get(name, name) + def define_argument_list(self, couplings=None): """define a list with the string of object required as incoming argument""" @@ -582,7 +606,11 @@ def get_declaration_txt(self): if type.startswith('list'): type = type[5:] #determine the size of the list - if name[0] in ['F', 'V', 'S', 'T', 'R']: + if name.startswith('FD'): + # FD gauge: 5-momentum and gauge direction of the inlined + # propagator factor (no wavefunction, despite the F) + out.write(' %s %s(0:4)\n' % (self.type2def[type], name)) + elif name[0] in ['F', 'V', 'S', 'T', 'R']: # All wavefunctions (inputs and outputs) are now passed and # built as type(aloha) / type(aloha2d), regardless of # loop_mode. This keeps the body code (which uses %W / %P @@ -718,8 +746,10 @@ def get_momenta_txt(self): if self.offshell and aloha.unitary_gauge == 3: # FD gauge type = self.particles[self.outgoing-1] - if type in ["S","V"]: - out.write(" %(type)s%(out)s %% W(:) = CZERO \n" % {'type': type, 'out':self.outgoing}) + if type in ["S","V"]: + out.write(" %(type)s%(out)s %% W(:) = CZERO \n" % {'type': type, 'out':self.outgoing}) + if self.declaration.is_used('FDQ'): + out.write(self.get_fd_gauge_txt()) # Returning result return out.getvalue() @@ -754,6 +784,19 @@ def get_coupling_def(self): return out.getvalue() + # a merged multiple coupling routine has one matrix coupling per Lorentz + # structure (MCOUP1, MCOUP2, ...). They all belong to the same vertex + # (same particles) so they share the flavour pairing, but a given + # flavour can be missing from some of them (vanishing coupling). The + # first coupling defines the flavour mapping of the routine, the others + # only contribute when they agree with it. + couplings = [name for ftype, name in self.declaration + if name.startswith('COUP')] + couplings.sort(key=lambda x: int(x[4:]) if x[4:] else 0) + if not couplings: + couplings = ['COUP'] + main = couplings[0] + if self.outgoing == 0 or self.particles[self.outgoing-1] not in ['F']: if not self.outgoing: fail = "VERTEX = (0d0,0d0)" @@ -763,32 +806,41 @@ def get_coupling_def(self): out.write(' flv_index1 = F1 %flv_index\n') out.write(' flv_index2 = F2 %flv_index\n') out.write(' if(flv_index1.eq.0.or.flv_index2.eq.0)then \n %s\n return\nendif\n' % fail) - out.write(' if(MCOUP %% PARTNER(flv_index1).ne.flv_index2)then \n %s\n return\n endif\n' %fail) - else: - incoming = [i+1 for i in range(len(self.particles)) if i+1 != self.outgoing and self.particles[self.outgoing-1] == 'F'][0] - if incoming %2 == 1: - outgoing = self.outgoing - out.write(' flv_index%i = F%i %%flv_index\n' % (incoming, incoming)) - out.write(' if(flv_index%i.eq.0)then\n' %(incoming)) - out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) - out.write(' flv_index2 = MCOUP %% PARTNER(FLV_INDEX%i)\n' %(incoming)) - out.write(' if(flv_index2.eq.0)then\n') - out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) - out.write(' F%i %% FLV_INDEX = FLV_INDEX2\n' % outgoing) - else: - outgoing = self.outgoing - out.write(' flv_index%i = F%i %%flv_index\n' % (incoming,incoming)) - out.write(' if(flv_index%i.eq.0)then\n' %(incoming)) - out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) - out.write(' flv_index1 = MCOUP %% PARTNER2(FLV_INDEX%i)\n' %(incoming)) - out.write(' if(flv_index1.eq.0)then\n') - out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) - out.write(' F%i %% FLV_INDEX = FLV_INDEX1\n' % outgoing) - - for ftype, name in self.declaration: - if name.startswith('COUP'): + out.write(' if(%s)then \n %s\n return\n endif\n' % + ('.and.'.join('M%s %% PARTNER(flv_index1).ne.flv_index2' % name + for name in couplings), fail)) + # the flavour of the outgoing particle is fixed by the incoming ones + for name in couplings: + if name == main: + out.write(' %s = M%s %% VAL(flv_index1) %% p \n' % (name, name)) + else: + out.write(' %s = (0d0,0d0)\n' % name) + out.write(' if(M%s %% PARTNER(flv_index1).eq.flv_index2) %s = M%s %% VAL(flv_index1) %% p \n' + % (name, name, name)) + return out.getvalue() + + incoming = [i+1 for i in range(len(self.particles)) if i+1 != self.outgoing and self.particles[self.outgoing-1] == 'F'][0] + outgoing = self.outgoing + # PARTNER maps the odd fermion onto the even one, PARTNER2 the opposite + partner = 'PARTNER' if incoming % 2 == 1 else 'PARTNER2' + # index of the fermion built by this routine + out_index = 2 if incoming % 2 == 1 else 1 + out.write(' flv_index%i = F%i %%flv_index\n' % (incoming, incoming)) + out.write(' if(flv_index%i.eq.0)then\n' %(incoming)) + out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) + out.write(' flv_index%i = M%s %% %s(FLV_INDEX%i)\n' % (out_index, main, partner, incoming)) + out.write(' if(flv_index%i.eq.0)then\n' % out_index) + out.write(' F%i %% W(:) = (0d0,0d0)\n F%i %% flv_index = 0 \n return\n endif\n' %(outgoing, outgoing)) + out.write(' F%i %% FLV_INDEX = FLV_INDEX%i\n' % (outgoing, out_index)) + + for name in couplings: + if name == main: out.write(' %s = M%s %% VAL(flv_index1) %% p \n' % (name, name)) - return out.getvalue() + else: + out.write(' %s = (0d0,0d0)\n' % name) + out.write(' if(M%s %% %s(flv_index%i).eq.flv_index%i) %s = M%s %% VAL(flv_index1) %% p \n' + % (name, partner, incoming, out_index, name, name)) + return out.getvalue() def get_one_momenta_def(self, i, strfile): @@ -816,7 +868,11 @@ def shift_indices(self, match): if aloha.unitary_gauge == 3 and match.group('var').startswith('S'): shift += 4 # In FD gauge Scalar indices goes to 5 (not 1) # to complement the vector 1-4 - return '%s %% W(%s)' % (match.group('var'), int(match.group('num'))+ shift) + # the slot is fixed by the spin this structure sees on the leg, the + # name by the routine the code is written in (rename_wf) + return '%s %% W(%s)' % (self.rename_wf(match.group('var')), + int(match.group('num'))+ shift) + def change_var_format(self, name): """Formatting the variable name to Fortran format""" @@ -841,7 +897,7 @@ def change_var_format(self, name): # the scratch variables stay real for performance. if decla.startswith('P'): vtype = 'complex' if aloha.loop_mode else 'double' - self.declaration.add(('list_%s' % vtype, decla)) + self.declaration.add(('list_%s' % vtype, self.rename_wf(decla))) else: self.declaration.add((name.type, name)) name = re.sub(r'(?P\w*)_(?P\d+)$', self.shift_indices , name) @@ -924,20 +980,30 @@ def sort_fct(a, b): numerator = self.routine.expr - if not 'Coup(1)' in self.routine.infostr: + if self.coup_name: + # one term of an assembled routine: this structure has its own + # coupling among COUP1, COUP2, ... + coup_name = self.coup_name + elif not 'Coup(1)' in self.routine.infostr: coup_name = 'COUP' else: + # the couplings are already inside the expression (merged routine) coup_name = '%s' % self.change_number_format(1) + has_coup = coup_name != self.change_number_format(1) if not self.offshell: - if coup_name == 'COUP': + if has_coup: formatted = self.write_obj(numerator.get_rep([0])) - if formatted.startswith(('+','-')): - out.write(' vertex = COUP*(%s)\n' % formatted) + if not formatted.startswith(('+','-')): + formatted = '%s*%s' % (coup_name, formatted) else: - out.write(' vertex = COUP*%s\n' % formatted) + formatted = '%s*(%s)' % (coup_name, formatted) else: - out.write(' vertex = %s\n' % self.write_obj(numerator.get_rep([0]))) + formatted = self.write_obj(numerator.get_rep([0])) + if self.combined_part: + out.write(' vertex = vertex + %s\n' % formatted) + else: + out.write(' vertex = %s\n' % formatted) else: OffShellParticle = '%s%d' % (self.particles[self.offshell-1],\ self.offshell) @@ -976,8 +1042,8 @@ def sort_fct(a, b): else: coeff = '%(COUP)s*' % {'COUP': coup_name} else: - if coup_name == 'COUP': - coeff = 'COUP*' + if has_coup: + coeff = '%s*' % coup_name else: coeff = '' to_order = {} @@ -999,9 +1065,14 @@ def sort_fct(a, b): # type(aloha) and we write into %W which is 1-indexed for # Lorentz components only. shift -= self.momentum_size + # the slot is fixed by the spin of this structure, the name by + # the routine the code is written in; a term of an assembled + # routine adds up into a slot several structures can feed + slot = '%s%%W(%d)' % (self.rename_wf(self.outname), + self.pass_to_HELAS(ind)+shift) to_order[self.pass_to_HELAS(ind)] = \ - ' %s%%W(%d)= %s%s\n' % (self.outname, self.pass_to_HELAS(ind)+shift, - coeff, formatted) + ' %s= %s%s%s\n' % (slot, '%s+' % slot if self.combined_part + else '', coeff, formatted) key = list(to_order.keys()) key.sort() for i in key: @@ -1031,6 +1102,13 @@ def sort_fct(a, b): # retry when removing the useless part. return self.define_expression() + if self.has_fd_propagator() and not self.combined_part: + # the FD gauge propagator factor is part of the wavefunction this + # routine builds, not a post-treatment. A term of an assembled + # routine gets it once, from the assembler, on the total. + self.declare_fd_propagator() + txt += self.get_fd_propagator_txt() + return txt def define_symmetry(self, new_nb, couplings=None): @@ -1042,27 +1120,210 @@ def define_symmetry(self, new_nb, couplings=None): # (self.get_header_txt(new_name, couplings), self.name, ','.join(arguments)) def get_foot_txt(self, combine=False): - text = ' ' - - if not combine and aloha.unitary_gauge == 3: # FD gauge - if self.outgoing and 'P1N' not in self.tag: - name = self.particles[self.outgoing-1] - if name.startswith(('V','S')): - # need to be smarter for Higgs - text += 'CALL MULTIPLY_PROPAGATOR_FACTOR(%(name)s%(i)s,%(mass)s%(i)s, %(name)s%(i)s)\n' %\ - {'name':name, 'mass': 'M%s' % name[1:], 'i': self.outgoing } + text = ' ' - - text += 'end\n\n' + text += 'end\n\n' return text + def get_fd_gauge_txt(self): + """FD gauge: the part of the propagator factor that only depends on the + momentum of the outgoing wavefunction, i.e. the 5-momentum q (the mass + sits in its 5th component) and the gauge direction n. + + This is the inlined counterpart of the head of + multiply_propagator_factor (aloha_functions_fd.f); it is written with + the momenta so that everything it defines stays out of the way of the + wavefunction dependent part (see get_fd_propagator_txt).""" + + out = StringIO() + outname = self.outname + out.write(' FDQ(0:3) = -%s %% P(:)\n' % outname) + out.write(' FDQ(4) = -CI*M%s\n' % self.outgoing) + out.write(' CALL DEFINE_GAUGE_DIR(FDQ, FDN)\n') + out.write(' FDNQ = %s\n' % '-'.join(['FDN(%d)*DBLE(FDQ(%d))' % (i,i) + for i in range(4)])) + return out.getvalue() + + def get_fd_propagator_txt(self): + """FD gauge: the wavefunction dependent part of the propagator factor. + + The 5 components built by the routine are projected on the physical + gauge: w -> w - q * js1 - n * js2. Inlining it (instead of calling + multiply_propagator_factor) keeps the routine a single expression: the + momentum only part is hoisted with the momenta, and what is left here + is linear in the wavefunction.""" + + out = StringIO() + outname = self.outname + # number of Lorentz components of the outgoing wavefunction: in FD + # gauge a vector and its Goldstone share the same 5 slots + size = self.type_to_size[self.particles[self.outgoing-1]] - 2 + # js1 and js2 contract the wavefunction with n and q (the 5th component + # of q carries the mass, hence the conjugation) + out.write(' FDJS1 = (%s)/FDNQ\n' % '-'.join( + ['FDN(%d)*%s%%W(%d)' % (i, outname, i+1) for i in range(4)])) + out.write(' FDJS2 = (%s-DCONJG(FDQ(4))*%s%%W(5))/FDNQ\n' % ( + '-'.join(['FDQ(%d)*%s%%W(%d)' % (i, outname, i+1) for i in range(4)]), + outname)) + for i in range(size): + out.write(' %(o)s%%W(%(k)d) = %(o)s%%W(%(k)d)-FDQ(%(i)d)*FDJS1' + '-FDN(%(i)d)*FDJS2\n' % {'o': outname, 'k': i+1, 'i': i}) + return out.getvalue() + + def declare_fd_propagator(self): + """Variables of the inlined FD gauge propagator factor. The mass it + needs is already an argument of any offshell routine.""" + + self.declaration.add(('list_complex', 'FDQ')) + self.declaration.add(('list_double', 'FDN')) + self.declaration.add(('double', 'FDNQ')) + self.declaration.add(('complex', 'FDJS1')) + self.declaration.add(('complex', 'FDJS2')) + def write_combined(self, lor_names, mode='self', offshell=None): - """Write routine for combine ALOHA call (more than one coupling)""" - + """Write routine for combine ALOHA call (more than one coupling). + + The routine is a genuine merge of the Lorentz structures: its body is + the sum of the structures, each multiplied by its own coupling, and it + is written exactly like a single coupling routine (contracted + temporaries, single propagator DENOM, one assignment per component). + Structures that do not act on the same spins can not share a single + expression; they are assembled term by term instead (write_combined_parts). + When neither is possible, fall back on a wrapper calling each single + structure routine in turn. + """ + + if offshell is None: + offshell = self.offshell + + merged = None + if offshell == self.offshell and not os.environ.get('MG_ALOHA_COMBINE_WRAPPER') \ + and hasattr(self.routine, 'get_combined_routine'): + merged = self.routine.get_combined_routine(lor_names) + if merged is None: + parts = self.routine.get_combined_routines(lor_names) + if parts: + return self.write_combined_parts(parts, lor_names, mode, + offshell) + if merged is None: + return self.write_combined_wrapper(lor_names, mode, offshell) + + name = combine_name(self.routine.name, lor_names, offshell, self.tag) + # the merged routine is written by a dedicated writer (fresh + # declarations) through the standard path: it is a routine like any + # other one, only its name and its file are those of the combination. + writer = self.__class__(merged, None, options=self.options) + writer.name = name + # the caller passes one coupling per Lorentz structure: declare them all + # up-front so that the signature does not depend on a structure + # surviving the algebra for this outgoing particle. + for i in range(len(lor_names) + 1): + writer.declaration.add(('complex', 'COUP%s' % (i+1))) + text = writer.write(mode=mode) + + if self.out_path: + fsock = self.writer(self.out_path, 'a') + commentstring = 'This File is Automatically generated by ALOHA \n' + commentstring += 'The process calculated in this file is: \n' + # indent: the merged expression starts with "Coup(1)" and a line + # starting with a 'C' is taken as an already formatted comment + commentstring += ' ' + merged.infostr + '\n' + fsock.write_comments(commentstring) + fsock.writelines(text) + return text + + def write_combined_parts(self, routines, lor_names, mode='self', offshell=None): + """Write a combined routine whose Lorentz structures do not act on the + same spins, by assembling them term by term. + + This is the FD gauge case: a leg is a massive vector in one structure + and its Goldstone in the next one, but both are the same wavefunction + (components 1-4 and 5 of one object), so no single ALOHA expression can + carry them. The routine is still a single one: the momenta, the + declarations and the propagator factor are shared, and each structure + adds its own coupling times its own expression into the slots it feeds. + """ + + name = combine_name(self.routine.name, lor_names, offshell, self.tag) + writers = [self.__class__(routine, None, options=self.options) + for routine in routines] + main = writers[0] + # one name per leg -- the one the first structure gives it + legs = ['%s%d' % (spin, i+1) for i, spin in enumerate(main.particles)] + + bodies = [] + for i, writer in enumerate(writers): + writer.name = name + writer.coup_name = 'COUP%s' % (i+1) + writer.combined_part = True + for j, spin in enumerate(writer.particles): + wf = '%s%d' % (spin, j+1) + if wf != legs[j]: + writer.wf_rename[wf] = legs[j] + bodies.append(writer.define_expression()) + + # everything the terms need has to be declared and, for the momenta, + # defined once for all of them + for writer in writers[1:]: + for entry in writer.declaration: + main.declaration.add(entry) + # one coupling per structure, whatever the algebra kept + for i in range(len(lor_names) + 1): + main.declaration.add(('complex', 'COUP%s' % (i+1))) + if main.has_fd_propagator(): + main.declare_fd_propagator() + + text = StringIO() + text.write(main.get_header_txt(name=name)) + text.write(main.get_declaration_txt()) + text.write(main.get_momenta_txt()) + text.write(main.get_coupling_def()) + # the terms add up into the output, so it has to start from zero (in FD + # gauge the momenta already reset the wavefunction) + zero = '(%s,%s)' % (self.change_number_format(0), + self.change_number_format(0)) + if not offshell: + text.write(' vertex = %s\n' % zero) + elif not (aloha.unitary_gauge == 3 and + main.particles[main.outgoing-1] in ['S', 'V']): + text.write(' %s%%W(:) = %s\n' % (main.outname, zero)) + # the structures share the contraction cache, so the same TMP/FCT is + # often built by several of them: keep the first definition only + seen = set() + for body in bodies: + for line in body.splitlines(True): + definition = line.strip() + if re.match(r'^(?:TMP|FCT)\d+\s*=', definition): + if definition in seen: + continue + seen.add(definition) + text.write(line) + if main.has_fd_propagator(): + text.write(main.get_fd_propagator_txt()) + text.write(main.get_foot_txt()) + + text = text.getvalue() + if self.out_path: + fsock = self.writer(self.out_path, 'a') + commentstring = 'This File is Automatically generated by ALOHA \n' + commentstring += 'The process calculated in this file is: \n' + commentstring += '\n'.join(' Coup(%s) * (%s)' % (i+1, routine.infostr) + for i, routine in enumerate(routines)) + commentstring += '\n' + fsock.write_comments(commentstring) + fsock.writelines(text) + return text + + def write_combined_wrapper(self, lor_names, mode='self', offshell=None): + """Write a combine ALOHA routine (more than one coupling) as a wrapper + calling the routine of each Lorentz structure and summing the results. + Only used when the merged routine can not be built (see write_combined). + """ + # Set some usefull command if offshell is None: sym = 1 - offshell = self.offshell + offshell = self.offshell else: sym = None name = combine_name(self.routine.name, lor_names, offshell, self.tag) diff --git a/aloha/create_aloha.py b/aloha/create_aloha.py index 1d992d22d..e5e90a92c 100755 --- a/aloha/create_aloha.py +++ b/aloha/create_aloha.py @@ -59,7 +59,11 @@ class AbstractRoutine(object): """ store the result of the computation of Helicity Routine this is use for storing and passing to writer """ - + + # builders of the merged multiple coupling routines, shared by all the + # outgoing of a given set of Lorentz structures (see get_combined_routine) + combined_builder = {} + def __init__(self, expr, outgoing, spins, name, infostr, model, denom=None): """ store the information """ @@ -84,10 +88,134 @@ def add_symmetry(self, outgoing): def add_combine(self, lor_list): """add a combine rule """ - + if lor_list not in self.combined: self.combined.append(lor_list) + def get_combined_routine(self, lor_names): + """Return the AbstractRoutine associated to the merged structure + Coup(1) * + Coup(i+1) * + i.e. the routine that a writer can output as a single subroutine (one + coupling argument per structure) instead of a wrapper calling each + single structure routine in turn. + + The expression is built exactly like the one of a single coupling + routine, so the momenta, the propagator denominator and the temporary + variables are shared by all the structures. + Return None if such a merge is not possible (unknown Lorentz structure, + structures acting on different spins, loop routine, ...). In that case + the caller has to fall back on the wrapper form. + """ + + if not hasattr(self, 'combined_routine'): + self.combined_routine = {} + + key = tuple(lor_names) + if key not in self.combined_routine: + try: + self.combined_routine[key] = self.compute_combined_routine(lor_names) + except Exception as error: + logger.debug('can not merge the routines %s (%s): %s', + self.name, ','.join(lor_names), error) + self.combined_routine[key] = None + + routine = self.combined_routine[key] + if routine is not None: + # the same merged routine is written once per tag set (the MP pass + # adds the 'MP' tag on the fly) -> keep the tag in sync. + routine.tag = list(self.tag) + return routine + + def compute_combined_routine(self, lor_names): + """Compute the merged routine of get_combined_routine (no caching)""" + + if self.model is None: + return None + if any(t.startswith('L') for t in self.tag): + # loop routines have their own (already explicit) combine mechanism + return None + + l_lorentz = [getattr(self.model.lorentz, name) + for name in (self.name,) + tuple(lor_names)] + if any(lor.spins != l_lorentz[0].spins for lor in l_lorentz[1:]): + return None + + conjg = tuple(int(t[1:]) for t in self.tag if t.startswith('C')) + # the model is part of the key: two models can define different + # structures under the same Lorentz name + key = (self.model, tuple(lor.name for lor in l_lorentz), conjg) + if key in self.combined_builder: + builder = self.combined_builder[key] + else: + builder = CombineRoutineBuilder(l_lorentz, self.model) + if conjg: + builder = builder.define_conjugate_builder(conjg) + # the kernel is computed once and re-used for each outgoing + self.combined_builder[key] = builder + + routine = builder.compute_routine(self.outgoing, list(self.tag)) + # the merge does not change which particles are identical + routine.symmetries = list(self.symmetries) + routine.tag = list(self.tag) + return routine + + def get_combined_routines(self, lor_names): + """Return the routines to sum -- one per Lorentz structure, this one + first -- when the structures do not act on the same spins, so that + get_combined_routine can not build a single expression for them. + + This is the FD gauge situation: a massive vector and its Goldstone are + the very same wavefunction (components 1-4 and 5 of one object), so the + structures of a vertex mix V and S on a given leg. They still write + into one output and read one wavefunction per leg, so a writer can + assemble them into a single routine. + Return None when even that is not possible. + """ + + if not hasattr(self, 'combined_routines'): + self.combined_routines = {} + + key = tuple(lor_names) + if key not in self.combined_routines: + try: + self.combined_routines[key] = self.compute_combined_routines(lor_names) + except Exception as error: + logger.debug('can not assemble the routines %s (%s): %s', + self.name, ','.join(lor_names), error) + self.combined_routines[key] = None + + routines = self.combined_routines[key] + if routines is not None: + for routine in routines: + routine.tag = list(self.tag) + return routines + + def compute_combined_routines(self, lor_names): + """Compute the routines of get_combined_routines (no caching)""" + + if self.model is None: + return None + if any(t.startswith('L') for t in self.tag): + return None + + conjg = tuple(int(t[1:]) for t in self.tag if t.startswith('C')) + routines = [self] + for name in lor_names: + lorentz = getattr(self.model.lorentz, name) + if len(lorentz.spins) != len(self.spins): + # not the same legs: nothing to assemble them into + return None + key = (self.model, name, conjg) + if key in self.combined_builder: + builder = self.combined_builder[key] + else: + builder = AbstractRoutineBuilder(lorentz, self.model) + if conjg: + builder = builder.define_conjugate_builder(conjg) + self.combined_builder[key] = builder + routines.append(builder.compute_routine(self.outgoing, list(self.tag))) + return routines + def write(self, output_dir, language='Fortran', mode='self', combine=True, options=None, **opt): """ write the content of the object """ # Set loop_mode based on the tag of this specific routine @@ -699,7 +827,15 @@ def __init__(self, l_lorentz, model=None): self.outgoing = None self.lorentz_expr = [] for i, lor in enumerate(l_lorentz): - self.lorentz_expr.append( 'Coup(%s) * (%s)' % (i+1, lor.structure)) + structure = lor.structure + # AbstractRoutineBuilder.__init__ only expanded the formfactors of + # l_lorentz[0], and that expansion is overwritten here -> redo it + # for each structure entering the combination. + if getattr(lor, 'formfactors', None): + for formf in lor.formfactors: + pat = re.compile(r'\b%s\b' % formf.name) + structure = pat.sub('(%s)' % formf.value, structure) + self.lorentz_expr.append( 'Coup(%s) * (%s)' % (i+1, structure)) self.lorentz_expr = ' + '.join(self.lorentz_expr) self.routine_kernel = None self.contracted = {} diff --git a/aloha/template_files/madmatrix/helas_fd.h b/aloha/template_files/madmatrix/helas_fd.h index 891e28b31..a56ad00de 100644 --- a/aloha/template_files/madmatrix/helas_fd.h +++ b/aloha/template_files/madmatrix/helas_fd.h @@ -493,14 +493,13 @@ // FD gauge const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1 + fptype_sv{ 0 } ); -#ifdef MGONGPU_CPPSIMD - fptype_sv nA[5]; - fptype_sv nB[5]; -#endif fptype_sv n[5]; fptype_sv nk; - const fptype_sv zero{0.}; - const fptype_sv one{1.}; + // NB: broadcast to every SIMD lane. 'fptype_sv one{1.}' sets the first + // lane only (the others are zero filled), which left n, and then nk, at + // zero in all lanes but the first one + const fptype_sv zero = 0. + fptype_sv{ 0 }; + const fptype_sv one = 1. + fptype_sv{ 0 }; if( vmass != 0. ) { @@ -569,7 +568,7 @@ w[1] = cxmake( -vmass/nk * n[1], zero ); w[2] = cxmake( -vmass/nk * n[2], zero ); w[3] = cxmake( -vmass/nk * n[3], zero ); - w[4] = static_cast(nsv)*cI; + w[4] = -static_cast(nsv)*cI; // as in fortran vxxxxx (vc%W(5) = -nsv*ci) and in the SIMD branch below } #else @@ -604,37 +603,33 @@ w[2] = cxternary( mask, vcA_4, cxternary( maskB, vcB1_4, vcB2_4 ) ); w[3] = cxternary( mask, vcA_5, vcB_5 ); - //FD gauge - //branch A - nA[0] = fpternary( pvec0 >= zero , one , -one); - nA[1] = -pvec1/pp; - nA[2] = -pvec2/pp; - nA[3] = -pvec3/pp; - nA[4] = zero; - - //branch B - nB[0] = nA[0]; - nB[1] = zero; - nB[2] = zero; - nB[3] = -nA[0]; - - const fptype_sv b_A = fpternary(pp > zero, one , zero); - const fptype_sv b_B = fpternary(pp <= zero , one , zero); - - n[0] = nA[0]*b_A + nB[0]*b_B; - n[1] = nA[1]*b_A + nB[1]*b_B; - n[2] = nA[2]*b_A + nB[2]*b_B; - n[3] = nA[3]*b_A + nB[3]*b_B; - n[4] = nA[4]; + //FD gauge: same two branches as the scalar code above, selected lane by + //lane. The division uses ppDENOM (=1 where pp==0) so that the lanes that + //do not take it are not poisoned: a select discards the other value, but + //a nan surviving an arithmetic combination (nan*0 is nan) would not be. + const bool_v maskFD = ( pp > zero ); + n[0] = fpternary( pvec0 >= zero , one , -one ); + n[1] = fpternary( maskFD, -pvec1 / ppDENOM, zero ); + n[2] = fpternary( maskFD, -pvec2 / ppDENOM, zero ); + n[3] = fpternary( maskFD, -pvec3 / ppDENOM, -n[0] ); + n[4] = zero; nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; - const bool_v mask3 = { (abs(nhel) == 1 ? 1 : 0) }; // first element replicated - w[0] = cxternary( mask3, w[0], cxmake( -vmass/nk * n[0], zero)); - w[1] = cxternary( mask3, w[1], cxmake( -vmass/nk * n[1], zero)); - w[2] = cxternary( mask3, w[2], cxmake( -vmass/nk * n[2], zero)); - w[3] = cxternary( mask3, w[3], cxmake( -vmass/nk * n[3], zero)); - w[4] = cxternary( mask3, cxzero_sv(), -static_cast(nsv)*cI); + // nhel is a scalar: no need for a per-lane mask here (and a bool_v built + // from a single value would only set the first lane) + if ( abs(nhel) == 1 ) + { + w[4] = cxzero_sv(); + } + else + { + w[0] = cxmake( -vmass/nk * n[0], zero ); + w[1] = cxmake( -vmass/nk * n[1], zero ); + w[2] = cxmake( -vmass/nk * n[2], zero ); + w[3] = cxmake( -vmass/nk * n[3], zero ); + w[4] = -static_cast(nsv)*cI; + } #endif } else @@ -1025,7 +1020,7 @@ n[0] = fpternary( q[0].real() >= 0.f , one , -one ); n[1] = zero; n[2] = zero; - n[3] = fpternary( q[0].real() >= 0.f , -one , one); //possible error in Fortran + n[3] = fpternary( q[0].real() >= 0.f , -one , one); // -sign(q0), as in fortran and python define_gauge_dir n[4] = zero; } #else @@ -1034,7 +1029,7 @@ n[0] = fpternary( q[0].real() >= 0.f , one , -one); n[1] = fpternary( qsign , -q[1].real() / qabs , zero ); n[2] = fpternary( qsign , -q[2].real() / qabs , zero ); - n[3] = fpternary( qsign , -q[3].real() / qabs , fpternary( q[0].real() >= 0.f , one , -one)); + n[3] = fpternary( qsign , -q[3].real() / qabs , fpternary( q[0].real() >= 0.f , -one , one)); // same gauge as the branch above n[4] = zero; #endif } diff --git a/madmatrix/model_handling.py b/madmatrix/model_handling.py index b4d1bac78..5805c7d4e 100644 --- a/madmatrix/model_handling.py +++ b/madmatrix/model_handling.py @@ -87,6 +87,9 @@ class MadMatrixALOHAWriter(aloha_writers.ALOHAWriterForGPU): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) + # a combined routine written as a wrapper needs a scratch wavefunction + # as an extra argument (see write_combined_cc) + self.combined_needs_tmp = True self.outname = 'w%s%s' % (self.particles[self.outgoing-1], self.outgoing) self.momentum_size = 0 # for ALOHAOBJ implementation the momentum is separated from the wavefunctions @@ -140,7 +143,10 @@ def get_header_txt(self, name=None, couplings=None,mode='', combined=False): """ if name is None: name = self.name - if fd_gauge and name.count("_") > 1: # FIXME ugly hack to get right header + if not self.combined_needs_tmp: + # assembled combined routine: it needs no scratch wavefunction + combined = False + elif fd_gauge and name.count("_") > 1: # FIXME ugly hack to get right header combined = True if mode=='': mode = self.mode @@ -229,12 +235,6 @@ def get_header_txt(self, name=None, couplings=None,mode='', combined=False): def get_foot_txt(self, combined=False): """Prototype for language specific footer""" text = ' ' - if not combined and fd_gauge: - if self.outgoing and 'P1N' not in self.tag: - name = self.particles[self.outgoing-1] - if name.startswith(('S','V')): - text += ' multiply_propagator_factor(%(name)s%(i)s,%(mass)s%(i)s, %(name)s%(i)s);\n' % \ - {'name':name, 'mass': 'M%s' % name[1:], 'i': self.outgoing } text +=' mgDebug( 1, __FUNCTION__ );\n' text +=' return;\n' text += ' }\n\n //--------------------------------------------------------------------------' # AV @@ -378,10 +378,50 @@ def get_momenta_txt(self): for i in range(1,6): cppindex = i -1 out.write(" w%(type)s%(out)s[%(ind)s] = CZERO ;\n" % {'type': type, 'out':self.outgoing, 'ind':cppindex}) + if self.has_fd_propagator(): + out.write(self.get_fd_gauge_txt()) # Returning result ###print('."' + out.getvalue() + '"') # AV - FOR DEBUGGING return out.getvalue() + # OM - FD gauge: the propagator factor is written in the routine (it used to + # be a call to multiply_propagator_factor), split in a part that only + # depends on the momentum and one that is linear in the wavefunction. + def get_fd_gauge_txt(self): + """FD gauge: the 5-momentum q (its 5th component carries the mass) and + the gauge direction n. Written with the momenta: it does not depend on + the wavefunction the routine builds.""" + + out = StringIO() + wf = '%s%s' % (self.particles[self.outgoing-1], self.outgoing) + # the 5th component is -i * m: built from broadcast reals, since a + # scalar complex has no conversion to the vector type (cxtype_sv) + out.write(' cxtype_sv FDQ[5] = { %s, cxmake( fptype_sv{ 0 }, -M%s + fptype_sv{ 0 } ) };\n' % + (', '.join('cxmake( -%s.pvec[%d], 0. )' % (wf, i) for i in range(4)), + self.outgoing)) + out.write(' fptype_sv FDN[5];\n') + out.write(' define_gauge_dir( FDQ, FDN );\n') + out.write(' const fptype_sv FDNQ = %s;\n' % + ' - '.join('FDN[%d] * FDQ[%d].real()' % (i, i) for i in range(4))) + return out.getvalue() + + def get_fd_propagator_txt(self): + """FD gauge: the part of the propagator factor that is linear in the + wavefunction, w -> w - q * js1 - n * js2, written after the components + have been built.""" + + out = StringIO() + w = self.outname + size = self.type_to_size[self.particles[self.outgoing-1]] - 2 + out.write(' const cxtype_sv FDJS1 = ( %s ) / FDNQ;\n' % + ' - '.join('FDN[%d] * %s[%d]' % (i, w, i) for i in range(4))) + out.write(' const cxtype_sv FDJS2 = ( %s - cxconj( FDQ[4] ) * %s[4] ) / FDNQ;\n' % + (' - '.join('FDQ[%d] * %s[%d]' % (i, w, i) for i in range(4)), w)) + for i in range(size): + out.write(' %(w)s[%(i)d] = %(w)s[%(i)d] - FDQ[%(i)d] * FDJS1 - FDN[%(i)d] * FDJS2;\n' + % {'w': w, 'i': i}) + return out.getvalue() + # AV - modify aloha_writers.ALOHAWriterForCPP method (improve formatting, add delayed declaration with initialisation) # This affects 'P1[0] = ' in HelAmps_sm.cc def get_one_momenta_def(self, i, strfile): @@ -565,12 +605,18 @@ def define_expression(self): format = ' %s = %s;\n' % (name, self.get_fct_format(fct)) out.write(format % ','.join([self.write_obj(obj) for obj in objs])) # AV not used in eemumu? numerator = self.routine.expr - if not 'Coup(1)' in self.routine.infostr: + if self.coup_name: + # one term of an assembled routine: its own coupling among COUP1, ... + coup_name = self.coup_name + elif not 'Coup(1)' in self.routine.infostr: coup_name = 'COUP' else: coup_name = '%s' % self.change_number_format(1) + has_coup = coup_name != self.change_number_format(1) + # OM the Ccoeff sign carrier goes with the coupling it multiplies + ccoeff = 'Ccoeff%s' % coup_name[4:] if not self.offshell: - if coup_name == 'COUP': + if has_coup: mydict = {'num': self.write_obj(numerator.get_rep([0]))} # '...(TMP4)-cI...' comes from here for c in ['coup', 'vertex']: if self.type2def['pointer_%s' %c] in ['*']: @@ -580,7 +626,10 @@ def define_expression(self): mydict['pre_%s' %c] = '' mydict['post_%s'%c] = '' # This affects '( *vertex ) = ' in HelAmps_sm.cc - out.write(' %(pre_vertex)svertex%(post_vertex)s = Ccoeff * %(pre_coup)sCOUP%(post_coup)s * %(num)s;\n' % mydict) # OM add Ccoeff (fix #825) + mydict['coup'] = coup_name + mydict['ccoeff'] = ccoeff + mydict['add'] = '%(pre_vertex)svertex%(post_vertex)s + ' % mydict if self.combined_part else '' + out.write(' %(pre_vertex)svertex%(post_vertex)s = %(add)s%(ccoeff)s * %(pre_coup)s%(coup)s%(post_coup)s * %(num)s;\n' % mydict) # OM add Ccoeff (fix #825) else: mydict= {} if self.type2def['pointer_vertex'] in ['*']: @@ -590,13 +639,16 @@ def define_expression(self): mydict['pre_vertex'] = '' mydict['post_vertex'] = '' mydict['data'] = self.write_obj(numerator.get_rep([0])) + mydict['add'] = '%(pre_vertex)svertex%(post_vertex)s + ' % mydict if self.combined_part else '' # This affects '( *vertex ) = ' in HelAmps_sm.cc - out.write(' %(pre_vertex)svertex%(post_vertex)s = %(data)s;\n' % mydict) + out.write(' %(pre_vertex)svertex%(post_vertex)s = %(add)s%(data)s;\n' % mydict) else: OffShellParticle = '%s%d' % (self.particles[self.offshell-1],\ self.offshell) if 'L' not in self.tag: - coeff = 'denom' + # each term of an assembled routine has its own denominator + # (they are 'const' declarations in the same scope) + coeff = 'denom%s' % (coup_name[4:] if self.combined_part else '') mydict = {} if self.type2def['pointer_coup'] in ['*']: mydict['pre_coup'] = '(*' @@ -607,13 +659,13 @@ def define_expression(self): mydict['coup'] = coup_name mydict['i'] = self.outgoing if self.nodeclare: - mydict['declnamedenom'] = 'const %s denom' % self.type2def['complex_v'] # AV + mydict['declnamedenom'] = 'const %s %s' % (self.type2def['complex_v'], coeff) # AV else: - mydict['declnamedenom'] = 'denom' # AV - self.declaration.add(('complex','denom')) + mydict['declnamedenom'] = coeff # AV + self.declaration.add(('complex', coeff)) # Need to add the unary operator before the coupling (OM fix for #825) - if mydict['coup'] != 'one': # but in case where the coupling is not used (one) - mydict['pre_coup'] = 'Ccoeff * %s' % mydict['pre_coup'] + if has_coup: # but in case where the coupling is not used (one) + mydict['pre_coup'] = '%s * %s' % (ccoeff, mydict['pre_coup']) if not aloha.complex_mass: # This affects 'denom = COUP' in HelAmps_sm.cc if self.routine.denominator: @@ -640,10 +692,22 @@ def define_expression(self): shift = 5 - 1 #to correspond to the shift in fortran indicies with -1 for C++ for ind in numerator.listindices(): # This affects 'V1[2] = ' and 'F1[2] = ' in HelAmps_sm.cc - ###out.write(' %s[%d]= %s*%s;\n' % (self.outname, - out.write(' %s[%d] = %s * %s;\n' % (self.outname, # AV - self.pass_to_HELAS(ind) + shift, coeff, + # the slot is fixed by the spin of this structure, the name by + # the routine the code is written in; a term of an assembled + # routine adds up into a slot several structures can feed + # outname carries the 'w' accessor prefix, rename_wf works on + # the wavefunction name itself + slot = 'w%s[%d]' % (self.rename_wf(self.outname[1:]), + self.pass_to_HELAS(ind) + shift) + out.write(' %s = %s%s * %s;\n' % (slot, # AV + '%s + ' % slot if self.combined_part else '', + coeff, self.write_obj(numerator.get_rep(ind)))) + if self.has_fd_propagator() and not self.combined_part: + # the FD gauge propagator factor is part of the wavefunction this + # routine builds, not a post-treatment (a term of an assembled + # routine gets it once, from the assembler, on the total) + out.write(self.get_fd_propagator_txt()) ###return out.getvalue() # AV # AV check if one, two, half or quarter are used and need to be defined (ugly hack for #291: can this be done better?) out2 = StringIO() @@ -695,7 +759,10 @@ def shift_indices(self, match): shift = -1 if fd_gauge and match.group('var').startswith('S'): shift += 4 - return 'w%s[%s]' % (match.group('var'), int(match.group('num')) + shift) + # the slot is fixed by the spin this structure sees on the leg, the + # name by the routine the code is written in (rename_wf) + return 'w%s[%s]' % (self.rename_wf(match.group('var')), + int(match.group('num')) + shift) # OM - overload aloha_writers.WriteALOHA and ALOHAWriterForCPP methods (handle 'unary minus' #628) def change_var_format(self, obj): @@ -781,6 +848,18 @@ def write_combined_cc(self, lor_names, offshell=None, sym=True, mode=''): if offshell is None: offshell = self.offshell + # structures acting on different spins (FD gauge: a leg is a vector in + # one structure and its Goldstone in the next one) can not share a + # single ALOHA expression, but they can still be assembled into a + # single routine instead of a wrapper calling each of them + self.combined_needs_tmp = True + if hasattr(self.routine, 'get_combined_routines') and \ + not os.environ.get('MG_ALOHA_COMBINE_WRAPPER'): + routines = self.routine.get_combined_routines(lor_names) + if routines: + self.combined_needs_tmp = False + return self.write_combined_parts_cc(routines, lor_names, offshell, mode) + name = combine_name(self.routine.name, lor_names, offshell, self.tag) self.name = name # write head - momenta - body - foot @@ -872,6 +951,78 @@ def write_combined_cc(self, lor_names, offshell=None, sym=True, mode=''): text = text.getvalue() return text + # OM - a combined routine whose structures do not act on the same spins is + # assembled term by term instead of calling each of them (see + # AbstractRoutine.get_combined_routines and the fortran writer) + def write_combined_parts_cc(self, routines, lor_names, offshell, mode=''): + """Assemble the structures of a combined call into a single routine: + one set of momenta, one set of declarations, one FD propagator factor, + and each structure adding its own coupling times its own expression + into the slots it feeds.""" + + name = combine_name(self.routine.name, lor_names, offshell, self.tag) + self.name = name + writers_l = [self.__class__(routine, None, options=self.options) + for routine in routines] + main = writers_l[0] + main.name = name + main.mode = mode if mode else self.mode + # assembled: no scratch wavefunction in the signature (the .h header is + # written from self, which write_combined_cc already flagged) + main.combined_needs_tmp = False + # one name per leg -- the one the first structure gives it + legs = ['%s%d' % (spin, i+1) for i, spin in enumerate(main.particles)] + + bodies = [] + for i, writer in enumerate(writers_l): + writer.name = name + writer.mode = main.mode + writer.coup_name = 'COUP%s' % (i+1) + writer.combined_part = True + for j, spin in enumerate(writer.particles): + wf = '%s%d' % (spin, j+1) + if wf != legs[j]: + writer.wf_rename[wf] = legs[j] + bodies.append(writer.define_expression()) + + new_couplings = ['COUP%s' % (i+1) for i in range(len(lor_names)+1)] + for writer in writers_l[1:]: + for entry in writer.declaration: + main.declaration.add(entry) + + text = StringIO() + text.write(main.get_header_txt(name=name, couplings=new_couplings, + mode=main.mode)) + # unlike the wrapper form, the assembled body needs cI + text.write(main.get_declaration_txt()) + text.write(main.get_momenta_txt()) + text.write(main.get_coupling_def()) + # the terms add up into the output, so it has to start from zero (in FD + # gauge the momenta already reset the wavefunction) + if not offshell: + text.write(' ( *vertex ) = cxzero_sv();\n') + elif not (fd_gauge and main.particles[main.outgoing-1] in ['S', 'V']): + for i in range(self.type_to_size[main.particles[main.outgoing-1]] - 2): + text.write(' %s[%d] = cxzero_sv();\n' % (main.outname, i)) + # the structures share the contraction cache, so the same TMP/FCT (and + # the same constexpr fptype) is built by several of them: a given name + # always stands for the same value, keep the first definition only + declared = re.compile(r'^(?:const|constexpr)\s+\S+\s+' + r'((?:TMP|FCT)\d+|one|two|half|quarter)\s*[=(]') + seen = set() + for body in bodies: + for line in body.splitlines(True): + found = declared.match(line.strip()) + if found: + if found.group(1) in seen: + continue + seen.add(found.group(1)) + text.write(line) + if main.has_fd_propagator(): + text.write(main.get_fd_propagator_txt()) + text.write(main.get_foot_txt()) + return text.getvalue() + from os.path import join as pjoin @@ -2822,8 +2973,6 @@ def generate_helas_call(self, argument): if isinstance(argument, helas_objects.HelasWavefunction): #arg['out'] = 'w_sv[%(out)d]' arg['out'] = 'aloha_obj[%(out)d]' - if fd_gauge and len(l) > 1: #FIXME: this is a hack to avoid a bug in the FD code - arg['out'] = arg['out'] + ", aloha_obj_tmp[0]" if aloha.complex_mass: arg['mass'] = 'm_pars->%(CM)s, ' else: @@ -2832,8 +2981,6 @@ def generate_helas_call(self, argument): #arg['out'] = '&_sv[%(out)d]' arg['out'] = '&_fp[%(out)d]' arg['out2'] = 'amp_sv[%(out)d]' - if fd_gauge and len(l) > 1: #FIXME: this is a hack to avoid a bug in the FD code - arg['out'] = arg['out'] + ", &_tmp_fp[0]" arg['mass'] = '' call = call % arg # Now we have a line correctly formatted diff --git a/tests/acceptance_tests/test_cmd.py b/tests/acceptance_tests/test_cmd.py index 0f9a64b0b..5ea141c5c 100755 --- a/tests/acceptance_tests/test_cmd.py +++ b/tests/acceptance_tests/test_cmd.py @@ -27,6 +27,7 @@ logger = logging.getLogger('test_cmd') import tests.unit_tests.iolibs.test_file_writers as test_file_writers +import tests.IOTests as IOTests import madgraph.interface.master_interface as Cmd import madgraph.interface.launch_ext_program as launch_ext @@ -1534,6 +1535,162 @@ def get_values(output_format): # gauge invariance: unitary-gauge values must also match the FD ones. self._assert_me_lists_close(standalone_mg7_no_fd, standalone, atol=1e-7) + def test_standalone_mg7_fd_vs_fortran(self): + """FD gauge: madmatrix and the Fortran standalone must agree on the + value, not merely both produce one. + + test_standalone_cpp_fd_output_consistency also compares the two + backends in FD gauge, but its process (h + 4 jets) has matrix elements + of order 1e-15 to 1e-21 while it asserts with atol=1e-7, so no wrong + value can make it fail. Two madmatrix FD bugs lived through it: the + Goldstone component of a longitudinal massive vector had the wrong sign + on the scalar backends (0.7% off), and on the SIMD ones the FD block of + vxxxxx broadcast its constants to the first lane only (nan elsewhere). + + u u~ > w+ w- is the smallest process that pins both down: the W's give + a longitudinal polarisation (hence the Goldstone component), the + massless initial state makes the two drivers evaluate the same RAMBO + point, and |M|^2 is O(1e-3), far above the noise floor -- so the + comparison is done with a relative tolerance and no absolute floor. + The magnitude is asserted too, to keep this test from silently decaying + into a vacuous one. + + madmatrix is checked once per backend: the FD wavefunctions are written + twice in helas_fd.h, once under #ifndef MGONGPU_CPPSIMD and once for the + vector types, and the two have drifted apart before. Which one a plain + 'make' builds depends on the host (cppauto), so neither is exercised + unless it is asked for by name. + """ + energy = '1000' + devnull = open(os.devnull, 'w') + me_re = re.compile(r'Matrix element\s*=\s*([\d.eE+-]+)\s*GeV', + re.IGNORECASE) + + def get_values(output_format, check_exe, build_source=False, + backend=None): + if os.path.isdir(self.out_dir): + shutil.rmtree(self.out_dir) + self.do('output %s %s -f' % (output_format, self.out_dir)) + if build_source: + subprocess.call(['make'], stdout=devnull, stderr=devnull, + cwd=pjoin(self.out_dir, 'Source')) + proc_root = pjoin(self.out_dir, 'SubProcesses') + dirs = sorted(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' % output_format) + values = [] + for d in dirs: + proc_dir = pjoin(proc_root, d) + target = ['make', 'check'] if output_format == 'standalone' \ + else ['make'] + (['BACKEND=%s' % backend] if backend else []) + subprocess.call(target, stdout=devnull, stderr=devnull, + cwd=proc_dir) + log = pjoin(proc_dir, 'check.log') + subprocess.call('%s %s' % (check_exe, energy), + stdout=open(log, 'w'), stderr=subprocess.STDOUT, + cwd=proc_dir, shell=True) + found = me_re.findall(open(log).read()) + self.assertTrue(found, '%s produced no matrix element (see %s)' + % (output_format, log)) + values.extend(float(v) for v in found) + return values + + self.do('import model sm') + self.do('set gauge FD') + try: + self.do('generate u u~ > w+ w-') + standalone = get_values('standalone', './check', build_source=True) + # cppnone is the scalar code path, cppsse4 the vector one (it maps + # to NEON on arm) + mg7 = dict((backend, + get_values('standalone_mg7', './check_sa.exe', + backend=backend)) + for backend in ('cppnone', 'cppsse4')) + finally: + self.do('set gauge unitary') + + # the comparison is only meaningful if the values are not noise + self.assertGreater(max(abs(v) for v in standalone), 1e-6, + 'matrix elements too small for a relative ' + 'comparison to mean anything: %s' % standalone) + for backend, values in sorted(mg7.items()): + # -ffast-math on the C++ side puts the backends ~3e-8 apart + self._assert_me_lists_close(values, standalone, rtol=1e-6) + + def test_standalone_mg7_fd_simd_lanes(self): + """FD gauge: the scalar and the vectorised madmatrix backends must + compute the same thing, over many events. + + The FD wavefunctions are written twice in helas_fd.h -- once under + #ifndef MGONGPU_CPPSIMD, once for the vector types -- and only the first + event of the first lane is ever printed by check_sa.exe 'matrix' mode. + A vector branch that is right in the first lane and wrong in the others + is therefore invisible to a single point: the FD block of vxxxxx built + its constants with a braced initialiser, which fills the first lane and + zeroes the rest, and that gave nan (and later, once the surrounding code + used selects, merely wrong numbers) in every other lane. + + 'perf' mode runs many events through every lane, so comparing the mean + matrix element of a scalar and a vectorised build covers them all. The + two agree to the printed precision when the branches agree, and differed + by 8e-5 relative with that bug in place. + """ + devnull = open(os.devnull, 'w') + if os.path.isdir(self.out_dir): + shutil.rmtree(self.out_dir) + + self.do('import model sm') + self.do('set gauge FD') + try: + self.do('generate u u~ > w+ w-') + self.do('output standalone_mg7 %s -f' % self.out_dir) + finally: + self.do('set gauge unitary') + + proc_root = pjoin(self.out_dir, 'SubProcesses') + dirs = sorted(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 standalone_mg7') + + def mean_me(proc_dir, backend): + """mean |M|^2 over a multi-event run of the given backend""" + exe = pjoin(proc_dir, 'check_sa.exe') + if os.path.exists(exe): + os.remove(exe) + # USEBUILDDIR keeps the two backends' objects apart + subprocess.call(['make', 'BACKEND=%s' % backend, 'USEBUILDDIR=1'], + stdout=devnull, stderr=devnull, cwd=proc_dir) + if not os.path.exists(exe): + return None # backend not available on this host + log = pjoin(proc_dir, 'perf_%s.log' % backend) + subprocess.call('./check_sa.exe perf 1 32 300', + stdout=open(log, 'w'), stderr=subprocess.STDOUT, + cwd=proc_dir, shell=True) + out = open(log).read() + found = re.search(r'MeanMatrixElemValue\s*=\s*\(\s*([^\s+]+)', out) + self.assertTrue(found, 'no mean matrix element (see %s)' % log) + value = found.group(1) + self.assertNotIn('nan', value.lower(), + 'nan matrix element with BACKEND=%s (see %s)' + % (backend, log)) + return float(value) + + for d in dirs: + proc_dir = pjoin(proc_root, d) + # cppnone is the scalar code path, cppsse4 the vector one (it maps + # to NEON on arm) + scalar = mean_me(proc_dir, 'cppnone') + self.assertTrue(scalar, 'standalone_mg7 did not build in %s' % proc_dir) + self.assertGreater(scalar, 0., + 'null mean matrix element in %s' % proc_dir) + vector = mean_me(proc_dir, 'cppsse4') + if vector is None: + continue # no vectorised backend here: nothing to compare + self.assertLessEqual(abs(vector - scalar), 1e-5 * scalar, + 'the scalar and vectorised backends disagree ' + 'in %s: %s vs %s' % (d, scalar, vector)) + def test_standalone_mg7_vs_cpp(self): """Cross-check that standalone_mg7 (madmatrix) reproduces the standalone_cpp matrix elements for p p > e+ e- QCD=0. @@ -3640,7 +3797,8 @@ def check_aloha_file(self): C This File is Automatically generated by ALOHA C The process calculated in this file is: -C Gamma(3,2,-1)*ProjM(-1,1) +C Coup(1) * (Gamma(3,2,-1)*ProjM(-1,1)) + Coup(2) * +C (Gamma(3,2,-1)*ProjM(-1,1) + 2*Gamma(3,2,-1)*ProjP(-1,1)) C SUBROUTINE FFV2_4_3(F1, F2, COUP1, COUP2, M3, W3,V3) USE ALOHA_OBJECT @@ -3656,18 +3814,48 @@ def check_aloha_file(self): REAL*8 M3 REAL*8 OM3 REAL*8 P3(0:3) + COMPLEX*16 TMP2 + COMPLEX*16 TMP5 TYPE(ALOHA) V3 - TYPE(ALOHA) VTMP REAL*8 W3 COMPLEX*16 DENOM - INTEGER*4 I - CALL FFV2_3(F1,F2,COUP1,M3,W3,V3) - CALL FFV4_3(F1,F2,COUP2,M3,W3,VTMP) - DO I = 1, 4 - V3 %W(I) = V3%W(I) + VTMP%W(I) - ENDDO + OM3 = 0D0 + IF (M3.NE.0D0) OM3=1D0/M3**2 + V3%P(:) = +F1%P(:)+F2%P(:) + P3(:) = -V3 % P (:) + FLV_INDEX1 = F1 %FLV_INDEX + FLV_INDEX2 = F2 %FLV_INDEX + IF(FLV_INDEX1.NE.FLV_INDEX2.OR.FLV_INDEX1.EQ.0)THEN + V3%W(:) = (0D0,0D0) + RETURN + ENDIF + TMP2 = (F1 % W(1)*(F2 % W(3)*(P3(0)+P3(3))+F2 % W(4)*(P3(1)+CI + $ *(P3(2))))+F1 % W(2)*(F2 % W(3)*(P3(1)-CI*(P3(2)))+F2 % W(4) + $ *(P3(0)-P3(3)))) + TMP5 = (F1 % W(3)*(F2 % W(1)*(P3(0)-P3(3))-F2 % W(2)*(P3(1)+CI + $ *(P3(2))))+F1 % W(4)*(F2 % W(1)*(-P3(1)+CI*(P3(2)))+F2 % W(2) + $ *(P3(0)+P3(3)))) + DENOM = 1D0/(P3(0)**2-P3(1)**2-P3(2)**2-P3(3)**2 - M3 * (M3 -CI* + $ W3)) + V3%W(1)= DENOM*(-2D0 * CI)*(COUP2*(OM3*-1D0/2D0 * P3(0)*(TMP2 + $ +2D0*(TMP5))+(+1D0/2D0*(F1 % W(1)*F2 % W(3)+F1 % W(2)*F2 % W(4)) + $ +F1 % W(3)*F2 % W(1)+F1 % W(4)*F2 % W(2)))+1D0/2D0*(COUP1*(F1 % + $ W(1)*F2 % W(3)+F1 % W(2)*F2 % W(4)-P3(0)*OM3*TMP2))) + V3%W(2)= DENOM*(-2D0 * CI)*(COUP2*(OM3*-1D0/2D0 * P3(1)*(TMP2 + $ +2D0*(TMP5))+(-1D0/2D0*(F1 % W(1)*F2 % W(4)+F1 % W(2)*F2 % W(3)) + $ +F1 % W(3)*F2 % W(2)+F1 % W(4)*F2 % W(1)))-1D0/2D0*(COUP1*(F1 % + $ W(1)*F2 % W(4)+F1 % W(2)*F2 % W(3)+P3(1)*OM3*TMP2))) + V3%W(3)= DENOM*CI*(COUP2*(OM3*P3(2)*(TMP2+2D0*(TMP5))+(+CI*(F1 % + $ W(1)*F2 % W(4))-CI*(F1 % W(2)*F2 % W(3))-2D0 * CI*(F1 % W(3) + $ *F2 % W(2))+2D0 * CI*(F1 % W(4)*F2 % W(1))))+COUP1*(+CI*(F1 % + $ W(1)*F2 % W(4))-CI*(F1 % W(2)*F2 % W(3))+P3(2)*OM3*TMP2)) + V3%W(4)= DENOM*2D0 * CI*(COUP2*(OM3*1D0/2D0 * P3(3)*(TMP2+2D0 + $ *(TMP5))+(+1D0/2D0*(F1 % W(1)*F2 % W(3))-1D0/2D0*(F1 % W(2)*F2 + $ % W(4))-F1 % W(3)*F2 % W(1)+F1 % W(4)*F2 % W(2)))+1D0/2D0 + $ *(COUP1*(F1 % W(1)*F2 % W(3)+P3(3)*OM3*TMP2-F1 % W(2)*F2 % W(4)) + $ )) END - + """ text = open(os.path.join(self.out_dir,'Source', 'DHELAS', 'FFV2_3.f')).read() @@ -4656,3 +4844,73 @@ def test_import_banner_command(self): self.assertIn("200 = nevents", run_card) os.chdir(cwd) + + +#=============================================================================== +# IOTestFDGauge +#=============================================================================== +class IOTestFDGauge(IOTests.IOTestManager): + """Reference files for the FD gauge output, on both backends. + + FD gauge writes wavefunctions and routines that no other gauge does: a + massive vector and its Goldstone share one 5 component object, the + propagator factor is part of the routine, and the Lorentz structures of a + vertex are assembled even though they act on different spins. None of that + was covered by a reference file, and three madmatrix bugs (a Goldstone sign, + a vector constant broadcast to the first lane only, and a gauge direction + that differed between the scalar and vector branches) lived in the released + templates because of it. Two of them show up as a wrong number, which + test_standalone_mg7_fd_vs_fortran and test_standalone_mg7_fd_simd_lanes now + catch; the third is a gauge choice, invisible to any matrix element, and + only a reference file can hold it still. + + u u~ > w+ w- is small and covers what is specific to FD: the merged flavour + routines, a same spin combination (FFV6_2M_3), two mixed spin ones + (VVV1_VVS1_VSV2_VSS1_0 and VVV1_VSV2_VSS2_SVV2_SVS2_SSV3_0) and the inlined + propagator factor of every offshell V/S routine. + """ + + def setUp(self): + super(IOTestFDGauge, self).setUp() + self.interface = Cmd.MasterCmd() + self.interface.no_notification() + + def generate_fd(self, output_format, path): + """u u~ > w+ w- in FD gauge, in the given output format""" + cmds = ['import model sm', + 'set gauge FD', + 'generate u u~ > w+ w-', + 'output %s %s -f' % (output_format, path)] + try: + for cmd in cmds: + self.interface.exec_cmd(cmd, errorhandling=False, printcmd=False, + precmd=True, postcmd=True) + finally: + # the gauge is global: do not leak it into the next test + self.interface.exec_cmd('set gauge unitary', errorhandling=False, + printcmd=False, precmd=True, postcmd=True) + + @IOTests.createIOTest() + def testIO_FDgauge_standalone_fortran(self): + r""" target: FD_fortran/Source/DHELAS/FFV6M_3.f + target: FD_fortran/Source/DHELAS/VVV1_0.f + """ + # Two files rather than the whole DHELAS, chosen to cover what FD gauge + # does that no other gauge does: + # - FFV6M_3.f: the propagator factor written into the body of an + # offshell routine (q, the gauge direction, js1/js2 and the update of + # the 5 components), and a combination merged into one expression + # (FFV6_2M_3); + # - VVV1_0.f: the combinations assembled out of structures that act on + # different spins (VVV1_VVS1_VSV2_VSS1_0 and + # VVV1_VSV2_VSS2_SVV2_SVS2_SSV3_0), which only exist because a + # massive vector and its Goldstone are the same wavefunction here. + self.generate_fd('standalone', pjoin(self.IOpath, 'FD_fortran')) + + @IOTests.createIOTest() + def testIO_FDgauge_madmatrix(self): + r""" target: FD_madmatrix/src/HelAmps_sm.h + """ + # one file holds all of it for madmatrix: the generated routines and, + # pasted above them, the FD helpers of helas_fd.h + self.generate_fd('standalone_mg7', pjoin(self.IOpath, 'FD_madmatrix')) diff --git a/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_madmatrix/%FD_madmatrix%src%HelAmps_sm.h b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_madmatrix/%FD_madmatrix%src%HelAmps_sm.h new file mode 100644 index 000000000..4aaeaa504 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_madmatrix/%FD_madmatrix%src%HelAmps_sm.h @@ -0,0 +1,2112 @@ +// Copyright (C) 2010 The ALOHA Development team and Contributors. +// Copyright (C) 2010 The MadGraph5_aMC@NLO development team and contributors. +// Created by: J. Alwall (Sep 2010) for the MG5aMC backend. +//========================================================================== +// Copyright (C) 2020-2026 CERN and UCLouvain. +// Licensed under the GNU Lesser General Public License (version 3 or later). +// Modified originally by: A. Valassi (Sep 2021) for the MG5aMC CUDACPP plugin. +// Further modified by: J. Teig, A. Valassi (2021-2024) for the MG5aMC CUDACPP plugin. +// Integrated with the MadGraph7 project in Feb 2026. +//========================================================================== +// This file has been automatically generated for CUDA/C++ standalone by +// MadGraph5_aMC@NLO v. %(version)s, %(date)s +// By the MadGraph5_aMC@NLO Development Team +// Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +//========================================================================== + +#ifndef HelAmps_sm_H +#define HelAmps_sm_H 1 + +#include "mgOnGpuConfig.h" + +#include "mgOnGpuVectors.h" + +#include "Parameters.h" + +#include +//#include +//#include +//#include +//#include + +#ifdef MGONGPUCPP_GPUIMPL +namespace mg5amcGpu +#else +namespace mg5amcCpu +#endif +{ + + // ALOHA-style object for easy flavor consolidation and non-template API + struct ALOHAOBJ { + + static constexpr int np4 = 4; // dimensions of 4-momenta (E,px,py,pz) + static constexpr int nw6 = 5; // dimensions of each wavefunction (notice, this is +1 in case of FD gauge) + fptype_sv * pvec; + fptype * w; + int flv_index; + + __host__ __device__ ALOHAOBJ() {} + __host__ __device__ ALOHAOBJ(fptype_sv * pvec_sv, cxtype_sv * w_sv, int flv = -1) + : pvec(pvec_sv), w(reinterpret_cast(w_sv)), flv_index(flv) {} + }; + + struct FLV_COUPLING_VIEW { + + const int* const partner1; + const int* const partner2; + const fptype* const value; + + __host__ __device__ + FLV_COUPLING_VIEW(const int* p1, const int* p2, const fptype* v) + : partner1(p1), partner2(p2), value(v) {} + }; + + // FSTRIDE is the number of fptype's used to store one flavor slot of the value buffer: + // - independent (fixed) flavored couplings: FSTRIDE = nx2 = 2 (a single scalar complex, broadcast across the SIMD vector) + // - dependent (event-by-event, running-alphas) flavored couplings: FSTRIDE = nx2*neppC (an AOSOA SIMD record) + // It must match C_ACCESS::flv_stride of the access type the consuming vertex routine is instantiated with. + template + class FLV_COUPLING_ARRAY { + + static_assert(SIZE >= 0, "flvCOUPs SIZE must be non-negative"); + static_assert(STRIDE > 0, "flvCOUPs STRIDE must be positive"); + static_assert(FSTRIDE > 0, "flvCOUPs FSTRIDE must be positive"); + const int* const partner1; + const int* const partner2; + const fptype* const value; + + public: + __host__ __device__ + FLV_COUPLING_ARRAY(const int* p1, const int* p2, const fptype* v) + : partner1(p1), partner2(p2), value(v) {} + + __host__ __device__ + FLV_COUPLING_VIEW operator[](const int i) const { + return FLV_COUPLING_VIEW{ + partner1 + i*STRIDE, + partner2 + i*STRIDE, + value + i*FSTRIDE*STRIDE + }; + } + }; + //-------------------------------------------------------------------------- + +#ifdef MGONGPU_INLINE_HELAMPS +#define INLINE inline +#define ALWAYS_INLINE __attribute__( ( always_inline ) ) +#else +#define INLINE +#define ALWAYS_INLINE +#endif + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ INLINE void + ixxxxx( const fptype momenta[], // input: momenta + const fptype fmass, // input: fermion mass + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fi, // output: aloha objects + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0) + template + __host__ __device__ INLINE void + ipzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fi, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0) + template + __host__ __device__ INLINE void + imzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fi, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PT > 0) + template + __host__ __device__ INLINE void + ixzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fi, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction vc[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ INLINE void + vxxxxx( const fptype momenta[], // input: momenta + const fptype vmass, // input: vector boson mass + const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const int nsv, // input: +1 (final) or -1 (initial) + const int flv, // input: flavor index + ALOHAOBJ & vc, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction sc[3] from the input momenta[npar*4*nevt] + template + __host__ __device__ INLINE void + sxxxxx( const fptype momenta[], // input: momenta + //const fptype, // WARNING: input "smass" unused (missing in Fortran) - scalar boson mass + //const int, // WARNING: input "nhel" unused (missing in Fortran) - scalar has no helicity! + const int nss, // input: +1 (final) or -1 (initial) + const int flv, // input: flavor index + ALOHAOBJ & sc, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ INLINE void + oxxxxx( const fptype momenta[], // input: momenta + const fptype fmass, // input: fermion mass + const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fo, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0) + template + __host__ __device__ INLINE void + opzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fo, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0) + template + __host__ __device__ INLINE void + omzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fo, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ INLINE void + oxzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavor index + ALOHAOBJ & fo, // output: wavefunctions + const int ipar // input: particle# out of npar + ) ALWAYS_INLINE; + + +//-------------------------------------------------------------------------- + + // Compute the direction n[5] of the gauge q[5] + __host__ __device__ INLINE void + define_gauge_dir( const fptype q[], // input: gauge + fptype n[] // output: direction + ) ALWAYS_INLINE; + + + //-------------------------------------------------------------------------- + // Compute a propagator factor d out of gauge q[5] and a mass + __host__ __device__ INLINE void + calculate_propagator_factor( const fptype_sv q[5], // input: gauge + const fptype_sv mass, // input: mass + fptype_sv *d // output: propagator factor + ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + // multiply by propagation factor from m and wawefunctionsin[] and output them + // as wavefunctionout[] + template< class W_ACCESS> + __host__ __device__ INLINE void + multiply_propagator_factor( const fptype wavefunctionsin[], // input: wavefunctions + const fptype m, // input: mass + fptype wavefunctionsout[] // output: wavefunctions + ) ALWAYS_INLINE; +//========================================================================== + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ void + ixxxxx( const fptype momenta[], // input: momenta + const fptype fmass, // input: fermion mass + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fi, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + // NEW IMPLEMENTATION FIXING FLOATING POINT EXCEPTIONS IN SIMD CODE (#701) + // Variables xxxDENOM are a hack to avoid division-by-0 FPE while preserving speed (#701 and #727) + // Variables xxxDENOM are declared as 'volatile' to make sure they are not optimized away on clang! (#724) + // A few additional variables are declared as 'volatile' to avoid sqrt-of-negative-number FPEs (#736) + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fi.w ); + fi.pvec[0] = -pvec0 * (fptype)nsf; + fi.pvec[1] = -pvec1 * (fptype)nsf; + fi.pvec[2] = -pvec2 * (fptype)nsf; + fi.pvec[3] = -pvec3 * (fptype)nsf; + fi.flv_index = flv; + const int nh = nhel * nsf; + if( fmass != 0. ) + { +#ifndef MGONGPU_CPPSIMD + const fptype_sv pp = fpmin( pvec0, fpsqrt( pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3 ) ); +#else + volatile fptype_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 + const fptype_sv pp = fpmin( pvec0, fpsqrt( p2 ) ); +#endif + // In C++ ixxxxx, use a single ip/im numbering that is valid both for pp==0 and pp>0, which have two numbering schemes in Fortran ixxxxx: + // for pp==0, Fortran sqm(0:1) has indexes 0,1 as in C++; but for Fortran pp>0, omega(2) has indexes 1,2 and not 0,1 + // NB: this is only possible in ixxxx, but in oxxxxx two different numbering schemes must be used + const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 + const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 +#ifndef MGONGPU_CPPSIMD + if( pp == 0. ) + { + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses + //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? + sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here + w[0] = cxmake( ip * sqm[ip], 0 ); + w[1] = cxmake( im * nsf * sqm[ip], 0 ); + w[2] = cxmake( ip * nsf * sqm[im], 0 ); + w[3] = cxmake( im * sqm[im], 0 ); + } + else + { + const fptype sf[2] = { fptype( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype)0.5, + fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; + fptype omega[2] = { fpsqrt( pvec0 + pp ), 0. }; + omega[1] = fmass / omega[0]; + const fptype sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype pp3 = fpmax( pp + pvec3, 0. ); + const cxtype chi[2] = { cxmake( fpsqrt( pp3 * (fptype)0.5 / pp ), 0. ), + ( pp3 == 0. ? cxmake( -nh, 0. ) : cxmake( nh * pvec1, pvec2 ) / fpsqrt( 2. * pp * pp3 ) ) }; + w[0] = sfomega[0] * chi[im]; + w[1] = sfomega[0] * chi[ip]; + w[2] = sfomega[1] * chi[im]; + w[3] = sfomega[1] * chi[ip]; + } +#else + // Branch A: pp == 0. + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses (NB: SCALAR!) + sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) + const cxtype fiA_2 = ip * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 + const cxtype fiA_3 = im * nsf * sqm[ip]; // scalar cxtype: real part initialised from fptype, imag part = 0 + const cxtype fiA_4 = ip * nsf * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 + const cxtype fiA_5 = im * sqm[im]; // scalar cxtype: real part initialised from fptype, imag part = 0 + // Branch B: pp != 0. + const fptype sf[2] = { fptype( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype)0.5, + fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; + fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; + omega[1] = fmass / omega[0]; + const fptype_v sfomega[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype_v pp3 = fpmax( pp + pvec3, 0 ); + volatile fptype_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 + volatile fptype_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 + volatile fptype_v chi0r2 = pp3 * 0.5 / ppDENOM; // volatile fixes #736 + const cxtype_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0 ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 + cxternary( ( pp3 == 0. ), + cxmake( -nh, 0 ), + cxmake( (fptype)nh * pvec1, pvec2 ) / fpsqrt( 2. * ppDENOM * pp3DENOM ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_v fiB_2 = sfomega[0] * chi[im]; + const cxtype_v fiB_3 = sfomega[0] * chi[ip]; + const cxtype_v fiB_4 = sfomega[1] * chi[im]; + const cxtype_v fiB_5 = sfomega[1] * chi[ip]; + // Choose between the results from branch A and branch B + const bool_v mask = ( pp == 0. ); + w[0] = cxternary( mask, fiA_2, fiB_2 ); + w[1] = cxternary( mask, fiA_3, fiB_3 ); + w[2] = cxternary( mask, fiA_4, fiB_4 ); + w[3] = cxternary( mask, fiA_5, fiB_5 ); +#endif + } + else + { +#ifdef MGONGPU_CPPSIMD + volatile fptype_sv p0p3 = fpmax( pvec0 + pvec3, 0 ); // volatile fixes #736 + volatile fptype_sv sqp0p3 = fpternary( ( pvec1 == 0. and pvec2 == 0. and pvec3 < 0. ), + fptype_sv{ 0 }, + fpsqrt( p0p3 ) * (fptype)nsf ); + volatile fptype_sv sqp0p3DENOM = fpternary( sqp0p3 != 0, (fptype_sv)sqp0p3, 1. ); // hack: dummy sqp0p3DENOM[ieppV]=1 if sqp0p3[ieppV]==0 + cxtype_sv chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), + cxternary( sqp0p3 == 0, + cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ), + cxmake( (fptype)nh * pvec1, pvec2 ) / (const fptype_v)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 +#else + const fptype_sv sqp0p3 = fpternary( ( pvec1 == 0. and pvec2 == 0. and pvec3 < 0. ), + fptype_sv{ 0 }, + fpsqrt( fpmax( pvec0 + pvec3, 0. ) ) * (fptype)nsf ); + const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), + ( sqp0p3 == 0. ? cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ) : cxmake( (fptype)nh * pvec1, pvec2 ) / sqp0p3 ) }; +#endif + if( nh == 1 ) + { + w[0] = cxzero_sv(); + w[1] = cxzero_sv(); + w[2] = chi[0]; + w[3] = chi[1]; + } + else + { + w[0] = chi[1]; + w[1] = chi[0]; + w[2] = cxzero_sv(); + w[3] = cxzero_sv(); + } + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0) + template + __host__ __device__ void + ipzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fi, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fi.w ); + fi.pvec[0] = -pvec3 * (fptype)nsf; + fi.pvec[1] = fptype_sv{ 0 }; + fi.pvec[2] = fptype_sv{ 0 }; + fi.pvec[3] = -pvec3 * (fptype)nsf; + fi.flv_index = flv; + const int nh = nhel * nsf; + const cxtype_sv sqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); + w[0] = cxmake( fi.pvec[1], fi.pvec[2] ); + if( nh == 1 ) + { + w[1] = cxmake( fi.pvec[1], fi.pvec[2] ); + w[2] = sqp0p3; + } + else + { + w[1] = sqp0p3; + w[2] = cxmake( fi.pvec[1], fi.pvec[2] ); + } + w[3] = cxmake( fi.pvec[1], fi.pvec[2] ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0) + template + __host__ __device__ void + imzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fi, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fi.w ); + fi.pvec[0] = pvec3 * (fptype)nsf; + fi.pvec[1] = fptype_sv{ 0 }; + fi.pvec[2] = fptype_sv{ 0 }; + fi.pvec[3] = -pvec3 * (fptype)nsf; + fi.flv_index = flv; + const int nh = nhel * nsf; + const cxtype_sv chi = cxmake( -(fptype)nhel * fpsqrt( -2. * pvec3 ), 0. ); + w[1] = cxzero_sv(); + w[2] = cxzero_sv(); + if( nh == 1 ) + { + w[0] = cxzero_sv(); + w[3] = chi; + } + else + { + w[0] = chi; + w[3] = cxzero_sv(); + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PT > 0) + template + __host__ __device__ void + ixzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fi, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fi.w ); + fi.pvec[0] = -pvec0 * (fptype)nsf; + fi.pvec[1] = -pvec1 * (fptype)nsf; + fi.pvec[2] = -pvec2 * (fptype)nsf; + fi.pvec[3] = -pvec3 * (fptype)nsf; + fi.flv_index = flv; + const int nh = nhel * nsf; + //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? + const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; + const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); + const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); + if( nh == 1 ) + { + w[0] = cxzero_sv(); + w[1] = cxzero_sv(); + w[2] = chi0; + w[3] = chi1; + } + else + { + w[0] = chi1; + w[1] = chi0; + w[2] = cxzero_sv(); + w[3] = cxzero_sv(); + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction vc[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ void + vxxxxx( const fptype momenta[], // input: momenta + const fptype vmass, // input: vector boson mass + const int nhel, // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson) + const int nsv, // input: +1 (final) or -1 (initial) + const int flv, // input: flavour + ALOHAOBJ & vc, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + // NEW IMPLEMENTATION FIXING FLOATING POINT EXCEPTIONS IN SIMD CODE (#701) + // Variables xxxDENOM are a hack to avoid division-by-0 FPE while preserving speed (#701 and #727) + // Variables xxxDENOM are declared as 'volatile' to make sure they are not optimized away on clang! (#724) + // A few additional variables are declared as 'volatile' to avoid sqrt-of-negative-number FPEs (#736) + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( vc.w ); + vc.pvec[0] = pvec0 * (fptype)nsv; + vc.pvec[1] = pvec1 * (fptype)nsv; + vc.pvec[2] = pvec2 * (fptype)nsv; + vc.pvec[3] = pvec3 * (fptype)nsv; + vc.flv_index = flv; + const fptype sqh = fpsqrt( 0.5 ); // AV this is > 0! + const fptype hel = nhel; + + // FD gauge + const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1 + fptype_sv{ 0 } ); + fptype_sv n[5]; + fptype_sv nk; + // NB: broadcast to every SIMD lane. 'fptype_sv one{1.}' sets the first + // lane only (the others are zero filled), which left n, and then nk, at + // zero in all lanes but the first one + const fptype_sv zero = 0. + fptype_sv{ 0 }; + const fptype_sv one = 1. + fptype_sv{ 0 }; + + if( vmass != 0. ) + { + const int nsvahl = nsv * std::abs( hel ); + const fptype hel0 = 1. - std::abs( hel ); +#ifndef MGONGPU_CPPSIMD + const fptype_sv pt2 = ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ); + const fptype_sv pp = fpmin( pvec0, fpsqrt( pt2 + ( pvec3 * pvec3 ) ) ); + const fptype_sv pt = fpmin( pp, fpsqrt( pt2 ) ); + if( pp == 0. ) + { + w[0] = cxmake( 0., 0. ); + w[1] = cxmake( -hel * sqh, 0. ); + w[2] = cxmake( 0., nsvahl * sqh ); + w[3] = cxmake( hel0, 0. ); + } + else + { + //printf( "DEBUG1011 (before emp): pvec0=%f vmass=%f pp=%f vmass*pp=%f\n", pvec0, vmass, pp, vmass * pp ); + //const fptype emp = pvec / ( vmass * pp ); // this may give a FPE #1011 (why?! maybe when vmass=+-epsilon?) + const fptype emp = pvec0 / vmass / pp; // workaround for FPE #1011 + //printf( "DEBUG1011 (after emp): emp=%f\n", emp ); + w[0] = cxmake( hel0 * pp / vmass, 0. ); + w[3] = cxmake( hel0 * pvec3 * emp + hel * pt / pp * sqh, 0. ); + if( pt != 0. ) + { + const fptype pzpt = pvec3 / ( pp * pt ) * sqh * hel; + w[1] = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -nsvahl * pvec2 / pt * sqh ); + w[2] = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, nsvahl * pvec1 / pt * sqh ); + } + else + { + w[1] = cxmake( -hel * sqh, 0. ); + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + //vc[4] = cxmake( 0., nsvahl * ( pvec3 < 0. ? -std::abs( sqh ) : std::abs( sqh ) ) ); // AV: why abs here? + w[2] = cxmake( 0., nsvahl * ( pvec3 < 0. ? -sqh : sqh ) ); // AV: removed an abs here + } + } + + //FD gauge + if( pp > 0. ) + { + n[0] = ( pvec0 >= zero) ? one : -one; + n[1] = -pvec1/pp; + n[2] = -pvec2/pp; + n[3] = -pvec3/pp; + n[4] = zero; + } + else + { + n[0] = ( pvec0 >= zero) ? one : -one; + n[1] = zero; + n[2] = zero; + n[3] = ( pvec0 >= zero) ? -one : one; + } + + + nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; + + if ( abs(nhel) == 1) + { + w[4] = cxzero_sv(); + } + else{ + w[0] = cxmake( -vmass/nk * n[0], zero ); + w[1] = cxmake( -vmass/nk * n[1], zero ); + w[2] = cxmake( -vmass/nk * n[2], zero ); + w[3] = cxmake( -vmass/nk * n[3], zero ); + w[4] = -static_cast(nsv)*cI; // as in fortran vxxxxx (vc%W(5) = -nsv*ci) and in the SIMD branch below + } + +#else + + volatile fptype_sv pt2 = ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ); + volatile fptype_sv p2 = pt2 + ( pvec3 * pvec3 ); // volatile fixes #736 + const fptype_sv pp = fpmin( pvec0, fpsqrt( p2 ) ); + const fptype_sv pt = fpmin( pp, fpsqrt( pt2 ) ); + // Branch A: pp == 0. + const cxtype vcA_2 = cxmake( 0, 0 ); + const cxtype vcA_3 = cxmake( -hel * sqh, 0 ); + const cxtype vcA_4 = cxmake( 0, nsvahl * sqh ); + const cxtype vcA_5 = cxmake( hel0, 0 ); + // Branch B: pp != 0. + volatile fptype_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 + const fptype_v emp = pvec0 / ( vmass * ppDENOM ); // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_v vcB_2 = cxmake( hel0 * pp / vmass, 0 ); + const cxtype_v vcB_5 = cxmake( hel0 * pvec3 * emp + hel * pt / ppDENOM * sqh, 0 ); // hack: dummy[ieppV] is not used if pp[ieppV]==0 + // Branch B1: pp != 0. and pt != 0. + volatile fptype_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 + const fptype_v pzpt = pvec3 / ( ppDENOM * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_v vcB1_3 = cxmake( hel0 * pvec1 * emp - pvec1 * pzpt, -(fptype)nsvahl * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_v vcB1_4 = cxmake( hel0 * pvec2 * emp - pvec2 * pzpt, (fptype)nsvahl * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + // Branch B2: pp != 0. and pt == 0. + const cxtype vcB2_3 = cxmake( -hel * sqh, 0. ); + const cxtype_v vcB2_4 = cxmake( 0., (fptype)nsvahl * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + // Choose between the results from branch A and branch B (and from branch B1 and branch B2) + const bool_v mask = ( pp == 0. ); + const bool_v maskB = ( pt != 0. ); + w[0] = cxternary( mask, vcA_2, vcB_2 ); + w[1] = cxternary( mask, vcA_3, cxternary( maskB, vcB1_3, vcB2_3 ) ); + w[2] = cxternary( mask, vcA_4, cxternary( maskB, vcB1_4, vcB2_4 ) ); + w[3] = cxternary( mask, vcA_5, vcB_5 ); + + //FD gauge: same two branches as the scalar code above, selected lane by + //lane. The division uses ppDENOM (=1 where pp==0) so that the lanes that + //do not take it are not poisoned: a select discards the other value, but + //a nan surviving an arithmetic combination (nan*0 is nan) would not be. + const bool_v maskFD = ( pp > zero ); + n[0] = fpternary( pvec0 >= zero , one , -one ); + n[1] = fpternary( maskFD, -pvec1 / ppDENOM, zero ); + n[2] = fpternary( maskFD, -pvec2 / ppDENOM, zero ); + n[3] = fpternary( maskFD, -pvec3 / ppDENOM, -n[0] ); + n[4] = zero; + + nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; + + // nhel is a scalar: no need for a per-lane mask here (and a bool_v built + // from a single value would only set the first lane) + if ( abs(nhel) == 1 ) + { + w[4] = cxzero_sv(); + } + else + { + w[0] = cxmake( -vmass/nk * n[0], zero ); + w[1] = cxmake( -vmass/nk * n[1], zero ); + w[2] = cxmake( -vmass/nk * n[2], zero ); + w[3] = cxmake( -vmass/nk * n[3], zero ); + w[4] = -static_cast(nsv)*cI; + } +#endif + } + else + { + const fptype_sv& pp = pvec0; // NB: rewrite the following as in Fortran, using pp instead of pvec0 +#ifndef MGONGPU_CPPSIMD + const fptype_sv pt = fpsqrt( ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ) ); +#else + volatile fptype_sv pt2 = pvec1 * pvec1 + pvec2 * pvec2; // volatile fixes #736 + const fptype_sv pt = fpsqrt( pt2 ); +#endif + w[0] = cxzero_sv(); + w[3] = cxmake( hel * pt / pp * sqh, 0. ); +#ifndef MGONGPU_CPPSIMD + if( pt != 0. ) + { + const fptype pzpt = pvec3 / ( pp * pt ) * sqh * hel; + w[1] = cxmake( -pvec1 * pzpt, -nsv * pvec2 / pt * sqh ); + w[2] = cxmake( -pvec2 * pzpt, nsv * pvec1 / pt * sqh ); + } + else + { + w[1] = cxmake( -hel * sqh, 0. ); + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + //w[2] = cxmake( 0, nsv * ( pvec3 < 0. ? -std::abs( sqh ) : std::abs( sqh ) ) ); // AV why abs here? + w[2] = cxmake( 0., nsv * ( pvec3 < 0. ? -sqh : sqh ) ); // AV: removed an abs here + } +#else + // Branch A: pt != 0. + volatile fptype_v ptDENOM = fpternary( pt != 0, pt, 1. ); // hack: ptDENOM[ieppV]=1 if pt[ieppV]==0 + const fptype_v pzpt = pvec3 / ( pp * ptDENOM ) * sqh * hel; // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_v vcA_3 = cxmake( -pvec1 * pzpt, -(fptype)nsv * pvec2 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + const cxtype_v vcA_4 = cxmake( -pvec2 * pzpt, (fptype)nsv * pvec1 / ptDENOM * sqh ); // hack: dummy[ieppV] is not used if pt[ieppV]==0 + // Branch B: pt == 0. + const cxtype vcB_3 = cxmake( -(fptype)hel * sqh, 0 ); + const cxtype_v vcB_4 = cxmake( 0, (fptype)nsv * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here + // Choose between the results from branch A and branch B + const bool_v mask = ( pt != 0. ); + w[1] = cxternary( mask, vcA_3, vcB_3 ); + w[2] = cxternary( mask, vcA_4, vcB_4 ); +#endif + //FD gauge + w[4] = cxzero_sv(); + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction sc[3] from the input momenta[npar*4*nevt] + template + __host__ __device__ void + sxxxxx( const fptype momenta[], // input: momenta + //const fptype, // WARNING: input "smass" unused (missing in Fortran) - scalar boson mass + //const int, // WARNING: input "nhel" unused (missing in Fortran) - scalar has no helicity! + const int nss, // input: +1 (final) or -1 (initial) + const int flv, // input: flavour + ALOHAOBJ &sc, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( sc.w ); + + sc.pvec[0] = pvec0 * (fptype)nss; + sc.pvec[1] = pvec1 * (fptype)nss; + sc.pvec[2] = pvec2 * (fptype)nss; + sc.pvec[3] = pvec3 * (fptype)nss; + + sc.flv_index = flv; + w[0] = cxmake( 1 + fptype_sv{ 0 }, 0 ); + //FD gauge + w[1] = cxmake( 0 + fptype_sv{ 0 }, 0 ); + w[2] = cxmake( 0 + fptype_sv{ 0 }, 0 ); + w[3] = cxmake( 0 + fptype_sv{ 0 }, 0 ); + w[4] = cxmake( 1 + fptype_sv{ 0 }, 0 ); + + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + template + __host__ __device__ void + oxxxxx( const fptype momenta[], // input: momenta + const fptype fmass, // input: fermion mass + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + int flv, // input: flavour + ALOHAOBJ & fo, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + // NEW IMPLEMENTATION FIXING FLOATING POINT EXCEPTIONS IN SIMD CODE (#701) + // Variables xxxDENOM are a hack to avoid division-by-0 FPE while preserving speed (#701 and #727) + // Variables xxxDENOM are declared as 'volatile' to make sure they are not optimized away on clang! (#724) + // A few additional variables are declared as 'volatile' to avoid sqrt-of-negative-number FPEs (#736) + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fo.w ); + fo.pvec[0] = pvec0 * (fptype)nsf; + fo.pvec[1] = pvec1 * (fptype)nsf; + fo.pvec[2] = pvec2 * (fptype)nsf; + fo.pvec[3] = pvec3 * (fptype)nsf; + fo.flv_index = flv; + const int nh = nhel * nsf; + if( fmass != 0. ) + { +#ifndef MGONGPU_CPPSIMD + const fptype_sv pp = fpmin( pvec0, fpsqrt( ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ) + ( pvec3 * pvec3 ) ) ); + if( pp == 0. ) + { + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses + //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? + sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here + const int ip = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int im = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + w[0] = cxmake( im * sqm[std::abs( ip )], 0 ); + w[1] = cxmake( ip * nsf * sqm[std::abs( ip )], 0 ); + w[2] = cxmake( im * nsf * sqm[std::abs( im )], 0 ); + w[3] = cxmake( ip * sqm[std::abs( im )], 0 ); + } + else + { + const fptype sf[2] = { fptype( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype)0.5, + fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; + fptype omega[2] = { fpsqrt( pvec0 + pp ), 0. }; + omega[1] = fmass / omega[0]; + const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 + const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 + const fptype sfomeg[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; + const fptype pp3 = fpmax( pp + pvec3, 0. ); + const cxtype chi[2] = { cxmake( fpsqrt( pp3 * (fptype)0.5 / pp ), 0. ), + ( ( pp3 == 0. ) ? cxmake( -nh, 0. ) + : cxmake( nh * pvec1, -pvec2 ) / fpsqrt( 2. * pp * pp3 ) ) }; + w[0] = sfomeg[1] * chi[im]; + w[1] = sfomeg[1] * chi[ip]; + w[2] = sfomeg[0] * chi[im]; + w[3] = sfomeg[0] * chi[ip]; + } +#else + volatile fptype_sv p2 = pvec1 * pvec1 + pvec2 * pvec2 + pvec3 * pvec3; // volatile fixes #736 + const fptype_sv pp = fpmin( pvec0, fpsqrt( p2 ) ); + // Branch A: pp == 0. + // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs! + fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0 }; // possibility of negative fermion masses + sqm[1] = ( fmass < 0 ? -sqm[0] : sqm[0] ); // AV: removed an abs here (as above) + const int ipA = -( ( 1 - nh ) / 2 ) * nhel; + const int imA = ( 1 + nh ) / 2 * nhel; + const cxtype foA_2 = imA * sqm[std::abs( ipA )]; + const cxtype foA_3 = ipA * nsf * sqm[std::abs( ipA )]; + const cxtype foA_4 = imA * nsf * sqm[std::abs( imA )]; + const cxtype foA_5 = ipA * sqm[std::abs( imA )]; + // Branch B: pp != 0. + const fptype sf[2] = { fptype( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype)0.5, + fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; + fptype_v omega[2] = { fpsqrt( pvec0 + pp ), 0 }; + omega[1] = fmass / omega[0]; + const int ipB = ( 1 + nh ) / 2; + const int imB = ( 1 - nh ) / 2; + const fptype_v sfomeg[2] = { sf[0] * omega[ipB], sf[1] * omega[imB] }; + const fptype_v pp3 = fpmax( pp + pvec3, 0. ); + volatile fptype_v ppDENOM = fpternary( pp != 0, pp, 1. ); // hack: ppDENOM[ieppV]=1 if pp[ieppV]==0 + volatile fptype_v pp3DENOM = fpternary( pp3 != 0, pp3, 1. ); // hack: pp3DENOM[ieppV]=1 if pp3[ieppV]==0 + volatile fptype_v chi0r2 = pp3 * 0.5 / ppDENOM; // volatile fixes #736 + const cxtype_v chi[2] = { cxmake( fpsqrt( chi0r2 ), 0. ), // hack: dummy[ieppV] is not used if pp[ieppV]==0 + ( cxternary( ( pp3 == 0. ), + cxmake( -nh, 0. ), + cxmake( (fptype)nh * pvec1, -pvec2 ) / fpsqrt( 2. * ppDENOM * pp3DENOM ) ) ) }; // hack: dummy[ieppV] is not used if pp[ieppV]==0 + const cxtype_v foB_2 = sfomeg[1] * chi[imB]; + const cxtype_v foB_3 = sfomeg[1] * chi[ipB]; + const cxtype_v foB_4 = sfomeg[0] * chi[imB]; + const cxtype_v foB_5 = sfomeg[0] * chi[ipB]; + // Choose between the results from branch A and branch B + const bool_v mask = ( pp == 0. ); + w[0] = cxternary( mask, foA_2, foB_2 ); + w[1] = cxternary( mask, foA_3, foB_3 ); + w[2] = cxternary( mask, foA_4, foB_4 ); + w[3] = cxternary( mask, foA_5, foB_5 ); +#endif + } + else + { +#ifdef MGONGPU_CPPSIMD + volatile fptype_sv p0p3 = fpmax( pvec0 + pvec3, 0 ); // volatile fixes #736 + volatile fptype_sv sqp0p3 = fpternary( ( pvec1 == 0. and pvec2 == 0. and pvec3 < 0. ), + fptype_sv{ 0 }, + fpsqrt( p0p3 ) * (fptype)nsf ); + volatile fptype_v sqp0p3DENOM = fpternary( sqp0p3 != 0, (fptype_sv)sqp0p3, 1. ); // hack: sqp0p3DENOM[ieppV]=1 if sqp0p3[ieppV]==0 + const cxtype_v chi[2] = { cxmake( (fptype_v)sqp0p3, 0. ), + cxternary( ( sqp0p3 == 0. ), + cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ), + cxmake( (fptype)nh * pvec1, -pvec2 ) / (const fptype_sv)sqp0p3DENOM ) }; // hack: dummy[ieppV] is not used if sqp0p3[ieppV]==0 +#else + const fptype_sv sqp0p3 = fpternary( ( pvec1 == 0. ) and ( pvec2 == 0. ) and ( pvec3 < 0. ), + 0, + fpsqrt( fpmax( pvec0 + pvec3, 0. ) ) * (fptype)nsf ); + const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), + ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ) : cxmake( (fptype)nh * pvec1, -pvec2 ) / sqp0p3 ) }; +#endif + if( nh == 1 ) + { + w[0] = chi[0]; + w[1] = chi[1]; + w[2] = cxzero_sv(); + w[3] = cxzero_sv(); + } + else + { + w[0] = cxzero_sv(); + w[1] = cxzero_sv(); + w[2] = chi[1]; + w[3] = chi[0]; + } + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0) + template + __host__ __device__ void + opzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fo, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fo.w ); + fo.pvec[0] = pvec3 * (fptype)nsf; + fo.pvec[1] = fptype_sv{ 0 }; + fo.pvec[2] = fptype_sv{ 0 }; + fo.pvec[3] = pvec3 * (fptype)nsf; + fo.flv_index = flv; + const int nh = nhel * nsf; + const cxtype_sv csqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); + w[1] = cxzero_sv(); + w[2] = cxzero_sv(); + if( nh == 1 ) + { + w[0] = csqp0p3; + w[3] = cxzero_sv(); + } + else + { + w[0] = cxzero_sv(); + w[3] = csqp0p3; + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0) + template + __host__ __device__ void + omzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fo, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fo.w ); + fo.pvec[0] = -pvec3 * (fptype)nsf; + fo.pvec[1] = fptype_sv{ 0 }; + fo.pvec[2] = fptype_sv{ 0 }; + fo.pvec[3] = pvec3 * (fptype)nsf; + fo.flv_index = flv; + const int nh = nhel * nsf; + const cxtype_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -2. * pvec3 ); + if( nh == 1 ) + { + w[0] = cxzero_sv(); + w[1] = chi1; + w[2] = cxzero_sv(); + w[3] = cxzero_sv(); + } + else + { + w[0] = cxzero_sv(); + w[1] = cxzero_sv(); + w[2] = chi1; + //w[3] = chi1; // AV: BUG! + w[3] = cxzero_sv(); // AV: BUG FIX + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt] + // ASSUMPTIONS: (FMASS == 0) and (PT > 0) + template + __host__ __device__ void + oxzxxx( const fptype momenta[], // input: momenta + //const fptype fmass, // [skip: ASSUME fermion mass==0] + const int nhel, // input: -1 or +1 (helicity of fermion) + const int nsf, // input: +1 (particle) or -1 (antiparticle) + const int flv, // input: flavour + ALOHAOBJ & fo, // output: wavefunctions + const int ipar ) // input: particle# out of npar + { + mgDebug( 0, __FUNCTION__ ); + const fptype_sv& pvec0 = M_ACCESS::kernelAccessIp4IparConst( momenta, 0, ipar ); + const fptype_sv& pvec1 = M_ACCESS::kernelAccessIp4IparConst( momenta, 1, ipar ); + const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); + const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); + cxtype_sv* w = W_ACCESS::kernelAccess( fo.w ); + fo.pvec[0] = pvec0 * (fptype)nsf; + fo.pvec[1] = pvec1 * (fptype)nsf; + fo.pvec[2] = pvec2 * (fptype)nsf; + fo.pvec[3] = pvec3 * (fptype)nsf; + fo.flv_index = flv; + const int nh = nhel * nsf; + //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? + const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; + const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); + const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); + if( nh == 1 ) + { + w[0] = chi0; + w[1] = chi1; + w[2] = cxzero_sv(); + w[3] = cxzero_sv(); + } + else + { + w[0] = cxzero_sv(); + w[1] = cxzero_sv(); + w[2] = chi1; + w[3] = chi0; + } + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + // Compute the direction n[5] of the gauge q[5] + // TODO: Utilise pvec instead of the whole q + __host__ __device__ INLINE void + define_gauge_dir( const cxtype_sv q[5], // input: gauge + fptype_sv n[5] ) // output: direction + { + const fptype_sv qabs2 = q[1].real()*q[1].real() + + q[2].real()*q[2].real() + + q[3].real()*q[3].real(); + + const fptype_sv one = 1. + fptype_sv{0}; + const fptype_sv zero = 0. + fptype_sv{0}; + +#ifndef MGONGPU_CPPSIMD + + if (qabs2 > 0.f) + { + const fptype_sv qabs = fpsqrt(qabs2); + + n[0] = fpternary( q[0].real() >= 0.f , one , -one); + n[1] = -q[1].real() / qabs; + n[2] = -q[2].real() / qabs; + n[3] = -q[3].real() / qabs; + n[4] = zero; + } + else + { + n[0] = fpternary( q[0].real() >= 0.f , one , -one ); + n[1] = zero; + n[2] = zero; + n[3] = fpternary( q[0].real() >= 0.f , -one , one); // -sign(q0), as in fortran and python define_gauge_dir + n[4] = zero; + } +#else + const fptype_sv qabs = fpsqrt(qabs2); + const bool_v qsign = (qabs2 > 0.f); + n[0] = fpternary( q[0].real() >= 0.f , one , -one); + n[1] = fpternary( qsign , -q[1].real() / qabs , zero ); + n[2] = fpternary( qsign , -q[2].real() / qabs , zero ); + n[3] = fpternary( qsign , -q[3].real() / qabs , fpternary( q[0].real() >= 0.f , -one , one)); // same gauge as the branch above + n[4] = zero; +#endif + } + +//-------------------------------------------------------------------------- +// Compute propagator factor d from the gauge q[5] and mass + __host__ __device__ INLINE void + calculate_propagator_factor( const cxtype_sv q[5], // input: gauge + const fptype mass, // input: mass + fptype_sv *d ) // output: propagator factor + { + const fptype_sv one = 1. + fptype_sv{0}; + const fptype_sv q2 = q[0].real()*q[0].real() - ( q[1].real()*q[1].real() + q[2].real()*q[2].real() + q[3].real()*q[3].real() ); + *d = one / (q2 - mass*mass); + } + +//-------------------------------------------------------------------------- +// Multiply the wavefunction by propagator factor from momenta and m +// TODO: check if d should not be used + template< class W_ACCESS> + __host__ __device__ INLINE void + multiply_propagator_factor( const ALOHAOBJ & Ain, // input: wavefunctions + const fptype m, // input: mass + ALOHAOBJ Aout ) // output: wavefunctions + { + + const cxtype_sv* win = W_ACCESS::kernelAccessConst( Ain.w ); + cxtype_sv* wout = W_ACCESS::kernelAccess( Aout.w ); + + cxtype_sv q[5]; + fptype_sv n[5]; + cxtype_sv w0[5], w1[5]; + + const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1. + fptype_sv{ 0 } ); + + // Construct q from momenta + q[0] = cxmake( -Ain.pvec[0], 0.); + q[1] = cxmake( -Ain.pvec[1], 0.); + q[2] = cxmake( -Ain.pvec[2], 0.); + q[3] = cxmake( -Ain.pvec[3], 0.); + q[4] = -cI*m; + + // Copy the momenta + Aout.pvec[0] = Ain.pvec[0]; + Aout.pvec[0] = Ain.pvec[0]; + Aout.pvec[0] = Ain.pvec[0]; + Aout.pvec[0] = Ain.pvec[0]; + + define_gauge_dir(q, n); + + w0[0] = win[0]; + w0[1] = win[1]; + w0[2] = win[2]; + w0[3] = win[3]; + w0[4] = win[4]; + + fptype_sv nq = + n[0]*q[0].real() + - n[1]*q[1].real() + - n[2]*q[2].real() + - n[3]*q[3].real(); + + cxtype_sv js1 = + ( n[0]*w0[0] + - n[1]*w0[1] + - n[2]*w0[2] + - n[3]*w0[3] ) / nq; + + cxtype_sv js2 = + ( q[0]*w0[0] + - q[1]*w0[1] + - q[2]*w0[2] + - q[3]*w0[3] + - cxconj(q[4]) * w0[4] ) / nq; + + w1[0] = w0[0] - q[0]*js1 - n[0]*js2; + w1[1] = w0[1] - q[1]*js1 - n[1]*js2; + w1[2] = w0[2] - q[2]*js1 - n[2]*js2; + w1[3] = w0[3] - q[3]*js1 - n[3]*js2; + w1[4] = w0[4] - q[4]*js1 - n[4]*js2; + + wout[0] = w1[0]; + wout[1] = w1[1]; + wout[2] = w1[2]; + wout[3] = w1[3]; + wout[4] = w1[4]; + } + //-------------------------------------------------------------------------- + //========================================================================== + + // Compute the output wavefunction 'V3[6]' from the input wavefunctions + template + __device__ INLINE void + FFV1MP0_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction 'V3[6]' from the input wavefunctions + template + __device__ INLINE void + FFV6M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + // Compute the output wavefunction 'V3[6]' from the input wavefunctions + template + __device__ INLINE void + FFV6_2M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP1, + const double Ccoeff1, + const FLV_COUPLING_VIEW &MCOUP2, + const double Ccoeff2, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + FFV2M_0( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const ALOHAOBJ & V3, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction 'F2[6]' from the input wavefunctions + template + __device__ INLINE void + FFV2M_2( const ALOHAOBJ & F1, + const ALOHAOBJ & V3, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M2, + const fptype & W2, + ALOHAOBJ & F2 ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction 'V3[6]' from the input wavefunctions + template + __device__ INLINE void + FFV2M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VVV1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VVV1_VVS1_VSV2_VSS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP1[], + const double Ccoeff1, + const fptype allCOUP2[], + const double Ccoeff2, + const fptype allCOUP3[], + const double Ccoeff3, + const fptype allCOUP4[], + const double Ccoeff4, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VVV1_VSV2_VSS2_SVV2_SVS2_SSV3_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP1[], + const double Ccoeff1, + const fptype allCOUP2[], + const double Ccoeff2, + const fptype allCOUP3[], + const double Ccoeff3, + const fptype allCOUP4[], + const double Ccoeff4, + const fptype allCOUP5[], + const double Ccoeff5, + const fptype allCOUP6[], + const double Ccoeff6, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VVS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VSV2_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VSS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + VSS2_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + SVV2_0( const ALOHAOBJ & S1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + SVS2_0( const ALOHAOBJ & S1, + const ALOHAOBJ & V2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ INLINE void + SSV3_0( const ALOHAOBJ & S1, + const ALOHAOBJ & S2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) ALWAYS_INLINE; + + //========================================================================== + + // Compute the output wavefunction 'V3[6]' from the input wavefunctions + template + __device__ void + FFV1MP0_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wF2 = W_ACCESS::kernelAccessConst( F2.w ); + cxtype_sv COUP; + cxtype_sv* wV3 = W_ACCESS::kernelAccess( V3.w ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + V3.pvec[0] = +F1.pvec[0] + F2.pvec[0]; + V3.pvec[1] = +F1.pvec[1] + F2.pvec[1]; + V3.pvec[2] = +F1.pvec[2] + F2.pvec[2]; + V3.pvec[3] = +F1.pvec[3] + F2.pvec[3]; + const fptype_sv P3[4] = { -V3.pvec[0], -V3.pvec[1], -V3.pvec[2], -V3.pvec[3] }; + wV3[0] = CZERO ; + wV3[1] = CZERO ; + wV3[2] = CZERO ; + wV3[3] = CZERO ; + wV3[4] = CZERO ; + cxtype_sv FDQ[5] = { cxmake( -V3.pvec[0], 0. ), cxmake( -V3.pvec[1], 0. ), cxmake( -V3.pvec[2], 0. ), cxmake( -V3.pvec[3], 0. ), cxmake( fptype_sv{ 0 }, -M3 + fptype_sv{ 0 } ) }; + fptype_sv FDN[5]; + define_gauge_dir( FDQ, FDN ); + const fptype_sv FDNQ = FDN[0] * FDQ[0].real() - FDN[1] * FDQ[1].real() - FDN[2] * FDQ[2].real() - FDN[3] * FDQ[3].real(); + const int & flv_index1 = F1.flv_index; + const int & flv_index2 = F2.flv_index; + if(flv_index1 == -1 || flv_index2 == -1) { + for(int i=0; i + __device__ void + FFV6M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wF2 = W_ACCESS::kernelAccessConst( F2.w ); + cxtype_sv COUP; + cxtype_sv* wV3 = W_ACCESS::kernelAccess( V3.w ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + V3.pvec[0] = +F1.pvec[0] + F2.pvec[0]; + V3.pvec[1] = +F1.pvec[1] + F2.pvec[1]; + V3.pvec[2] = +F1.pvec[2] + F2.pvec[2]; + V3.pvec[3] = +F1.pvec[3] + F2.pvec[3]; + const fptype_sv P3[4] = { -V3.pvec[0], -V3.pvec[1], -V3.pvec[2], -V3.pvec[3] }; + wV3[0] = CZERO ; + wV3[1] = CZERO ; + wV3[2] = CZERO ; + wV3[3] = CZERO ; + wV3[4] = CZERO ; + cxtype_sv FDQ[5] = { cxmake( -V3.pvec[0], 0. ), cxmake( -V3.pvec[1], 0. ), cxmake( -V3.pvec[2], 0. ), cxmake( -V3.pvec[3], 0. ), cxmake( fptype_sv{ 0 }, -M3 + fptype_sv{ 0 } ) }; + fptype_sv FDN[5]; + define_gauge_dir( FDQ, FDN ); + const fptype_sv FDNQ = FDN[0] * FDQ[0].real() - FDN[1] * FDQ[1].real() - FDN[2] * FDQ[2].real() - FDN[3] * FDQ[3].real(); + const int & flv_index1 = F1.flv_index; + const int & flv_index2 = F2.flv_index; + if(flv_index1 == -1 || flv_index2 == -1) { + for(int i=0; i + __device__ void + FFV6_2M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP1, + const double Ccoeff1, + const FLV_COUPLING_VIEW &MCOUP2, + const double Ccoeff2, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wF2 = W_ACCESS::kernelAccessConst( F2.w ); + cxtype_sv COUP1; + cxtype_sv COUP2; + cxtype_sv* wV3 = W_ACCESS::kernelAccess( V3.w ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + V3.pvec[0] = +F1.pvec[0] + F2.pvec[0]; + V3.pvec[1] = +F1.pvec[1] + F2.pvec[1]; + V3.pvec[2] = +F1.pvec[2] + F2.pvec[2]; + V3.pvec[3] = +F1.pvec[3] + F2.pvec[3]; + const fptype_sv P3[4] = { -V3.pvec[0], -V3.pvec[1], -V3.pvec[2], -V3.pvec[3] }; + wV3[0] = CZERO ; + wV3[1] = CZERO ; + wV3[2] = CZERO ; + wV3[3] = CZERO ; + wV3[4] = CZERO ; + cxtype_sv FDQ[5] = { cxmake( -V3.pvec[0], 0. ), cxmake( -V3.pvec[1], 0. ), cxmake( -V3.pvec[2], 0. ), cxmake( -V3.pvec[3], 0. ), cxmake( fptype_sv{ 0 }, -M3 + fptype_sv{ 0 } ) }; + fptype_sv FDN[5]; + define_gauge_dir( FDQ, FDN ); + const fptype_sv FDNQ = FDN[0] * FDQ[0].real() - FDN[1] * FDQ[1].real() - FDN[2] * FDQ[2].real() - FDN[3] * FDQ[3].real(); + const int & flv_index1 = F1.flv_index; + const int & flv_index2 = F2.flv_index; + int zero_coup1 = 0; + int zero_coup2 = 0; + if(flv_index1 != flv_index2 || flv_index1 == -1) { + for(int i=0; i + __device__ void + FFV2M_0( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const ALOHAOBJ & V3, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wF2 = W_ACCESS::kernelAccessConst( F2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + cxtype_sv COUP; + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const int & flv_index1 = F1.flv_index; + const int & flv_index2 = F2.flv_index; + if(flv_index1 == -1 || flv_index2 == -1) { + *vertex = cxzero_sv(); + return; + } + int flv_sel = -1; + if(MCOUP.partner1[flv_index1] == flv_index2) flv_sel = flv_index1; + else if(MCOUP.partner1[flv_index2] == flv_index1) flv_sel = flv_index2; + if(flv_sel == -1) { + *vertex = cxzero_sv(); + return; + } + COUP = C_ACCESS::kernelAccessConst( MCOUP.value + C_ACCESS::flv_stride*flv_sel ); + const cxtype_sv TMP0 = ( wF1[0] * ( wF2[2] * ( wV3[0] + wV3[3] ) + wF2[3] * ( wV3[1] + cI * wV3[2] ) ) + wF1[1] * ( wF2[2] * ( wV3[1] - cI * wV3[2] ) + wF2[3] * ( wV3[0] - wV3[3] ) ) ); + ( *vertex ) = Ccoeff * COUP * -cI * TMP0; + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output wavefunction 'F2[6]' from the input wavefunctions + template + __device__ void + FFV2M_2( const ALOHAOBJ & F1, + const ALOHAOBJ & V3, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M2, + const fptype & W2, + ALOHAOBJ & F2 ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + cxtype_sv COUP; + cxtype_sv* wF2 = W_ACCESS::kernelAccess( F2.w ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + F2.pvec[0] = +F1.pvec[0] + V3.pvec[0]; + F2.pvec[1] = +F1.pvec[1] + V3.pvec[1]; + F2.pvec[2] = +F1.pvec[2] + V3.pvec[2]; + F2.pvec[3] = +F1.pvec[3] + V3.pvec[3]; + const fptype_sv P2[4] = { -F2.pvec[0], -F2.pvec[1], -F2.pvec[2], -F2.pvec[3] }; + int flv_index1 = F1.flv_index; + if(flv_index1 == -1) { + for(int i=0; i + __device__ void + FFV2M_3( const ALOHAOBJ & F1, + const ALOHAOBJ & F2, + const FLV_COUPLING_VIEW &MCOUP, + const double Ccoeff, + const fptype & M3, + const fptype & W3, + ALOHAOBJ & V3 ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wF1 = W_ACCESS::kernelAccessConst( F1.w ); + const cxtype_sv* wF2 = W_ACCESS::kernelAccessConst( F2.w ); + cxtype_sv COUP; + cxtype_sv* wV3 = W_ACCESS::kernelAccess( V3.w ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + V3.pvec[0] = +F1.pvec[0] + F2.pvec[0]; + V3.pvec[1] = +F1.pvec[1] + F2.pvec[1]; + V3.pvec[2] = +F1.pvec[2] + F2.pvec[2]; + V3.pvec[3] = +F1.pvec[3] + F2.pvec[3]; + const fptype_sv P3[4] = { -V3.pvec[0], -V3.pvec[1], -V3.pvec[2], -V3.pvec[3] }; + wV3[0] = CZERO ; + wV3[1] = CZERO ; + wV3[2] = CZERO ; + wV3[3] = CZERO ; + wV3[4] = CZERO ; + cxtype_sv FDQ[5] = { cxmake( -V3.pvec[0], 0. ), cxmake( -V3.pvec[1], 0. ), cxmake( -V3.pvec[2], 0. ), cxmake( -V3.pvec[3], 0. ), cxmake( fptype_sv{ 0 }, -M3 + fptype_sv{ 0 } ) }; + fptype_sv FDN[5]; + define_gauge_dir( FDQ, FDN ); + const fptype_sv FDNQ = FDN[0] * FDQ[0].real() - FDN[1] * FDQ[1].real() - FDN[2] * FDQ[2].real() - FDN[3] * FDQ[3].real(); + const int & flv_index1 = F1.flv_index; + const int & flv_index2 = F2.flv_index; + if(flv_index1 == -1 || flv_index2 == -1) { + for(int i=0; i + __device__ void + VVV1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P1[4] = { +V1.pvec[0], +V1.pvec[1], +V1.pvec[2], +V1.pvec[3] }; + const fptype_sv P2[4] = { +V2.pvec[0], +V2.pvec[1], +V2.pvec[2], +V2.pvec[3] }; + const fptype_sv P3[4] = { +V3.pvec[0], +V3.pvec[1], +V3.pvec[2], +V3.pvec[3] }; + const cxtype_sv TMP1 = ( wV2[0] * wV1[0] - wV2[1] * wV1[1] - wV2[2] * wV1[2] - wV2[3] * wV1[3] ); + const cxtype_sv TMP2 = ( wV3[0] * P1[0] - wV3[1] * P1[1] - wV3[2] * P1[2] - wV3[3] * P1[3] ); + const cxtype_sv TMP3 = ( wV3[0] * P2[0] - wV3[1] * P2[1] - wV3[2] * P2[2] - wV3[3] * P2[3] ); + const cxtype_sv TMP4 = ( wV2[0] * P1[0] - wV2[1] * P1[1] - wV2[2] * P1[2] - wV2[3] * P1[3] ); + const cxtype_sv TMP5 = ( wV3[0] * wV1[0] - wV3[1] * wV1[1] - wV3[2] * wV1[2] - wV3[3] * wV1[3] ); + const cxtype_sv TMP6 = ( wV2[0] * P3[0] - wV2[1] * P3[1] - wV2[2] * P3[2] - wV2[3] * P3[3] ); + const cxtype_sv TMP7 = ( wV3[0] * wV2[0] - wV3[1] * wV2[1] - wV3[2] * wV2[2] - wV3[3] * wV2[3] ); + const cxtype_sv TMP8 = ( P2[0] * wV1[0] - P2[1] * wV1[1] - P2[2] * wV1[2] - P2[3] * wV1[3] ); + const cxtype_sv TMP9 = ( wV1[0] * P3[0] - wV1[1] * P3[1] - wV1[2] * P3[2] - wV1[3] * P3[3] ); + ( *vertex ) = Ccoeff * COUP * ( TMP1 * ( -cI * TMP2 + cI * TMP3 ) + ( TMP5 * ( +cI * TMP4 - cI * TMP6 ) + TMP7 * ( -cI * TMP8 + cI * TMP9 ) ) ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VVV1_VVS1_VSV2_VSS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP1[], + const double Ccoeff1, + const fptype allCOUP2[], + const double Ccoeff2, + const fptype allCOUP3[], + const double Ccoeff3, + const fptype allCOUP4[], + const double Ccoeff4, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP1 = C_ACCESS::kernelAccessConst( allCOUP1 ); + const cxtype_sv COUP2 = C_ACCESS::kernelAccessConst( allCOUP2 ); + const cxtype_sv COUP3 = C_ACCESS::kernelAccessConst( allCOUP3 ); + const cxtype_sv COUP4 = C_ACCESS::kernelAccessConst( allCOUP4 ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P1[4] = { +V1.pvec[0], +V1.pvec[1], +V1.pvec[2], +V1.pvec[3] }; + const fptype_sv P2[4] = { +V2.pvec[0], +V2.pvec[1], +V2.pvec[2], +V2.pvec[3] }; + const fptype_sv P3[4] = { +V3.pvec[0], +V3.pvec[1], +V3.pvec[2], +V3.pvec[3] }; + ( *vertex ) = cxzero_sv(); + const cxtype_sv TMP1 = ( wV2[0] * wV1[0] - wV2[1] * wV1[1] - wV2[2] * wV1[2] - wV2[3] * wV1[3] ); + const cxtype_sv TMP2 = ( wV3[0] * P1[0] - wV3[1] * P1[1] - wV3[2] * P1[2] - wV3[3] * P1[3] ); + const cxtype_sv TMP3 = ( wV3[0] * P2[0] - wV3[1] * P2[1] - wV3[2] * P2[2] - wV3[3] * P2[3] ); + const cxtype_sv TMP4 = ( wV2[0] * P1[0] - wV2[1] * P1[1] - wV2[2] * P1[2] - wV2[3] * P1[3] ); + const cxtype_sv TMP5 = ( wV3[0] * wV1[0] - wV3[1] * wV1[1] - wV3[2] * wV1[2] - wV3[3] * wV1[3] ); + const cxtype_sv TMP6 = ( wV2[0] * P3[0] - wV2[1] * P3[1] - wV2[2] * P3[2] - wV2[3] * P3[3] ); + const cxtype_sv TMP7 = ( wV3[0] * wV2[0] - wV3[1] * wV2[1] - wV3[2] * wV2[2] - wV3[3] * wV2[3] ); + const cxtype_sv TMP8 = ( P2[0] * wV1[0] - P2[1] * wV1[1] - P2[2] * wV1[2] - P2[3] * wV1[3] ); + const cxtype_sv TMP9 = ( wV1[0] * P3[0] - wV1[1] * P3[1] - wV1[2] * P3[2] - wV1[3] * P3[3] ); + ( *vertex ) = ( *vertex ) + Ccoeff1 * COUP1 * ( TMP1 * ( -cI * TMP2 + cI * TMP3 ) + ( TMP5 * ( +cI * TMP4 - cI * TMP6 ) + TMP7 * ( -cI * TMP8 + cI * TMP9 ) ) ); + ( *vertex ) = ( *vertex ) + Ccoeff2 * COUP2 * -cI * TMP1 * wV3[4]; + ( *vertex ) = ( *vertex ) + Ccoeff3 * COUP3 * -cI * TMP5 * wV2[4]; + ( *vertex ) = ( *vertex ) + Ccoeff4 * COUP4 * wV2[4] * wV3[4] * ( -cI * TMP8 + cI * TMP9 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VVV1_VSV2_VSS2_SVV2_SVS2_SSV3_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP1[], + const double Ccoeff1, + const fptype allCOUP2[], + const double Ccoeff2, + const fptype allCOUP3[], + const double Ccoeff3, + const fptype allCOUP4[], + const double Ccoeff4, + const fptype allCOUP5[], + const double Ccoeff5, + const fptype allCOUP6[], + const double Ccoeff6, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP1 = C_ACCESS::kernelAccessConst( allCOUP1 ); + const cxtype_sv COUP2 = C_ACCESS::kernelAccessConst( allCOUP2 ); + const cxtype_sv COUP3 = C_ACCESS::kernelAccessConst( allCOUP3 ); + const cxtype_sv COUP4 = C_ACCESS::kernelAccessConst( allCOUP4 ); + const cxtype_sv COUP5 = C_ACCESS::kernelAccessConst( allCOUP5 ); + const cxtype_sv COUP6 = C_ACCESS::kernelAccessConst( allCOUP6 ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P1[4] = { +V1.pvec[0], +V1.pvec[1], +V1.pvec[2], +V1.pvec[3] }; + const fptype_sv P2[4] = { +V2.pvec[0], +V2.pvec[1], +V2.pvec[2], +V2.pvec[3] }; + const fptype_sv P3[4] = { +V3.pvec[0], +V3.pvec[1], +V3.pvec[2], +V3.pvec[3] }; + ( *vertex ) = cxzero_sv(); + const cxtype_sv TMP1 = ( wV2[0] * wV1[0] - wV2[1] * wV1[1] - wV2[2] * wV1[2] - wV2[3] * wV1[3] ); + const cxtype_sv TMP2 = ( wV3[0] * P1[0] - wV3[1] * P1[1] - wV3[2] * P1[2] - wV3[3] * P1[3] ); + const cxtype_sv TMP3 = ( wV3[0] * P2[0] - wV3[1] * P2[1] - wV3[2] * P2[2] - wV3[3] * P2[3] ); + const cxtype_sv TMP4 = ( wV2[0] * P1[0] - wV2[1] * P1[1] - wV2[2] * P1[2] - wV2[3] * P1[3] ); + const cxtype_sv TMP5 = ( wV3[0] * wV1[0] - wV3[1] * wV1[1] - wV3[2] * wV1[2] - wV3[3] * wV1[3] ); + const cxtype_sv TMP6 = ( wV2[0] * P3[0] - wV2[1] * P3[1] - wV2[2] * P3[2] - wV2[3] * P3[3] ); + const cxtype_sv TMP7 = ( wV3[0] * wV2[0] - wV3[1] * wV2[1] - wV3[2] * wV2[2] - wV3[3] * wV2[3] ); + const cxtype_sv TMP8 = ( P2[0] * wV1[0] - P2[1] * wV1[1] - P2[2] * wV1[2] - P2[3] * wV1[3] ); + const cxtype_sv TMP9 = ( wV1[0] * P3[0] - wV1[1] * P3[1] - wV1[2] * P3[2] - wV1[3] * P3[3] ); + ( *vertex ) = ( *vertex ) + Ccoeff1 * COUP1 * ( TMP1 * ( -cI * TMP2 + cI * TMP3 ) + ( TMP5 * ( +cI * TMP4 - cI * TMP6 ) + TMP7 * ( -cI * TMP8 + cI * TMP9 ) ) ); + ( *vertex ) = ( *vertex ) + Ccoeff2 * COUP2 * -cI * TMP5 * wV2[4]; + ( *vertex ) = ( *vertex ) + Ccoeff3 * COUP3 * wV2[4] * wV3[4] * ( -cI * TMP9 + cI * TMP8 ); + ( *vertex ) = ( *vertex ) + Ccoeff4 * COUP4 * -cI * TMP7 * wV1[4]; + ( *vertex ) = ( *vertex ) + Ccoeff5 * COUP5 * wV1[4] * wV3[4] * ( -cI * TMP6 + cI * TMP4 ); + ( *vertex ) = ( *vertex ) + Ccoeff6 * COUP6 * wV1[4] * wV2[4] * ( -cI * TMP2 + cI * TMP3 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VVS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & V2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wS3 = W_ACCESS::kernelAccessConst( S3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const cxtype_sv TMP1 = ( wV2[0] * wV1[0] - wV2[1] * wV1[1] - wV2[2] * wV1[2] - wV2[3] * wV1[3] ); + ( *vertex ) = Ccoeff * COUP * -cI * TMP1 * wS3[4]; + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VSV2_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wS2 = W_ACCESS::kernelAccessConst( S2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const cxtype_sv TMP5 = ( wV3[0] * wV1[0] - wV3[1] * wV1[1] - wV3[2] * wV1[2] - wV3[3] * wV1[3] ); + ( *vertex ) = Ccoeff * COUP * -cI * TMP5 * wS2[4]; + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VSS1_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wS2 = W_ACCESS::kernelAccessConst( S2.w ); + const cxtype_sv* wS3 = W_ACCESS::kernelAccessConst( S3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P2[4] = { +S2.pvec[0], +S2.pvec[1], +S2.pvec[2], +S2.pvec[3] }; + const fptype_sv P3[4] = { +S3.pvec[0], +S3.pvec[1], +S3.pvec[2], +S3.pvec[3] }; + const cxtype_sv TMP8 = ( P2[0] * wV1[0] - P2[1] * wV1[1] - P2[2] * wV1[2] - P2[3] * wV1[3] ); + const cxtype_sv TMP9 = ( wV1[0] * P3[0] - wV1[1] * P3[1] - wV1[2] * P3[2] - wV1[3] * P3[3] ); + ( *vertex ) = Ccoeff * COUP * wS2[4] * wS3[4] * ( -cI * TMP8 + cI * TMP9 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + VSS2_0( const ALOHAOBJ & V1, + const ALOHAOBJ & S2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wV1 = W_ACCESS::kernelAccessConst( V1.w ); + const cxtype_sv* wS2 = W_ACCESS::kernelAccessConst( S2.w ); + const cxtype_sv* wS3 = W_ACCESS::kernelAccessConst( S3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P2[4] = { +S2.pvec[0], +S2.pvec[1], +S2.pvec[2], +S2.pvec[3] }; + const fptype_sv P3[4] = { +S3.pvec[0], +S3.pvec[1], +S3.pvec[2], +S3.pvec[3] }; + const cxtype_sv TMP8 = ( P2[0] * wV1[0] - P2[1] * wV1[1] - P2[2] * wV1[2] - P2[3] * wV1[3] ); + const cxtype_sv TMP9 = ( wV1[0] * P3[0] - wV1[1] * P3[1] - wV1[2] * P3[2] - wV1[3] * P3[3] ); + ( *vertex ) = Ccoeff * COUP * wS2[4] * wS3[4] * ( -cI * TMP9 + cI * TMP8 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + SVV2_0( const ALOHAOBJ & S1, + const ALOHAOBJ & V2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wS1 = W_ACCESS::kernelAccessConst( S1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const cxtype_sv TMP7 = ( wV3[0] * wV2[0] - wV3[1] * wV2[1] - wV3[2] * wV2[2] - wV3[3] * wV2[3] ); + ( *vertex ) = Ccoeff * COUP * -cI * TMP7 * wS1[4]; + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + SVS2_0( const ALOHAOBJ & S1, + const ALOHAOBJ & V2, + const ALOHAOBJ & S3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wS1 = W_ACCESS::kernelAccessConst( S1.w ); + const cxtype_sv* wV2 = W_ACCESS::kernelAccessConst( V2.w ); + const cxtype_sv* wS3 = W_ACCESS::kernelAccessConst( S3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P1[4] = { +S1.pvec[0], +S1.pvec[1], +S1.pvec[2], +S1.pvec[3] }; + const fptype_sv P3[4] = { +S3.pvec[0], +S3.pvec[1], +S3.pvec[2], +S3.pvec[3] }; + const cxtype_sv TMP4 = ( wV2[0] * P1[0] - wV2[1] * P1[1] - wV2[2] * P1[2] - wV2[3] * P1[3] ); + const cxtype_sv TMP6 = ( wV2[0] * P3[0] - wV2[1] * P3[1] - wV2[2] * P3[2] - wV2[3] * P3[3] ); + ( *vertex ) = Ccoeff * COUP * wS1[4] * wS3[4] * ( -cI * TMP6 + cI * TMP4 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + + // Compute the output amplitude 'vertex' from the input wavefunctions + template + __device__ void + SSV3_0( const ALOHAOBJ & S1, + const ALOHAOBJ & S2, + const ALOHAOBJ & V3, + const fptype allCOUP[], + const double Ccoeff, + fptype allvertexes[] ) + { + mgDebug( 0, __FUNCTION__ ); + const cxtype_sv* wS1 = W_ACCESS::kernelAccessConst( S1.w ); + const cxtype_sv* wS2 = W_ACCESS::kernelAccessConst( S2.w ); + const cxtype_sv* wV3 = W_ACCESS::kernelAccessConst( V3.w ); + const cxtype_sv COUP = C_ACCESS::kernelAccessConst( allCOUP ); + cxtype_sv* vertex = A_ACCESS::kernelAccess( allvertexes ); + cxtype_sv CZERO=cxzero_sv(); + const cxtype cI = cxmake( 0., 1. ); + const fptype_sv P1[4] = { +S1.pvec[0], +S1.pvec[1], +S1.pvec[2], +S1.pvec[3] }; + const fptype_sv P2[4] = { +S2.pvec[0], +S2.pvec[1], +S2.pvec[2], +S2.pvec[3] }; + const cxtype_sv TMP2 = ( wV3[0] * P1[0] - wV3[1] * P1[1] - wV3[2] * P1[2] - wV3[3] * P1[3] ); + const cxtype_sv TMP3 = ( wV3[0] * P2[0] - wV3[1] * P2[1] - wV3[2] * P2[2] - wV3[3] * P2[3] ); + ( *vertex ) = Ccoeff * COUP * wS1[4] * wS2[4] * ( -cI * TMP2 + cI * TMP3 ); + mgDebug( 1, __FUNCTION__ ); + return; + } + + //-------------------------------------------------------------------------- + +} // end namespace + +#endif // HelAmps_sm_H diff --git a/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%FFV6M_3.f b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%FFV6M_3.f new file mode 100644 index 000000000..d1d174302 --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%FFV6M_3.f @@ -0,0 +1,148 @@ +C This File is Automatically generated by ALOHA +C The process calculated in this file is: +C Gamma(3,2,-1)*ProjP(-1,1) +C + SUBROUTINE FFV6M_3(F1, F2, MCOUP, M3, W3,V3) + USE ALOHA_OBJECT + USE MODEL_OBJECT + IMPLICIT NONE + COMPLEX*16 CI + PARAMETER (CI=(0D0,1D0)) + COMPLEX*16 CZERO + PARAMETER (CZERO=(0D0,0D0)) + TYPE(FLV_COUPLING) MCOUP + DOUBLE COMPLEX COUP + INTEGER FLV_INDEX + TYPE(ALOHA) F1 + INTEGER FLV_INDEX1 + TYPE(ALOHA) F2 + INTEGER FLV_INDEX2 + COMPLEX*16 FDJS1 + COMPLEX*16 FDJS2 + REAL*8 FDN(0:4) + REAL*8 FDNQ + COMPLEX*16 FDQ(0:4) + REAL*8 M3 + REAL*8 P3(0:3) + TYPE(ALOHA) V3 + REAL*8 W3 + COMPLEX*16 DENOM + INTEGER*4 I + V3%P(:) = +F1%P(:)+F2%P(:) + P3(:) = -V3 % P (:) + V3 % W(:) = CZERO + FDQ(0:3) = -V3 % P(:) + FDQ(4) = -CI*M3 + CALL DEFINE_GAUGE_DIR(FDQ, FDN) + FDNQ = FDN(0)*DBLE(FDQ(0))-FDN(1)*DBLE(FDQ(1))-FDN(2)*DBLE(FDQ(2) + $ )-FDN(3)*DBLE(FDQ(3)) + FLV_INDEX1 = F1 %FLV_INDEX + FLV_INDEX2 = F2 %FLV_INDEX + IF(FLV_INDEX1.EQ.0.OR.FLV_INDEX2.EQ.0)THEN + V3%W(:) = (0D0,0D0) + RETURN + ENDIF + IF(MCOUP % PARTNER(FLV_INDEX1).NE.FLV_INDEX2)THEN + V3%W(:) = (0D0,0D0) + RETURN + ENDIF + COUP = MCOUP % VAL(FLV_INDEX1) % P + DENOM = COUP/(P3(0)**2-P3(1)**2-P3(2)**2-P3(3)**2 - M3 * (M3 -CI + $ * W3)) + V3%W(1)= DENOM*(-CI)*(F2 % W(1)*F1 % W(3)+F2 % W(2)*F1 % W(4)) + V3%W(2)= DENOM*(-CI)*(F2 % W(2)*F1 % W(3)+F2 % W(1)*F1 % W(4)) + V3%W(3)= DENOM*(-CI)*(+CI*(F2 % W(2)*F1 % W(3))-CI*(F2 % W(1)*F1 + $ % W(4))) + V3%W(4)= DENOM*(-CI)*(F2 % W(1)*F1 % W(3)-F2 % W(2)*F1 % W(4)) + FDJS1 = (FDN(0)*V3%W(1)-FDN(1)*V3%W(2)-FDN(2)*V3%W(3)-FDN(3) + $ *V3%W(4))/FDNQ + FDJS2 = (FDQ(0)*V3%W(1)-FDQ(1)*V3%W(2)-FDQ(2)*V3%W(3)-FDQ(3) + $ *V3%W(4)-DCONJG(FDQ(4))*V3%W(5))/FDNQ + V3%W(1) = V3%W(1)-FDQ(0)*FDJS1-FDN(0)*FDJS2 + V3%W(2) = V3%W(2)-FDQ(1)*FDJS1-FDN(1)*FDJS2 + V3%W(3) = V3%W(3)-FDQ(2)*FDJS1-FDN(2)*FDJS2 + V3%W(4) = V3%W(4)-FDQ(3)*FDJS1-FDN(3)*FDJS2 + V3%W(5) = V3%W(5)-FDQ(4)*FDJS1-FDN(4)*FDJS2 + END + + +C This File is Automatically generated by ALOHA +C The process calculated in this file is: +C Coup(1) * (Gamma(3,2,-1)*ProjP(-1,1)) + Coup(2) * +C (Gamma(3,2,-1)*ProjM(-1,1)) +C + SUBROUTINE FFV6_2M_3(F1, F2, MCOUP1, MCOUP2, M3, W3,V3) + USE ALOHA_OBJECT + USE MODEL_OBJECT + IMPLICIT NONE + COMPLEX*16 CI + PARAMETER (CI=(0D0,1D0)) + COMPLEX*16 CZERO + PARAMETER (CZERO=(0D0,0D0)) + TYPE(FLV_COUPLING) MCOUP1 + DOUBLE COMPLEX COUP1 + INTEGER FLV_INDEX + TYPE(FLV_COUPLING) MCOUP2 + DOUBLE COMPLEX COUP2 + TYPE(ALOHA) F1 + INTEGER FLV_INDEX1 + TYPE(ALOHA) F2 + INTEGER FLV_INDEX2 + COMPLEX*16 FDJS1 + COMPLEX*16 FDJS2 + REAL*8 FDN(0:4) + REAL*8 FDNQ + COMPLEX*16 FDQ(0:4) + REAL*8 M3 + REAL*8 P3(0:3) + TYPE(ALOHA) V3 + REAL*8 W3 + COMPLEX*16 DENOM + INTEGER*4 I + V3%P(:) = +F1%P(:)+F2%P(:) + P3(:) = -V3 % P (:) + V3 % W(:) = CZERO + FDQ(0:3) = -V3 % P(:) + FDQ(4) = -CI*M3 + CALL DEFINE_GAUGE_DIR(FDQ, FDN) + FDNQ = FDN(0)*DBLE(FDQ(0))-FDN(1)*DBLE(FDQ(1))-FDN(2)*DBLE(FDQ(2) + $ )-FDN(3)*DBLE(FDQ(3)) + FLV_INDEX1 = F1 %FLV_INDEX + FLV_INDEX2 = F2 %FLV_INDEX + IF(FLV_INDEX1.EQ.0.OR.FLV_INDEX2.EQ.0)THEN + V3%W(:) = (0D0,0D0) + RETURN + ENDIF + IF(MCOUP1 % PARTNER(FLV_INDEX1).NE.FLV_INDEX2.AND.MCOUP2 % + $ PARTNER(FLV_INDEX1).NE.FLV_INDEX2)THEN + V3%W(:) = (0D0,0D0) + RETURN + ENDIF + COUP1 = MCOUP1 % VAL(FLV_INDEX1) % P + COUP2 = (0D0,0D0) + IF(MCOUP2 % PARTNER(FLV_INDEX1).EQ.FLV_INDEX2) COUP2 = MCOUP2 % + $ VAL(FLV_INDEX1) % P + DENOM = 1D0/(P3(0)**2-P3(1)**2-P3(2)**2-P3(3)**2 - M3 * (M3 -CI* + $ W3)) + V3%W(1)= DENOM*(-CI)*(COUP1*(F2 % W(1)*F1 % W(3)+F2 % W(2)*F1 % + $ W(4))+COUP2*(F2 % W(3)*F1 % W(1)+F2 % W(4)*F1 % W(2))) + V3%W(2)= DENOM*(-CI)*(COUP1*(F2 % W(2)*F1 % W(3)+F2 % W(1)*F1 % + $ W(4))-COUP2*(F2 % W(4)*F1 % W(1)+F2 % W(3)*F1 % W(2))) + V3%W(3)= DENOM*(-CI)*(COUP1*(+CI*(F2 % W(2)*F1 % W(3))-CI*(F2 % + $ W(1)*F1 % W(4)))+COUP2*(-CI*(F2 % W(4)*F1 % W(1))+CI*(F2 % W(3) + $ *F1 % W(2)))) + V3%W(4)= DENOM*(-CI)*(COUP1*(F2 % W(1)*F1 % W(3)-F2 % W(2)*F1 % + $ W(4))+COUP2*(-F2 % W(3)*F1 % W(1)+F2 % W(4)*F1 % W(2))) + FDJS1 = (FDN(0)*V3%W(1)-FDN(1)*V3%W(2)-FDN(2)*V3%W(3)-FDN(3) + $ *V3%W(4))/FDNQ + FDJS2 = (FDQ(0)*V3%W(1)-FDQ(1)*V3%W(2)-FDQ(2)*V3%W(3)-FDQ(3) + $ *V3%W(4)-DCONJG(FDQ(4))*V3%W(5))/FDNQ + V3%W(1) = V3%W(1)-FDQ(0)*FDJS1-FDN(0)*FDJS2 + V3%W(2) = V3%W(2)-FDQ(1)*FDJS1-FDN(1)*FDJS2 + V3%W(3) = V3%W(3)-FDQ(2)*FDJS1-FDN(2)*FDJS2 + V3%W(4) = V3%W(4)-FDQ(3)*FDJS1-FDN(3)*FDJS2 + V3%W(5) = V3%W(5)-FDQ(4)*FDJS1-FDN(4)*FDJS2 + END + + + diff --git a/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%VVV1_0.f b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%VVV1_0.f new file mode 100644 index 000000000..6969057bb --- /dev/null +++ b/tests/input_files/IOTestsComparison/IOTestFDGauge/FDgauge_standalone_fortran/%FD_fortran%Source%DHELAS%VVV1_0.f @@ -0,0 +1,202 @@ +C This File is Automatically generated by ALOHA +C The process calculated in this file is: +C P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - P(2,1)*Metric(1,3) + +C P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - P(1,3)*Metric(2,3) +C + SUBROUTINE VVV1_0(V1, V2, V3, COUP,VERTEX) + USE ALOHA_OBJECT + IMPLICIT NONE + COMPLEX*16 CI + PARAMETER (CI=(0D0,1D0)) + COMPLEX*16 CZERO + PARAMETER (CZERO=(0D0,0D0)) + COMPLEX*16 COUP + REAL*8 P1(0:3) + REAL*8 P2(0:3) + REAL*8 P3(0:3) + COMPLEX*16 TMP1 + COMPLEX*16 TMP2 + COMPLEX*16 TMP3 + COMPLEX*16 TMP4 + COMPLEX*16 TMP5 + COMPLEX*16 TMP6 + COMPLEX*16 TMP7 + COMPLEX*16 TMP8 + COMPLEX*16 TMP9 + TYPE(ALOHA) V1 + TYPE(ALOHA) V2 + TYPE(ALOHA) V3 + INTEGER*4 I + COMPLEX*16 VERTEX + P1(:) = V1 % P (:) + P2(:) = V2 % P (:) + P3(:) = V3 % P (:) + TMP1 = (V2 % W(1)*V1 % W(1)-V2 % W(2)*V1 % W(2)-V2 % W(3)*V1 % + $ W(3)-V2 % W(4)*V1 % W(4)) + TMP2 = (V3 % W(1)*P1(0)-V3 % W(2)*P1(1)-V3 % W(3)*P1(2)-V3 % W(4) + $ *P1(3)) + TMP3 = (V3 % W(1)*P2(0)-V3 % W(2)*P2(1)-V3 % W(3)*P2(2)-V3 % W(4) + $ *P2(3)) + TMP4 = (V2 % W(1)*P1(0)-V2 % W(2)*P1(1)-V2 % W(3)*P1(2)-V2 % W(4) + $ *P1(3)) + TMP5 = (V3 % W(1)*V1 % W(1)-V3 % W(2)*V1 % W(2)-V3 % W(3)*V1 % + $ W(3)-V3 % W(4)*V1 % W(4)) + TMP6 = (V2 % W(1)*P3(0)-V2 % W(2)*P3(1)-V2 % W(3)*P3(2)-V2 % W(4) + $ *P3(3)) + TMP7 = (V3 % W(1)*V2 % W(1)-V3 % W(2)*V2 % W(2)-V3 % W(3)*V2 % + $ W(3)-V3 % W(4)*V2 % W(4)) + TMP8 = (P2(0)*V1 % W(1)-P2(1)*V1 % W(2)-P2(2)*V1 % W(3)-P2(3)*V1 + $ % W(4)) + TMP9 = (V1 % W(1)*P3(0)-V1 % W(2)*P3(1)-V1 % W(3)*P3(2)-V1 % W(4) + $ *P3(3)) + VERTEX = COUP*(TMP1*(-CI*(TMP2)+CI*(TMP3))+(TMP5*(+CI*(TMP4)-CI + $ *(TMP6))+TMP7*(-CI*(TMP8)+CI*(TMP9)))) + END + + +C This File is Automatically generated by ALOHA +C The process calculated in this file is: +C Coup(1) * (P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - +C P(2,1)*Metric(1,3) + P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - +C P(1,3)*Metric(2,3)) +C Coup(2) * (Metric(1,2)) +C Coup(3) * (Metric(1,3)) +C Coup(4) * (P(1,2) - P(1,3)) +C + SUBROUTINE VVV1_VVS1_VSV2_VSS1_0(V1, V2, V3, COUP1, COUP2, COUP3 + $ , COUP4,VERTEX) + USE ALOHA_OBJECT + IMPLICIT NONE + COMPLEX*16 CI + PARAMETER (CI=(0D0,1D0)) + COMPLEX*16 CZERO + PARAMETER (CZERO=(0D0,0D0)) + COMPLEX*16 COUP1 + COMPLEX*16 COUP2 + COMPLEX*16 COUP3 + COMPLEX*16 COUP4 + REAL*8 P1(0:3) + REAL*8 P2(0:3) + REAL*8 P3(0:3) + COMPLEX*16 TMP1 + COMPLEX*16 TMP2 + COMPLEX*16 TMP3 + COMPLEX*16 TMP4 + COMPLEX*16 TMP5 + COMPLEX*16 TMP6 + COMPLEX*16 TMP7 + COMPLEX*16 TMP8 + COMPLEX*16 TMP9 + TYPE(ALOHA) V1 + TYPE(ALOHA) V2 + TYPE(ALOHA) V3 + INTEGER*4 I + COMPLEX*16 VERTEX + P1(:) = V1 % P (:) + P2(:) = V2 % P (:) + P3(:) = V3 % P (:) + VERTEX = (0D0,0D0) + TMP1 = (V2 % W(1)*V1 % W(1)-V2 % W(2)*V1 % W(2)-V2 % W(3)*V1 % + $ W(3)-V2 % W(4)*V1 % W(4)) + TMP2 = (V3 % W(1)*P1(0)-V3 % W(2)*P1(1)-V3 % W(3)*P1(2)-V3 % W(4) + $ *P1(3)) + TMP3 = (V3 % W(1)*P2(0)-V3 % W(2)*P2(1)-V3 % W(3)*P2(2)-V3 % W(4) + $ *P2(3)) + TMP4 = (V2 % W(1)*P1(0)-V2 % W(2)*P1(1)-V2 % W(3)*P1(2)-V2 % W(4) + $ *P1(3)) + TMP5 = (V3 % W(1)*V1 % W(1)-V3 % W(2)*V1 % W(2)-V3 % W(3)*V1 % + $ W(3)-V3 % W(4)*V1 % W(4)) + TMP6 = (V2 % W(1)*P3(0)-V2 % W(2)*P3(1)-V2 % W(3)*P3(2)-V2 % W(4) + $ *P3(3)) + TMP7 = (V3 % W(1)*V2 % W(1)-V3 % W(2)*V2 % W(2)-V3 % W(3)*V2 % + $ W(3)-V3 % W(4)*V2 % W(4)) + TMP8 = (P2(0)*V1 % W(1)-P2(1)*V1 % W(2)-P2(2)*V1 % W(3)-P2(3)*V1 + $ % W(4)) + TMP9 = (V1 % W(1)*P3(0)-V1 % W(2)*P3(1)-V1 % W(3)*P3(2)-V1 % W(4) + $ *P3(3)) + VERTEX = VERTEX + COUP1*(TMP1*(-CI*(TMP2)+CI*(TMP3))+(TMP5*(+CI + $ *(TMP4)-CI*(TMP6))+TMP7*(-CI*(TMP8)+CI*(TMP9)))) + VERTEX = VERTEX + COUP2*(-CI * TMP1*V3 % W(5)) + VERTEX = VERTEX + COUP3*(-CI * TMP5*V2 % W(5)) + VERTEX = VERTEX + COUP4*V2 % W(5)*V3 % W(5)*(-CI*(TMP8)+CI*(TMP9) + $ ) + END + + +C This File is Automatically generated by ALOHA +C The process calculated in this file is: +C Coup(1) * (P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - +C P(2,1)*Metric(1,3) + P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - +C P(1,3)*Metric(2,3)) +C Coup(2) * (Metric(1,3)) +C Coup(3) * (P(1,3) - P(1,2)) +C Coup(4) * (Metric(2,3)) +C Coup(5) * (P(2,3) - P(2,1)) +C Coup(6) * (P(3,1) - P(3,2)) +C + SUBROUTINE VVV1_VSV2_VSS2_SVV2_SVS2_SSV3_0(V1, V2, V3, COUP1, + $ COUP2, COUP3, COUP4, COUP5, COUP6,VERTEX) + USE ALOHA_OBJECT + IMPLICIT NONE + COMPLEX*16 CI + PARAMETER (CI=(0D0,1D0)) + COMPLEX*16 CZERO + PARAMETER (CZERO=(0D0,0D0)) + COMPLEX*16 COUP1 + COMPLEX*16 COUP2 + COMPLEX*16 COUP3 + COMPLEX*16 COUP4 + COMPLEX*16 COUP5 + COMPLEX*16 COUP6 + REAL*8 P1(0:3) + REAL*8 P2(0:3) + REAL*8 P3(0:3) + COMPLEX*16 TMP1 + COMPLEX*16 TMP2 + COMPLEX*16 TMP3 + COMPLEX*16 TMP4 + COMPLEX*16 TMP5 + COMPLEX*16 TMP6 + COMPLEX*16 TMP7 + COMPLEX*16 TMP8 + COMPLEX*16 TMP9 + TYPE(ALOHA) V1 + TYPE(ALOHA) V2 + TYPE(ALOHA) V3 + INTEGER*4 I + COMPLEX*16 VERTEX + P1(:) = V1 % P (:) + P2(:) = V2 % P (:) + P3(:) = V3 % P (:) + VERTEX = (0D0,0D0) + TMP1 = (V2 % W(1)*V1 % W(1)-V2 % W(2)*V1 % W(2)-V2 % W(3)*V1 % + $ W(3)-V2 % W(4)*V1 % W(4)) + TMP2 = (V3 % W(1)*P1(0)-V3 % W(2)*P1(1)-V3 % W(3)*P1(2)-V3 % W(4) + $ *P1(3)) + TMP3 = (V3 % W(1)*P2(0)-V3 % W(2)*P2(1)-V3 % W(3)*P2(2)-V3 % W(4) + $ *P2(3)) + TMP4 = (V2 % W(1)*P1(0)-V2 % W(2)*P1(1)-V2 % W(3)*P1(2)-V2 % W(4) + $ *P1(3)) + TMP5 = (V3 % W(1)*V1 % W(1)-V3 % W(2)*V1 % W(2)-V3 % W(3)*V1 % + $ W(3)-V3 % W(4)*V1 % W(4)) + TMP6 = (V2 % W(1)*P3(0)-V2 % W(2)*P3(1)-V2 % W(3)*P3(2)-V2 % W(4) + $ *P3(3)) + TMP7 = (V3 % W(1)*V2 % W(1)-V3 % W(2)*V2 % W(2)-V3 % W(3)*V2 % + $ W(3)-V3 % W(4)*V2 % W(4)) + TMP8 = (P2(0)*V1 % W(1)-P2(1)*V1 % W(2)-P2(2)*V1 % W(3)-P2(3)*V1 + $ % W(4)) + TMP9 = (V1 % W(1)*P3(0)-V1 % W(2)*P3(1)-V1 % W(3)*P3(2)-V1 % W(4) + $ *P3(3)) + VERTEX = VERTEX + COUP1*(TMP1*(-CI*(TMP2)+CI*(TMP3))+(TMP5*(+CI + $ *(TMP4)-CI*(TMP6))+TMP7*(-CI*(TMP8)+CI*(TMP9)))) + VERTEX = VERTEX + COUP2*(-CI * TMP5*V2 % W(5)) + VERTEX = VERTEX + COUP3*V2 % W(5)*V3 % W(5)*(-CI*(TMP9)+CI*(TMP8) + $ ) + VERTEX = VERTEX + COUP4*(-CI * TMP7*V1 % W(5)) + VERTEX = VERTEX + COUP5*V1 % W(5)*V3 % W(5)*(-CI*(TMP6)+CI*(TMP4) + $ ) + VERTEX = VERTEX + COUP6*V1 % W(5)*V2 % W(5)*(-CI*(TMP2)+CI*(TMP3) + $ ) + END + + diff --git a/tests/parallel_tests/test_aloha.py b/tests/parallel_tests/test_aloha.py index 7a6c9f1a7..96ace7dce 100755 --- a/tests/parallel_tests/test_aloha.py +++ b/tests/parallel_tests/test_aloha.py @@ -3374,16 +3374,19 @@ def test_short_multiple_lorentz_subset(self): TYPE(ALOHA) F1 TYPE(ALOHA) F2 TYPE(ALOHA) V3 -COMPLEX*16 TMP COMPLEX*16 VERTEX -CALL FFV1C1_0(F2,F1,V3,COUP1,VERTEX) -CALL FFV2C1_0(F2,F1,V3,COUP2,TMP) -VERTEX = VERTEX + TMP END""" - + data = [ l.strip() for l in fsock.read().split('\n')] for line in goal.split('\n'): self.assertIn(line.strip(), data) + # the combined routine is a merged expression, not a wrapper + # calling the routine of each single Lorentz structure: no call, + # and the two couplings enter the same expression. The names of + # the temporaries depend on what the kernel computed before, so + # they are not part of the comparison. + self.assertFalse([l for l in data if l.startswith('CALL ')]) + self.assertTrue([l for l in data if 'COUP1' in l and 'COUP2' in l])