Index: Splutter/GPU/SplutterGPU.py
===================================================================
--- Splutter/GPU/SplutterGPU.py	(revision 65)
+++ Splutter/GPU/SplutterGPU.py	(revision 66)
@@ -7,12 +7,15 @@
 #
 # Thanks to Andreas Klockner for PyCUDA:
 # http://mathema.tician.de/software/pycuda
+# http://mathema.tician.de/software/pyopencl
 # 
 
 # 2013-01-01 : problems with launch timeout
 # http://stackoverflow.com/questions/497685/how-do-you-get-around-the-maximum-cuda-run-time
 # Option "Interactive" "0" in /etc/X11/xorg.conf
 
+# Marsaglia elements about RNG 
+
 # Common tools
 import numpy
 from numpy.random import randint as nprnd
@@ -64,6 +67,11 @@
 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
@@ -82,7 +90,30 @@
 #define KISSfp KISS * 2.328306435454494e-10f
 
 #define MAX (ulong)4294967296
+#define UMAX (uint)2147483648
 
+__global__ void SplutterGlobal(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+{
+    const ulong id=(ulong)(blockIdx.x);
+   
+    uint z=seed_z-(uint)id;
+    uint w=seed_w+(uint)id;
+
+    uint jsr=seed_z;
+    uint jcong=seed_w;
+
+   for ( ulong i=0;i<iterations;i++) {
+
+      // All version
+      uint position=(uint)( ((ulong)MWC*(ulong)space)/MAX );
+
+      // UMAX is set to avoid round over overflow
+      atomicInc(&s[position],UMAX);
+   }
+
+   __syncthreads();
+}
+
 __global__ void SplutterGlobalDense(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
     const ulong id=(ulong)(threadIdx.x+blockIdx.x*blockDim.x);
@@ -229,6 +260,8 @@
 """
 
 KERNEL_CODE_OPENCL="""
+#pragma OPENCL EXTENSION cl_khr_global_int32_base_atomics : enable
+
 // Marsaglia RNG very simple implementation
 #define znew  ((z=36969*(z&65535)+(z>>16))<<16)
 #define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
@@ -253,6 +286,30 @@
     return (value >> shift) | (value << (sizeof(value) * CHAR_BIT - shift));
 }
 
+__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;
+
+   __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*(ulong)space)/MAX );
+
+      atomic_inc(&s[position]);
+   }
+
+   barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
+
+}
+
 __kernel void SplutterGlobalDense(__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);
