Index: Splutter/GPU/SplutterGPU.py
===================================================================
--- Splutter/GPU/SplutterGPU.py	(revision 67)
+++ Splutter/GPU/SplutterGPU.py	(revision 68)
@@ -25,6 +25,20 @@
 import math
 from socket import gethostname
 
+# find prime factors of a number
+# Get for WWW :
+# http://pythonism.wordpress.com/2008/05/17/looking-at-factorisation-in-python/
+def PrimeFactors(x):
+  factorlist=numpy.array([]).astype('uint32')
+  loop=2
+  while loop<=x:
+    if x%loop==0:
+      x/=loop
+      factorlist=numpy.append(factorlist,[loop])
+    else:
+      loop+=1
+  return factorlist
+  
 # Try to find the best thread number in Hybrid approach (Blocks&Threads)
 # output is thread number
 def BestThreadsNumber(jobs):
@@ -53,11 +67,6 @@
 def Mylq2(N, T1,s,c1,c2,p):
   return (T1*(s+p/N)+c1*N+c2*N*N)
 
-prout="""
-
-"""
-
-
 KERNEL_CODE_CUDA="""
 
 // Marsaglia RNG very simple implementation
@@ -275,8 +284,6 @@
 __kernel void SplutterGlobal(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
    __private const ulong id=(ulong)get_global_id(0);
-   __private const ulong size=(ulong)get_global_size(0);
-   __private const ulong block=(ulong)space/(ulong)size;
    
    __private uint z=seed_z-(uint)id;
    __private uint w=seed_w+(uint)id;
@@ -287,7 +294,7 @@
    for (__private ulong i=0;i<iterations;i++) {
 
       // Dense version
-      __private size_t position=(size_t)( ((ulong)MWC*(ulong)space)/MAX );
+      __private size_t position=(size_t)( MWC%space );
 
       atomic_inc(&s[position]);
    }
@@ -296,59 +303,9 @@
 
 }
 
-__kernel void SplutterGlobalDense(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+__kernel void SplutterLocal(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   __private const ulong id=(ulong)get_global_id(0);
-   __private const ulong size=(ulong)get_global_size(0);
-   __private const ulong block=(ulong)space/(ulong)size;
-   
-   __private uint z=seed_z-(uint)id;
-   __private uint w=seed_w+(uint)id;
-
-   __private uint jsr=seed_z;
-   __private uint jcong=seed_w;
-
-   for (__private ulong i=0;i<iterations;i++) {
-
-      // Dense version
-      __private size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
-
-      s[position]++;
-   }
-
-   barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
-
-}
-
-__kernel void SplutterGlobalSparse(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
-{
-   __private const ulong id=(ulong)get_global_id(0);
-   __private const ulong size=(ulong)get_global_size(0);
-   __private const ulong block=(ulong)space/(ulong)size;
-   
-   __private uint z=seed_z-(uint)id;
-   __private uint w=seed_w+(uint)id;
-
-   __private uint jsr=seed_z;
-   __private uint jcong=seed_w;
-
-   for (__private ulong i=0;i<iterations;i++) {
-
-      // Sparse version
-      __private size_t position=(size_t)( (ulong)MWC*block/MAX*size+id );
-
-      s[position]++;
-   }
-
-   barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
-
-}
-
-__kernel void SplutterLocalDense(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
-{
    __private const ulong id=(ulong)get_local_id(0);
-   __private const ulong size=(ulong)get_local_size(0);
-   __private const ulong block=(ulong)space/(ulong)size;
    
    __private uint z=seed_z-(uint)id;
    __private uint w=seed_w+(uint)id;
@@ -359,20 +316,19 @@
    for (__private ulong i=0;i<iterations;i++) {
 
       // Dense version
-      __private size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
+      //__private size_t position=(size_t)( (MWC+id*block)%space );
+      __private size_t position=(size_t)( MWC%space );
 
-      s[position]++;
+      atomic_inc(&s[position]);
    }
 
    barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
 
 }
 
-__kernel void SplutterLocalSparse(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+__kernel void SplutterHybrid(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   __private const ulong id=(ulong)get_local_id(0);
-   __private const ulong size=(ulong)get_local_size(0);
-   __private const ulong block=(ulong)space/(ulong)size;
+   __private const ulong id=(ulong)(get_global_id(0)+get_local_id(0));
    
    __private uint z=seed_z-(uint)id;
    __private uint w=seed_w+(uint)id;
@@ -382,63 +338,17 @@
 
    for (__private ulong i=0;i<iterations;i++) {
 
-      // Sparse version
-      __private size_t position=(size_t)( (ulong)MWC*block/MAX*size+id );
-
-      s[position]++;
-   }
-
-   barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
-
-}
-
-__kernel void SplutterHybridDense(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
-{
-   __private const ulong id=(ulong)(get_global_id(0));
-   __private const ulong size=(ulong)(get_local_size(0)*get_num_groups(0));
-   __private const ulong block=(ulong)space/(ulong)size;
-   
-   __private uint z=seed_z-(uint)id;
-   __private uint w=seed_w+(uint)id;
-
-   __private uint jsr=seed_z;
-   __private uint jcong=seed_w;
-
-   for (__private ulong i=0;i<iterations;i++) {
-
       // Dense version
-      __private size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
+      __private size_t position=(size_t)( MWC%space );
 
-      s[position]++;
+      atomic_inc(&s[position]);
    }
       
 }
 
-__kernel void SplutterHybridSparse(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
-{
-   __private const ulong id=(ulong)(get_global_id(0));
-   __private const ulong size=(ulong)(get_local_size(0)*get_num_groups(0));
-   __private const ulong block=(ulong)space/(ulong)size;
-   
-   __private uint z=seed_z-(uint)id;
-   __private uint w=seed_w+(uint)id;
-
-   __private uint jsr=seed_z;
-   __private uint jcong=seed_w;
-
-   for (__private ulong i=0;i<iterations;i++) {
-
-      // Sparse version
-      __private size_t position=(size_t)( (ulong)MWC*block/MAX*size+id );
-
-      s[position]++;
-   }
-      
-}
-
 """
 
