Cheaper colour sum: DDM basis, Kleiss-Kuijf flows, reflection folding, BLAS, and the same for madmatrix

Template/LO/Source/run.inc

@@ -113,3 +113,8 @@ c                          2 means approximation by the denominator of the propa
 c
        double precision limhel
        common/to_limhel/limhel
+c
+c      whether the color sum may go through the batched BLAS-3 call, where
+c      the matrix element was generated with one
+       logical blas_color_sum
+       common/to_blas_color_sum/blas_color_sum

Template/LO/SubProcesses/makefile

@@ -26,6 +26,11 @@ endif
 
 LINKLIBS = $(LINK_MADLOOP_LIB) $(LINK_LOOP_LIBS) -L../../lib/ -ldhelas -ldsample -lmodel -lgeneric -lpdf -lgammaUPC -lcernlib $(llhapdf) -lbias 
 
+# what the batched color sum is linked against, empty unless a matrix
+# element was generated with one
+BLASLIBS =
+LINKLIBS += $(BLASLIBS)
+
 LIBS = $(LIBDIR)libbias.$(libext) $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libdsample.$(libext) $(LIBDIR)libgeneric.$(libext) $(LIBDIR)libpdf.$(libext) $(LIBDIR)libgammaUPC.$(libext) $(LIBDIR)libmodel.$(libext) $(LIBDIR)libcernlib.$(libext) $(MADLOOP_LIB) $(LOOP_LIBS)
 
 ifneq ("$(wildcard ../../Source/RUNNING)","")

madgraph/core/color_algebra.py

@@ -332,6 +332,13 @@ def simplify(self):
 
         return ColorFactor([col_str1, col_str2])
 
+    def complex_conjugate(self):
+        """f (and d) are real, so complex conjugation leaves them untouched.
+        The default behaviour of reversing the indices would introduce a
+        spurious sign, since f is totally antisymmetric."""
+
+        return self
+
 #===============================================================================
 # d
 #===============================================================================
@@ -1128,18 +1135,39 @@ def extend_str(self, new_col_fact):
         for col_str in new_col_fact:
             self.append_str(col_str)
 
+    @staticmethod
+    def similarity_key(col_str):
+        """Hashable key which is equal for two color strings exactly when
+        ColorString.is_similar says they are, i.e. same Nc power, same
+        imaginary character and same canonical representation (both the
+        canonical structure and the index replacement dictionary)."""
+
+        canonical, repl_dict = col_str.to_canonical()
+        return (col_str.Nc_power, col_str.is_imaginary, canonical,
+                tuple(sorted(repl_dict.items())))
+
     def simplify(self):
         """Returns a new color factor where each color string has been
-        simplified once and similar strings have been added."""
+        simplified once and similar strings have been added.
+
+        Similar strings are looked up through a dictionary rather than by
+        scanning the strings accumulated so far: the color factors appearing
+        for high multiplicity processes have thousands of terms, and the linear
+        scan of append_str made this quadratic. Insertion order is preserved,
+        so the resulting color factor is identical to the one a scan produces."""
 
         new_col_factor = ColorFactor()
+        similar = {}
         # Simplify
         for col_str in self:
             res = col_str.simplify()
-            if res:
-                new_col_factor.extend_str(res)
-            else:
-                new_col_factor.append_str(col_str)
+            for new_str in (res if res else [col_str]):
+                key = self.similarity_key(new_str)
+                try:
+                    similar[key].add(new_str)
+                except KeyError:
+                    similar[key] = new_str
+                    new_col_factor.append(new_str)
 
         # Only returns non zero elements
         return ColorFactor([col_str for col_str in \