Merge pull request #3 from zeniheisser/gpucpp

merging upstream gpucpp changes into rex

.gitignore

@@ -11,6 +11,7 @@ ptest*
 *.pkl
 p_ME_*
 *.log
+.DS_Store
 input/default_run_card_lo.dat
 input/default_run_card_nlo.dat
 input/mg5_configuration.txt
@@ -71,4 +72,5 @@ tmp
 check_sa_born_splitOrders.f
 models/*/
 models/*.tgz
+tests/status
 nsqso_born.inc

Template/LO/Source/dsample.f

@@ -37,8 +37,12 @@ subroutine sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs,VECSIZE_USED)
       integer jmax,i,j,ipole
       integer itmax_adjust
 
-      integer imirror, iproc, iconf
-      integer ivec ! position of the event in the vector (max is VECSIZE_MEMMAX, loops go over VECSIZE_USED)
+c      integer imirror, iproc, iconf
+      integer imirror_vec(NB_WARP), iproc, ICONF_VEC(NB_WARP) 
+      integer ivec              ! position of the event in the vector (max is VECSIZE_MEMMAX, loops go over VECSIZE_USED)
+      integer ilock             !  position of the event in the current warp (max is WARP_SIZE)
+      integer iwarp               ! position of the current warp (max is NB_WARP)
+c     NOTE THAT IVEC = (IWARP-1)*NB_WARP + ILOCK      
 
 c
 c     External
@@ -161,6 +165,8 @@ subroutine sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs,VECSIZE_USED)
       ievent = 0
       iter = 1
       ivec = 0
+      ilock = 0
+      iwarp = 1
       do while(iter .le. itmax)
 c
 c     Get integration point
@@ -174,6 +180,11 @@ subroutine sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs,VECSIZE_USED)
             CUTSPASSED=.FALSE.
             if (passcuts(p,VECSIZE_USED)) then
                ivec=ivec+1
+               ilock = ilock+1
+               if (ilock.gt.WARP_SIZE)then
+                  ilock = 1
+                  iwarp = iwarp +1
+               endif
 c              write(*,*) 'pass_point ivec is ', ivec
                all_p(:,ivec) = p(:)
                all_wgt(ivec) = wgt
@@ -187,25 +198,35 @@ subroutine sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs,VECSIZE_USED)
 c               bckp(i) = fx
 c               write(*,*) i, all_wgt(i), fx, all_wgt(i)*fx
 c               all_wgt(i) = all_wgt(i)*fx
-               if (ivec.lt.VECSIZE_USED)then
+               if (ilock.ne.WARP_SIZE)then
                   cycle
                endif
-               ivec=0
+
                if (VECSIZE_USED.le.1) then
                   all_fx(1) = dsig(all_p, all_wgt,0)
                else
-               do i=1, VECSIZE_USED
+c                 Here "i" is the position in the full grid of the event                  
+                  do i=(iwarp-1)*WARP_SIZE+1, iwarp*warp_size
+
 c                 need to restore common block                  
                   xbk(:) = all_xbk(:, i)
                   cm_rap = all_cm_rap(i)
                   q2fact(:) = all_q2fact(:,i)
                   CUTSDONE=.TRUE.
                   CUTSPASSED=.TRUE.
-                  call prepare_grouping_choice(all_p(1,i), all_wgt(i), i.eq.1)
+                  call prepare_grouping_choice(all_p(1,i), all_wgt(i),i.eq.(iwarp-1)*WARP_SIZE+1)
                enddo
-               call select_grouping(imirror, iproc, iconf, all_wgt, VECSIZE_USED)
+               call select_grouping(imirror_vec(iwarp), iproc, iconf_vec(iwarp), all_wgt, iwarp)
+               if (ivec.lt.VECSIZE_USED)then
+                  cycle
+               endif
+c              reset variable for the next grid               
+               ivec = 0
+               ilock = 0
+               iwarp =1
+
                call dsig_vec(all_p, all_wgt, all_xbk, all_q2fact, all_cm_rap,
-     &                          iconf, iproc, imirror, all_fx,VECSIZE_USED)
+     &                          iconf_vec, iproc, imirror_vec, all_fx,VECSIZE_USED)
 
                 do i=1, VECSIZE_USED
 c                 need to restore common block                  

Template/LO/SubProcesses/reweight.f

@@ -620,6 +620,7 @@ logical function setclscales(p, keepq2bck, ivec)
 
 c     Variables for keeping track of jets
       logical goodjet(n_max_cl)
+c111      logical set_goodjet(n_max_cl) ! variable for debugging unset call 
       integer fsnum(2),ida(2),imo,jcode
       logical chclusold,fail,increasecode
       save chclusold
@@ -633,6 +634,13 @@ logical function setclscales(p, keepq2bck, ivec)
       external is_octet
       setclscales=.true.
 
+c    WARNING: goodjet() is sometimes accessed for variable which are not
+c      initialised by the algorithm (the line where this was happening is
+c      flagged). To be on a safe side, we initialise all entry to False.
+c      In case of weird behavior, uncomment line c111 to debug/investigate 
+      goodjet(:) = .false.
+c111     set_goodjet(:) = .false.
+
       if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.q2fact(2).gt.0.and.scale.gt.0) then
          if(use_syst)then
             s_scale(ivec)=scale
@@ -737,11 +745,13 @@ logical function setclscales(p, keepq2bck, ivec)
 c        goodjet gives whether this cluster line is considered a jet
 c        i.e. if all related/previous clustering are jet
          goodjet(ibeam(i))=partonline(i)
+c111         set_goodjet(ibeam(i))= .true.
       enddo
 
       do i=3,nexternal
          j=ishft(1,i-1)
          goodjet(j)=isjet(ipdgcl(j,igraphs(1),iproc))
+c111         set_goodjet(j)= .true.
       enddo
 
 c     Go through clusterings and set factorization scale points for use in dsig
@@ -782,8 +792,10 @@ logical function setclscales(p, keepq2bck, ivec)
               else if (jfirst(j).eq.0) then
                  jfirst(j) = n
                  goodjet(imo)=.false.
+c111                 set_goodjet(imo)= .true.
               else
                  goodjet(imo)=.false.
+c111                 set_goodjet(imo)= .true.
               endif
 c             If not jet vertex, increase jcode. This is needed
 c             e.g. in VBF if we pass over to the other side and hit
@@ -799,6 +811,7 @@ logical function setclscales(p, keepq2bck, ivec)
                endif
 c             Consider t-channel jet radiations as jets only if
 c             FS line is a jet line
+c111               if(.not.set_goodjet(ida(3-i))) stop 1 
               if(goodjet(ida(3-i))) then
                  if(partonline(j).or.
      $ ipdgcl(ida(3-i),igraphs(1),iproc).eq.21)then
@@ -862,6 +875,7 @@ logical function setclscales(p, keepq2bck, ivec)
            endif
 c          Set goodjet to false for mother
               goodjet(imocl(n))=.false.
+c111              set_goodjet(imocl(n))= .true.
               cycle
            endif
 
@@ -881,6 +895,7 @@ logical function setclscales(p, keepq2bck, ivec)
            goodjet(imocl(n))=
      $          (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
      $          goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
+c111           set_goodjet(imocl(n))= .true.
 
 c       check case with g > g g
 c       where the hardest gluon is not goodjet but the other is.
@@ -907,6 +922,7 @@ logical function setclscales(p, keepq2bck, ivec)
 c     $             ipart(1,idacl(n,2)),'/', ipart(2,idacl(n,2))
 c              write(*,*) 'googjet', goodjet(idacl(n,1)),goodjet(idacl(n,2))
               if (ipart(1,imocl(n)).eq.ipart(1, idacl(n,1))) then
+c111                 if(.not.set_goodjet(idacl(n,1)).or..not.set_goodjet(idacl(n,2))) stop 1                  
                   if (.not.goodjet(idacl(n,1)).and.goodjet(idacl(n,2))) then
 c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,2)), ipart(2,idacl(n,2))
                      do m =n_max_cl,n,-1
@@ -917,6 +933,7 @@ logical function setclscales(p, keepq2bck, ivec)
                      enddo
                   endif
                else
+c111                 if(.not.set_goodjet(idacl(n,1)).or..not.set_goodjet(idacl(n,2))) stop 1
                   if (.not.goodjet(idacl(n,2)).and.goodjet(idacl(n,1))) then
 c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,1)), ipart(2,idacl(n,1))
                      do m =n_max_cl,n,-1
@@ -1227,8 +1244,12 @@ logical function setclscales(p, keepq2bck, ivec)
          do i=1,2
             do j=1,2
 c              First adjust goodjet based on iqjets
-               if(goodjet(ida(i)).and.ipart(j,ida(i)).gt.2)then
-                  if(iqjets(ipart(j,ida(i))).eq.0) goodjet(ida(i))=.false.
+c111              if(.not.set_goodjet(ida(i)).and.ipart(j,ida(i)).gt.2) stop 1            
+               if(goodjet(ida(i)).and.ipart(j,ida(i)).gt.2)then ! This is where  goodjet can be used uninitialised #111
+                  if(iqjets(ipart(j,ida(i))).eq.0)then
+                      goodjet(ida(i))=.false.
+c111                      set_goodjet(ida(i))=.true.
+                  endif
                endif
 c              Now reset ptclus if jet vertex
                if(ipart(j,ida(i)).gt.2) then

Template/LO/bin/internal/restore_data

@@ -45,30 +45,24 @@ for i in `cat subproc.mg` ; do
     else
 	cp  results.dat $1_results.dat  >& /dev/null
     fi
-    for k in G* ; do
-	if [[ ! -d $k ]]; then
-	    continue
-	fi
-	cd $k
-	for j in $1_results.dat ; do
-	    if [[ -e $j ]] ; then
-		cp  $j results.dat
-	    else
-		cp  results.dat $j
-	    fi
-	done
-	for j in $1_ftn26.gz ; do
-	    if [[ -e $j ]]; then
-		rm -f ftn26 >& /dev/null
-		rm -f $1_ftn26 >& /dev/null
-		gunzip $j
-		cp $1_ftn26 ftn26
-		gzip $1_ftn26
-	    fi
-	done
-	cd ../
-    done
     cd ../
 done
 
-
+find . -mindepth 2 -maxdepth 2 -type d -name 'G*' -print0 \
+    | xargs --null -P "$(nproc --all)" -I{} bash -c "
+cd {}
+for j in $1_results.dat ; do
+    if [[ -e \$j ]] ; then
+        cp \$j results.dat
+    else
+        cp results.dat \$j
+    fi
+done
+for j in $1_ftn26.gz ; do
+    if [[ -e \$j ]]; then
+        rm -f ftn26 >& /dev/null
+        rm -f $1_ftn26 >& /dev/null
+        gunzip --keep \$j
+        mv $1_ftn26 ftn26
+    fi
+done"

aloha/template_files/aloha_functions.f

@@ -1201,7 +1201,7 @@ subroutine rotxxx(p,q , prot)
 c
       implicit none
       double precision p(0:3),q(0:3),prot(0:3),qt2,qt,psgn,qq,p1
-
+      volatile qt, p1, qq ! prevent optimizations with -O3 (workaround for SIGFPE crashes in rotxxx: madgraph5/madgraph4gpu#855)
       double precision rZero, rOne
       parameter( rZero = 0.0d0, rOne = 1.0d0 )
 

madgraph/core/base_objects.py

@@ -946,6 +946,7 @@ def get_canonical_couplings_keys_order(self):
         '%s_%s_%s'%(self['color'][k[0]],self['lorentz'][k[1]],self['couplings'][k]))
 
 
+
 #===============================================================================
 # InteractionList
 #===============================================================================

madgraph/interface/madevent_interface.py

@@ -3967,6 +3967,9 @@ def do_create_gridpack(self, line):
         for P in Pdir: 
             allG = misc.glob('G*', path=P)
             for G in allG:
+                # avoid case where some file starts with G (madgraph5/madgraph4gpu#947)
+                if not os.path.isdir(G): 
+                    continue
                 if pjoin(P, G) not in Gdir:
                     logger.debug('removing %s', pjoin(P,G))
                     shutil.rmtree(pjoin(P,G))

madgraph/interface/madgraph_interface.py

@@ -1534,7 +1534,7 @@ def check_set(self, args, log=True):
 
         if args[0] in ['group_subprocesses']:
             if args[1].lower() not in ['false', 'true', 'auto', 'gpu']:
-                raise self.InvalidCmd('%s needs argument False, True or Auto, got %s' % \
+                raise self.InvalidCmd('%s needs argument False, True, gpu or Auto, got %s' % \
                                       (args[0], args[1]))
         if args[0] in ['ignore_six_quark_processes']:
             if args[1] not in list(self._multiparticles.keys()) and args[1].lower() != 'false':
@@ -1636,7 +1636,6 @@ def check_output(self, args, default='madevent'):
                 args.pop(0)
                 if hasattr(output_cls, 'change_output_args'):
                     args[:] = output_cls.change_output_args(args, self) 
-                self.export
             else:
                 self._export_format = default
         else:
@@ -2499,7 +2498,7 @@ def complete_open(self, text, line, begidx, endidx):
     def complete_output(self, text, line, begidx, endidx,
                         possible_options = ['f', 'noclean', 'nojpeg'],
                         possible_options_full = ['-f', '-noclean', '-nojpeg', '--noeps=True','--hel_recycling=False',
-                                                 '--jamp_optim=', '--t_strategy=', '--vector_size=4']):
+                                                 '--jamp_optim=', '--t_strategy=', '--vector_size=4', '--nb_wrap=1']):
         "Complete the output command"
 
         possible_format = self._export_formats
@@ -2597,11 +2596,13 @@ def complete_set(self, text, line, begidx, endidx):
             return self.list_completion(text, opts)
 
         if len(args) == 2:
-            if args[1] in ['group_subprocesses', 'complex_mass_scheme',\
+            if args[1] in ['complex_mass_scheme',\
                            'loop_optimized_output', 'loop_color_flows',\
                            'include_lepton_initiated_processes',\
                            'low_mem_multicore_nlo_generation', 'nlo_mixed_expansion']:
                 return self.list_completion(text, ['False', 'True', 'default'])
+            elif args[1] in ['group_subprocesses']:
+                return self.list_completion(text, ['False', 'True', 'Auto', 'gpu', 'nlo'])
             elif args[1] in ['ignore_six_quark_processes']:
                 return self.list_completion(text, list(self._multiparticles.keys()))
             elif args[1].lower() == 'ewscheme':
@@ -6147,7 +6148,6 @@ def advanced_install(self, tool_to_install,
             return_code = misc.call([sys.executable, pjoin(MG5DIR,'HEPTools',
               'HEPToolsInstallers', 'HEPToolInstaller.py'), tool,'--prefix=%s'%
               prefix] + compiler_options + add_options)
-        misc.sprint(return_code)
         if return_code in [0,11]:
             logger.info("%s successfully installed in %s."%(
                    tool_to_install, prefix),'$MG:color:GREEN')
@@ -8127,6 +8127,7 @@ def do_output(self, line):
 
         args = self.split_arg(line)
         # Check Argument validity
+        self._export_plugin = None
         self.check_output(args)
 
         noclean = '-noclean' in args

madgraph/iolibs/export_cpp.py

@@ -1517,8 +1517,8 @@ def get_sigmaKin_lines(self, color_amplitudes, write=True):
             """
             replace_dict['madE_caclwfcts_call'] = '&multi_chanel_num, &multi_chanel_denom'
             replace_dict['madE_update_answer'] = '   allMEs[iproc*nprocesses + ievt] *= multi_chanel_num/multi_chanel_denom;'
+
             multi_channel = self.get_multi_channel_dictionary(self.matrix_elements[0].get('diagrams'), self.include_multi_channel)
-            replace_dict['is_LC'] = self.get_icolamp_lines(multi_channel, self.matrix_elements[0], 1)
             replace_dict['nb_channel'] = len(multi_channel)
             replace_dict['nb_color'] = max(1, len(self.matrix_elements[0].get('color_basis')))
 
@@ -1760,7 +1760,6 @@ def write_process_cc_file(self, writer):
             return replace_dict
 
 
-
     def get_icolamp_lines(self, mapconfigs, matrix_element, num_matrix_element):
         """Return the ICOLAMP matrix, showing which JAMPs contribute to
         which configs (diagrams)."""