-def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,Density):
+def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,Density,Memory):
 
   # Avec PyCUDA autoinit, rien a faire !
 
@@ -536,7 +446,7 @@
 
 
 def MetropolisOpenCL(circle,iterations,steps,jobs,
-                     ParaStyle,Alu,Device,Density):
+                     ParaStyle,Alu,Device,Memory):
 	
   # Initialisation des variables en les CASTant correctement
 
@@ -616,7 +526,8 @@
   print MaxWorks,2**(int)(numpy.log2(MemoryXPU))
   
   #Splutter=numpy.zeros((MaxWorks/jobs)*jobs).astype(numpy.uint32)
-  Splutter=numpy.zeros(jobs*16).astype(numpy.uint32)
+  #Splutter=numpy.zeros(jobs*16).astype(numpy.uint32)
+  Splutter=numpy.zeros(Memory).astype(numpy.uint32)
 
   for i in range(steps):
 		
@@ -645,67 +556,37 @@
       #                                      numpy.uint64(iterationsCL),
       #                                      numpy.uint32(nprnd(2**30/jobs)),
       #                                      numpy.uint32(nprnd(2**30/jobs)))
-      if Density=='Dense':
-        CLLaunch=MetropolisCL.SplutterGlobalDense(queue,(jobs,),None,
-                                                  SplutterCL,
-                                                  numpy.uint32(len(Splutter)),
-                                                  numpy.uint64(iterationsCL),
-                                                  numpy.uint32(521288629),
-                                                  numpy.uint32(362436069))
-      elif Density=='Sparse':
-        CLLaunch=MetropolisCL.SplutterGlobalSparse(queue,(jobs,),None,
-                                                   SplutterCL,
-                                                   numpy.uint32(len(Splutter)),
-                                                   numpy.uint64(iterationsCL),
-                                                   numpy.uint32(521288629),
-                                                   numpy.uint32(362436069))
+      CLLaunch=MetropolisCL.SplutterGlobal(queue,(jobs,),None,
+                                           SplutterCL,
+                                           numpy.uint32(len(Splutter)),
+                                           numpy.uint64(iterationsCL),
+                                           numpy.uint32(nprnd(2**30/jobs)),
+                                           numpy.uint32(nprnd(2**30/jobs)))
         