@@ -395,29 +452,25 @@
 
 """
 
-def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,Dense):
+def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,Density):
 
   # Avec PyCUDA autoinit, rien a faire !
-  
+
   circleCU = cuda.InOut(circle)
   
-  print "prout"
-
   mod = SourceModule(KERNEL_CODE_CUDA)
 
-  print "prout 2"
-
-  if Dense:
+  if Density=='Dense':
     MetropolisBlocksCU=mod.get_function("SplutterGlobalDense")
     MetropolisThreadsCU=mod.get_function("SplutterLocalDense")
     MetropolisHybridCU=mod.get_function("SplutterHybridDense")
-  else:
+  elif Density=='Sparse':
     MetropolisBlocksCU=mod.get_function("SplutterGlobalSparse")
     MetropolisThreadsCU=mod.get_function("SplutterLocalSparse")
     MetropolisHybridCU=mod.get_function("SplutterHybridSparse")
-  
-  print "prout 3"
-
+  else:
+    MetropolisBlocksCU=mod.get_function("SplutterGlobal")
+    
   start = pycuda.driver.Event()
   stop = pycuda.driver.Event()
   
@@ -496,7 +549,8 @@
   return(numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration))
 
 
-def MetropolisOpenCL(circle,iterations,steps,jobs,ParaStyle,Alu,Device,Dense):
+def MetropolisOpenCL(circle,iterations,steps,jobs,
+                     ParaStyle,Alu,Device,Density):
 	
   # Initialisation des variables en les CASTant correctement
 
@@ -605,14 +659,14 @@
       #                                      numpy.uint64(iterationsCL),
       #                                      numpy.uint32(nprnd(2**30/jobs)),
       #                                      numpy.uint32(nprnd(2**30/jobs)))
-      if Dense:
+      if Density=='Dense':
         CLLaunch=MetropolisCL.SplutterGlobalDense(queue,(jobs,),None,
                                                   SplutterCL,
                                                   numpy.uint32(len(Splutter)),
                                                   numpy.uint64(iterationsCL),
                                                   numpy.uint32(521288629),
                                                   numpy.uint32(362436069))
-      else:
+      elif Density=='Sparse':
         CLLaunch=MetropolisCL.SplutterGlobalSparse(queue,(jobs,),None,
                                                    SplutterCL,
                                                    numpy.uint32(len(Splutter)),
@@ -620,19 +674,27 @@
                                                    numpy.uint32(521288629),
                                                    numpy.uint32(362436069))
         
+      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)
       # en OpenCL, necessaire de mettre un Global_id identique au local_id
-      if Dense:
+      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)))
-      else:
+      elif Density=='Sparse':
         CLLaunch=MetropolisCL.SplutterHybridSparse(queue,(jobs,),(threads,),
                                                    SplutterCL,
                                                    numpy.uint32(len(Splutter)),
@@ -643,15 +705,15 @@
       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 Dense:
+      # 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)))
-      else:
+      elif Density=='Sparse':
         CLLaunch=MetropolisCL.SplutterLocalSparse(queue,(jobs,),(jobs,),
                                                   SplutterCL,
                                                   numpy.uint32(len(Splutter)),
@@ -659,6 +721,7 @@
                                                   numpy.uint32(nprnd(2**30/jobs)),
                                                   numpy.uint32(nprnd(2**30/jobs)))
         
+        
       print "%s with %i %s done" % (Alu,jobs,ParaStyle)
 
     CLLaunch.wait()
@@ -765,7 +828,7 @@
   # Parallel distribution can be on Threads or Blocks
   ParaStyle='Blocks'
   # Iterations is integer
-  Iterations=100000000
+  Iterations=10000000
   # JobStart in first number of Jobs to explore
   JobStart=1
   # JobEnd is last number of Jobs to explore
@@ -782,17 +845,17 @@
   # Fit is True to print the curves
   Fit=False
   # Spluttering is Dense by default
-  Dense=True
+  Density='All'
 
   try:
-    opts, args = getopt.getopt(sys.argv[1:],"hocfvwa:g:p:i:s:e:t:r:d:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","jobstep=","redo=","device="])
+    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="])
   except getopt.GetoptError:
-    print '%s -o (Out of Core Metrology) -c (Print Curves) -f (Fit to Amdahl Law) -v (Dense Spluttering) -w (Sparse Spluttering) -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) -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]
     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)  -v (Dense Spluttering) -w (Sparse Spluttering) -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)  -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 "\nInformations about devices detected under OpenCL:"
       # For PyOpenCL import
@@ -816,10 +879,8 @@
       Metrology='OutMetro'
     elif opt == '-c':
       Curves=True
-    elif opt == '-v':
-      Dense=True
-    elif opt == '-w':
-      Dense=False
+    elif opt in ("-y", "--density"):
+      Density = arg
     elif opt == '-f':
       Fit=True
     elif opt in ("-a", "--alu"):
@@ -841,19 +902,21 @@
     elif opt in ("-r", "--redo"):
       Redo = int(arg)
 
+  print "Toto %s" % Alu
+
   if Alu=='CPU' and GpuStyle=='CUDA':
     print "Alu can't be CPU for CUDA, set Alu to GPU"
     Alu='GPU'
 
   if ParaStyle not in ('Blocks','Threads','Hybrid'):
     print "%s not exists, ParaStyle set as Threads !" % ParaStyle
-    ParaStyle='Threads'
+    ParaStyle='Blocks'
 
   print "Compute unit : %s" % Alu
   print "Device Identification : %s" % Device
   print "GpuStyle used : %s" % GpuStyle
   print "Parallel Style used : %s" % ParaStyle
-  print "Dense (or Sparse) Spluttering : %r" % Dense
+  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
@@ -905,14 +968,13 @@
         start=time.time()
         if GpuStyle=='CUDA':
           try:
-            print "toto"
-            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,Dense)
+            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,Density)
           except:
             print "Problem with %i // computations on Cuda" % Jobs
         elif GpuStyle=='OpenCL':
           try:
             a,m,s=MetropolisOpenCL(circle,Iterations,1,Jobs,ParaStyle,
-                                   Alu,Device,Dense)
+                                   Alu,Device,Density)
           except:
             print "Problem with %i // computations on OpenCL" % Jobs            
         duration=numpy.append(duration,time.time()-start)
@@ -925,12 +987,14 @@
     else:
       if GpuStyle=='CUDA':
         try:
-          avg,med,std=MetropolisCuda(circle,Iterations,Redo,Jobs,ParaStyle,Dense)
+          avg,med,std=MetropolisCuda(circle,Iterations,Redo,
+                                     Jobs,ParaStyle,Density)
         except:
           print "Problem with %i // computations on Cuda" % Jobs
       elif GpuStyle=='OpenCL':
         try:
-          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device,Dense)
+          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,
+                                       Jobs,ParaStyle,Alu,Device,Density)
         except:
           print "Problem with %i // computations on OpenCL" % Jobs            
 