-      else:
-        CLLaunch=MetropolisCL.SplutterGlobal(queue,(jobs,),None,
-                                             SplutterCL,
-                                             numpy.uint32(len(Splutter)),
-                                             numpy.uint64(iterationsCL),
-                                             numpy.uint32(521288629),
-                                             numpy.uint32(362436069))
-        
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs,1,ParaStyle)
     elif ParaStyle=='Hybrid':
-      threads=BestThreadsNumber(jobs)
+      #threads=BestThreadsNumber(jobs)
+      threads=BestThreadsNumber(256)
+      print "print",threads      
       # en OpenCL, necessaire de mettre un Global_id identique au local_id
-      if Density=='Dense':
-        CLLaunch=MetropolisCL.SplutterHybridDense(queue,(jobs,),(threads,),
-                                                  SplutterCL,
-                                                  numpy.uint32(len(Splutter)),
-                                                  numpy.uint64(iterationsCL),
-                                                  numpy.uint32(nprnd(2**30/jobs)),
-                                                  numpy.uint32(nprnd(2**30/jobs)))
-      elif Density=='Sparse':
-        CLLaunch=MetropolisCL.SplutterHybridSparse(queue,(jobs,),(threads,),
-                                                   SplutterCL,
-                                                   numpy.uint32(len(Splutter)),
-                                                   numpy.uint64(iterationsCL),
-                                                   numpy.uint32(nprnd(2**30/jobs)),
-                                                   numpy.uint32(nprnd(2**30/jobs)))
+      CLLaunch=MetropolisCL.SplutterHybrid(queue,(jobs,),(threads,),
+                                           SplutterCL,
+                                           numpy.uint32(len(Splutter)),
+                                           numpy.uint64(iterationsCL),
+                                           numpy.uint32(nprnd(2**30/jobs)),
+                                           numpy.uint32(nprnd(2**30/jobs)))
         
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs/threads,threads,ParaStyle)
     else:
       # en OpenCL, necessaire de mettre un global_id identique au local_id
-      if Density=='Dense':
-        CLLaunch=MetropolisCL.SplutterLocalDense(queue,(jobs,),(jobs,),
-                                                 SplutterCL,
-                                                 numpy.uint32(len(Splutter)),
-                                                 numpy.uint64(iterationsCL),
-                                                 numpy.uint32(nprnd(2**30/jobs)),
-                                                 numpy.uint32(nprnd(2**30/jobs)))
-      elif Density=='Sparse':
-        CLLaunch=MetropolisCL.SplutterLocalSparse(queue,(jobs,),(jobs,),
-                                                  SplutterCL,
-                                                  numpy.uint32(len(Splutter)),
-                                                  numpy.uint64(iterationsCL),
-                                                  numpy.uint32(nprnd(2**30/jobs)),
-                                                  numpy.uint32(nprnd(2**30/jobs)))
+      CLLaunch=MetropolisCL.SplutterLocal(queue,(jobs,),(jobs,),
+                                          SplutterCL,
+                                          numpy.uint32(len(Splutter)),
+                                          numpy.uint64(iterationsCL),
+                                          numpy.uint32(nprnd(2**30/jobs)),
+                                          numpy.uint32(nprnd(2**30/jobs)))
         
         
       print "%s with %i %s done" % (Alu,jobs,ParaStyle)
@@ -719,9 +600,9 @@
     print Splutter,sum(Splutter)
     #MySplutter[i]=numpy.median(Splutter)
     #print numpy.mean(Splutter)*len(Splutter),MySplutter[i]*len(Splutter),numpy.std(Splutter)
+    
+  SplutterCL.release()
 
-    SplutterCL.release()
-
   print jobs,numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration)
 	
   return(numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration))
@@ -830,18 +711,18 @@
   Curves=False
   # Fit is True to print the curves
   Fit=False
-  # Spluttering is Dense by default
-  Density='All'
+  # Memory of vector explored
+  Memory=1024
 
   try:
-    opts, args = getopt.getopt(sys.argv[1:],"hocfa:g:p:i:s:e:t:r:d:y:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","jobstep=","redo=","device=","density="])
+    opts, args = getopt.getopt(sys.argv[1:],"hocfa:g:p:i:s:e:t:r:d:m:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","jobstep=","redo=","device="])
   except getopt.GetoptError:
-    print '%s -o (Out of Core Metrology) -c (Print Curves) -f (Fit to Amdahl Law) -y <Dense/Sparse/All> -a <CPU/GPU/ACCELERATOR> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Hybrid/Blocks> -i <Iterations> -s <JobStart> -e <JobEnd> -t <JobStep> -r <RedoToImproveStats> ' % sys.argv[0]
+    print '%s -o (Out of Core Metrology) -c (Print Curves) -f (Fit to Amdahl Law) -a <CPU/GPU/ACCELERATOR> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Hybrid/Blocks> -i <Iterations> -s <JobStart> -e <JobEnd> -t <JobStep> -r <RedoToImproveStats> -m <MemoryRaw>' % sys.argv[0]
     sys.exit(2)
     
   for opt, arg in opts:
     if opt == '-h':
-      print '%s -o (Out of Core Metrology) -c (Print Curves) -f (Fit to Amdahl Law)  -y <Dense/Sparse/All> -a <CPU/GPU/ACCELERATOR> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Hybrid/Blocks> -i <Iterations> -s <JobStart> -e <JobEnd> -t <JobStep> -r <RedoToImproveStats>' % sys.argv[0]
+      print '%s -o (Out of Core Metrology) -c (Print Curves) -f (Fit to Amdahl Law) -a <CPU/GPU/ACCELERATOR> -d <DeviceId> -g <CUDA/OpenCL> -p <Threads/Hybrid/Blocks> -i <Iterations> -s <JobStart> -e <JobEnd> -t <JobStep> -r <RedoToImproveStats> -m <MemoryRaw>' % sys.argv[0]
 
       print "\nInformations about devices detected under OpenCL:"
       # For PyOpenCL import
@@ -865,8 +746,6 @@
       Metrology='OutMetro'
     elif opt == '-c':
       Curves=True
-    elif opt in ("-y", "--density"):
-      Density = arg
     elif opt == '-f':
       Fit=True
     elif opt in ("-a", "--alu"):
@@ -887,6 +766,8 @@
       JobStep = int(arg)
     elif opt in ("-r", "--redo"):
       Redo = int(arg)
+    elif opt in ("-m", "--memory"):
+      Memory = int(arg)
 
   print "Toto %s" % Alu
 
@@ -902,11 +783,11 @@
   print "Device Identification : %s" % Device
   print "GpuStyle used : %s" % GpuStyle
   print "Parallel Style used : %s" % ParaStyle
-  print "Density Spluttering : %s" % Density
   print "Iterations : %s" % Iterations
   print "Number of threads on start : %s" % JobStart
   print "Number of threads on end : %s" % JobEnd
   print "Number of redo : %s" % Redo
+  print "Memory  : %s" % Memory
   print "Metrology done out of CPU/GPU : %r" % OutMetrology
 
   if GpuStyle=='CUDA':
@@ -954,13 +835,14 @@
         start=time.time()
         if GpuStyle=='CUDA':
           try:
-            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,Density)
+            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,
+                                 Memory)
           except:
             print "Problem with %i // computations on Cuda" % Jobs
         elif GpuStyle=='OpenCL':
           try:
             a,m,s=MetropolisOpenCL(circle,Iterations,1,Jobs,ParaStyle,
-                                   Alu,Device,Density)
+                                   Alu,Device,Memory)
           except:
             print "Problem with %i // computations on OpenCL" % Jobs            
         duration=numpy.append(duration,time.time()-start)
@@ -974,13 +856,13 @@
       if GpuStyle=='CUDA':
         try:
           avg,med,std=MetropolisCuda(circle,Iterations,Redo,
-                                     Jobs,ParaStyle,Density)
+                                     Jobs,ParaStyle,Memory)
         except:
           print "Problem with %i // computations on Cuda" % Jobs
       elif GpuStyle=='OpenCL':
         try:
-          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,
-                                       Jobs,ParaStyle,Alu,Device,Density)
+          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,
+                                       ParaStyle,Alu,Device,Memory)
         except:
           print "Problem with %i // computations on OpenCL" % Jobs            
 
