Index: Splutter/GPU/SplutterGPU.py
===================================================================
--- Splutter/GPU/SplutterGPU.py	(revision 62)
+++ Splutter/GPU/SplutterGPU.py	(revision 63)
@@ -70,84 +70,162 @@
 
 #define znew  ((z=36969*(z&65535)+(z>>16))<<16)
 #define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
+
 #define MWC   (znew+wnew)
 #define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
 #define CONG  (jcong=69069*jcong+1234567)
 #define KISS  ((MWC^CONG)+SHR3)
 
+#define CONGfp CONG * 2.328306435454494e-10f
+#define SHR3fp SHR3 * 2.328306435454494e-10f
 #define MWCfp MWC * 2.328306435454494e-10f
 #define KISSfp KISS * 2.328306435454494e-10f
 
-#define MAX 4294967296
+#define MAX (ulong)4294967296
 
-uint rotl(uint value, int shift) {
-   return (value << shift) | (value >> (sizeof(value) * 8 - shift));
+__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);
+    const ulong size=(ulong)(gridDim.x*blockDim.x);
+    const ulong block=(ulong)space/(ulong)size;
+   
+    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++) {
+
+      // Dense version
+       uint position=(uint)( ((ulong)MWC+id*MAX)*block/MAX );
+
+      s[position]++;
+   }
+
+   __syncthreads();
 }
- 
-uint rotr(uint value, int shift) {
-   return (value >> shift) | (value << (sizeof(value) * 8 - shift));
+
+__global__ void SplutterGlobalSparse(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);
+    const ulong size=(ulong)(gridDim.x*blockDim.x);
+    const ulong block=(ulong)space/(ulong)size;
+   
+    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++) {
+
+      // Sparse version
+       uint position=(uint)( (ulong)MWC*block/MAX*size+id );
+
+      s[position]++;
+   }
+
+   __syncthreads();
 }
 
-__global__ void MainLoopBlocks(uint *s,uint size,ulong iterations,uint seed_w,uint seed_z)
+__global__ void SplutterLocalDense(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   // uint z=rotl(seed_z,blockIdx.x);
-   // uint w=rotr(seed_w,blockIdx.x);
+    const ulong id=(ulong)(threadIdx.x);
+    const ulong size=(ulong)(blockDim.x);
+    const ulong block=(ulong)space/(ulong)size;
+   
+    uint z=seed_z-(uint)id;
+    uint w=seed_w+(uint)id;
 
-   // uint jsr=rotl(seed_z,blockIdx.x);
-   // uint jcong=rotr(seed_w,blockIdx.x);
+    uint jsr=seed_z;
+    uint jcong=seed_w;
 
-   uint z=seed_z/(blockIdx.x+1);
-   uint w=seed_w%(blockIdx.x+1);
+   for ( ulong i=0;i<iterations;i++) {
 
-   uint jsr=seed_z/(blockIdx.x+1);
-   uint jcong=seed_w%(blockIdx.x+1);
+      // Dense version
+       size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
 
-   for (ulong i=0;i<iterations;i++) {
+      s[position]++;
+   }
 
-      s[(uint)(((ulong)size*(ulong)CONG)/(ulong)MAX)]+=1;
 
-   }
-   __threadfence_block();
+   __syncthreads();
+
 }
 
-__global__ void MainLoopThreads(uint *s,uint size,ulong iterations,uint seed_w,uint seed_z)
-{ 
-   // uint z=rotl(seed_z,threadIdx.x);
-   // uint w=rotr(seed_w,threadIdx.x);
+__global__ void SplutterLocalSparse(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+{
+    const ulong id=(ulong)threadIdx.x;
+    const ulong size=(ulong)blockDim.x;
+    const ulong block=(ulong)space/(ulong)size;
+   
+    uint z=seed_z-(uint)id;
+    uint w=seed_w+(uint)id;
 
-   // uint jsr=rotl(seed_z,threadIdx.x);
-   // uint jcong=rotr(seed_w,threadIdx.x);
+    uint jsr=seed_z;
+    uint jcong=seed_w;
 
-   uint z=seed_z;
-   uint w=seed_w;
+   for ( ulong i=0;i<iterations;i++) {
 
-   uint jsr=seed_z;
-   uint jcong=seed_w;
+      // Sparse version
+       size_t position=(size_t)( (ulong)MWC*block/MAX*size+id );
 
-   for (ulong i=0;i<iterations;i++) {
-
-      s[(uint)(((ulong)size*(ulong)CONG)/(ulong)MAX)]+=1;
+      s[position]++;
    }
 
    __syncthreads();
+
 }
 
-__global__ void MainLoopHybrid(uint *s,uint size,ulong iterations,uint seed_w,uint seed_z)
+__global__ void SplutterHybridDense(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   uint z=seed_z;
-   uint w=seed_w;
+    const ulong id=(ulong)(blockIdx.x);
+    const ulong size=(ulong)(gridDim.x);
+    const ulong block=(ulong)space/(ulong)size;
+   
+    uint z=seed_z-(uint)id;
+    uint w=seed_w+(uint)id;
 
-   uint jsr=seed_z;
-   uint jcong=seed_w;
+    uint jsr=seed_z;
+    uint jcong=seed_w;
 
-   for (ulong i=0;i<iterations;i++) {
+   for ( ulong i=0;i<iterations;i++) {
 
-      s[(uint)(((ulong)size*(ulong)CONG)/(ulong)MAX)]+=1;
+      // Dense version
+      size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
+
+      s[position]++;
    }
+      
+   __syncthreads();
+}
 
+__global__ void SplutterHybridSparse(uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+{
+    const ulong id=(ulong)(blockIdx.x);
+    const ulong size=(ulong)(gridDim.x);
+    const ulong block=(ulong)space/(ulong)size;
+   
+    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++) {
+
+      // Sparse version
+      size_t position=(size_t)( (((ulong)MWC*block)/MAX)*size+id );
+
+      s[position]++;
+
+   }
+
+   //s[blockIdx.x]=blockIdx.x;
    __syncthreads();
+}
 
-}
 """
 
 KERNEL_CODE_OPENCL="""
@@ -175,12 +253,8 @@
     return (value >> shift) | (value << (sizeof(value) * CHAR_BIT - shift));
 }
 
-__kernel void MainLoopGlobal(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
+__kernel void SplutterGlobalDense(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   //__private const float id=(float)get_global_id(0);
-   //__private const float size=(float)get_global_size(0);
-   //__private const float block=space/size;
-
    __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;
@@ -193,25 +267,57 @@
 
    for (__private ulong i=0;i<iterations;i++) {
 
-      // Standard version does not work for several processes (some lost!) memory unconsistent
-      //__private size_t position=(size_t)(((ulong)space*(ulong)MWC)/(ulong)MAX);
-      
       // Dense version
-      //__private size_t position=(size_t)( ((ulong)MWC+(ulong)id*(ulong)MAX)*(ulong)block/(ulong)MAX );
+      __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+(ulong)id*(ulong)MAX)*(ulong)block/(ulong)MAX );
-      //__private size_t position=(size_t)( ((ulong)MWC*(block)+(ulong)id*(ulong)MAX)/(ulong)MAX );
-      // First
-      //__private size_t position=(size_t)( (ulong)(0)*(ulong)size+(ulong)id );
-      // Last
-      //__private size_t position=(size_t)( (ulong)(block-1)*(ulong)size+(ulong)id );
-      // General
-      __private size_t position=(size_t)( (ulong)MWC*(ulong)(block)/(ulong)MAX*(ulong)size+(ulong)id );
+      __private size_t position=(size_t)( (ulong)MWC*block/MAX*size+id );
 
-      // Float version seems to be the best...
-      //__private uint position=(uint)( block*(CONGfp+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;
+
+   __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]++;
    }
 
@@ -219,54 +325,99 @@
 
 }
 
-__kernel void MainLoopLocal(__global uint *s,uint size,ulong iterations,uint seed_w,uint seed_z)
+__kernel void SplutterLocalSparse(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   uint z=rotl(seed_z,get_local_id(0));
-   uint w=rotr(seed_w,get_local_id(0));
+   __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;
 
-   uint jsr=rotl(seed_z,get_local_id(0));
-   uint jcong=rotr(seed_w,get_local_id(0));
+   __private uint jsr=seed_z;
+   __private uint jcong=seed_w;
 
-   for (ulong i=0;i<iterations;i++) {
+   for (__private ulong i=0;i<iterations;i++) {
 
-      s[(int)(((ulong)size*(ulong)CONG)/(ulong)MAX)]+=(uint)1;
+      // 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);
 
-   barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
 }
 
-__kernel void MainLoopHybrid(__global uint *s,uint size,ulong iterations,uint seed_w,uint seed_z)
+__kernel void SplutterHybridDense(__global uint *s,const uint space,const ulong iterations,const uint seed_w,const uint seed_z)
 {
-   uint z=rotl(seed_z,get_group_id(0)*get_num_groups(0)+get_local_id(0));
-   uint w=rotr(seed_w,get_group_id(0)*get_num_groups(0)+get_local_id(0));
+   __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;
 
-   uint jsr=rotl(seed_z,get_group_id(0)*get_num_groups(0)+get_local_id(0));
-   uint jcong=rotr(seed_w,get_group_id(0)*get_num_groups(0)+get_local_id(0));
+   __private uint jsr=seed_z;
+   __private uint jcong=seed_w;
 
-   for (ulong i=0;i<iterations;i++) {
+   for (__private ulong i=0;i<iterations;i++) {
 
-      s[(int)(((ulong)size*(ulong)CONG)/(ulong)MAX)]+=1;
+      // Dense version
+      __private size_t position=(size_t)( ((ulong)MWC+id*MAX)*block/MAX );
 
-      barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
+      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):
+def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,Dense):
 
   # Avec PyCUDA autoinit, rien a faire !
   
   circleCU = cuda.InOut(circle)
   
+  print "prout"
+
   mod = SourceModule(KERNEL_CODE_CUDA)
 
-  MetropolisBlocksCU=mod.get_function("MainLoopBlocks")
-  MetropolisJobsCU=mod.get_function("MainLoopThreads")
-  MetropolisHybridCU=mod.get_function("MainLoopHybrid")
+  print "prout 2"
+
+  if Dense:
+    MetropolisBlocksCU=mod.get_function("SplutterGlobalDense")
+    MetropolisThreadsCU=mod.get_function("SplutterLocalDense")
+    MetropolisHybridCU=mod.get_function("SplutterHybridDense")
+  else:
+    MetropolisBlocksCU=mod.get_function("SplutterGlobalSparse")
+    MetropolisThreadsCU=mod.get_function("SplutterLocalSparse")
+    MetropolisHybridCU=mod.get_function("SplutterHybridSparse")
   
+  print "prout 3"
+
   start = pycuda.driver.Event()
   stop = pycuda.driver.Event()
   
@@ -280,11 +431,13 @@
     
   iterationsNew=iterationsCL*jobs
 
+  Splutter=numpy.zeros(jobs*16).astype(numpy.uint32)
+
   for i in range(steps):
 
-    Splutter=numpy.zeros(1024).astype(numpy.uint32)
+    Splutter[:]=0
     
-    print Splutter
+    print Splutter,len(Splutter)
 
     SplutterCU = cuda.InOut(Splutter)
 
@@ -313,7 +466,7 @@
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs/threads,threads,ParaStyle)
     else:
-      MetropolisJobsCU(SplutterCU,
+      MetropolisThreadsCU(SplutterCU,
                        numpy.uint32(len(Splutter)),
                        numpy.uint64(iterationsCL),
                        numpy.uint32(nprnd(2**30/jobs)),
@@ -321,7 +474,7 @@
                        grid=(1,1),
                        block=(jobs,1,1))
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
-            (Alu,jobs,1,ParaStyle)
+            (Alu,1,jobs,ParaStyle)
     stop.record()
     stop.synchronize()
                 
@@ -343,7 +496,7 @@
   return(numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration))
 
 
-def MetropolisOpenCL(circle,iterations,steps,jobs,ParaStyle,Alu,Device):
+def MetropolisOpenCL(circle,iterations,steps,jobs,ParaStyle,Alu,Device,Dense):
 	
   # Initialisation des variables en les CASTant correctement
 
@@ -401,21 +554,20 @@
 
   # Je cree le contexte et la queue pour son execution
   ctx = cl.Context([XPU])
-  queue = cl.CommandQueue(ctx,
-                          properties=cl.command_queue_properties.PROFILING_ENABLE)
-
+  queue = cl.CommandQueue(ctx,properties=cl.command_queue_properties.PROFILING_ENABLE)
+  
   # Je recupere les flag possibles pour les buffers
   mf = cl.mem_flags
-	
+
   MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build(options = "-cl-mad-enable -cl-fast-relaxed-math")
-  
+      
   MyDuration=numpy.zeros(steps)
   
   if iterations%jobs==0:
     iterationsCL=numpy.uint64(iterations/jobs)
   else:
     iterationsCL=numpy.uint64(iterations/jobs+1)
-
+    
   iterationsNew=numpy.uint64(iterationsCL*jobs)
 
   MySplutter=numpy.zeros(steps)
@@ -453,34 +605,60 @@
       #                                      numpy.uint64(iterationsCL),
       #                                      numpy.uint32(nprnd(2**30/jobs)),
       #                                      numpy.uint32(nprnd(2**30/jobs)))
-      CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
-                                           SplutterCL,
-                                           numpy.uint32(len(Splutter)),
-                                           numpy.uint64(iterationsCL),
-                                           numpy.uint32(521288629),
-                                           numpy.uint32(362436069))
+      if Dense:
+        CLLaunch=MetropolisCL.SplutterGlobalDense(queue,(jobs,),None,
+                                                  SplutterCL,
+                                                  numpy.uint32(len(Splutter)),
+                                                  numpy.uint64(iterationsCL),
+                                                  numpy.uint32(521288629),
+                                                  numpy.uint32(362436069))
+      else:
+        CLLaunch=MetropolisCL.SplutterGlobalSparse(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
-      CLLaunch=MetropolisCL.MainLoopHybrid(queue,(jobs,),(threads,),
-                                           SplutterCL,
-                                           numpy.uint32(len(Splutter)),
-                                           numpy.uint64(iterationsCL),
-                                           numpy.uint32(nprnd(2**30/jobs)),
-                                           numpy.uint32(nprnd(2**30/jobs)))
+      if 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:
+        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)))
         
       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
-      CLLaunch=MetropolisCL.MainLoopLocal(queue,(jobs,),(jobs,),
-                                          SplutterCL,
-                                          numpy.uint32(len(Splutter)),
-                                          numpy.uint64(iterationsCL),
-                                          numpy.uint32(nprnd(2**30/jobs)),
-                                          numpy.uint32(nprnd(2**30/jobs)))
+      if 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:
+        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)))
+        
       print "%s with %i %s done" % (Alu,jobs,ParaStyle)
 
     CLLaunch.wait()
@@ -603,16 +781,18 @@
   Curves=False
   # Fit is True to print the curves
   Fit=False
+  # Spluttering is Dense by default
+  Dense=True
 
   try:
-    opts, args = getopt.getopt(sys.argv[1:],"hoclfa:g:p:i:s:e:t:r:d:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","jobstep=","redo=","device="])
+    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="])
   except getopt.GetoptError:
-    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> ' % sys.argv[0]
+    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]
     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) -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)  -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 "\nInformations about devices detected under OpenCL:"
       # For PyOpenCL import
@@ -636,6 +816,10 @@
       Metrology='OutMetro'
     elif opt == '-c':
       Curves=True
+    elif opt == '-v':
+      Dense=True
+    elif opt == '-w':
+      Dense=False
     elif opt == '-f':
       Fit=True
     elif opt in ("-a", "--alu"):
@@ -669,6 +853,7 @@
   print "Device Identification : %s" % Device
   print "GpuStyle used : %s" % GpuStyle
   print "Parallel Style used : %s" % ParaStyle
+  print "Dense (or Sparse) Spluttering : %r" % Dense
   print "Iterations : %s" % Iterations
   print "Number of threads on start : %s" % JobStart
   print "Number of threads on end : %s" % JobEnd
@@ -720,13 +905,14 @@
         start=time.time()
         if GpuStyle=='CUDA':
           try:
-            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle)
+            print "toto"
+            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,Dense)
           except:
             print "Problem with %i // computations on Cuda" % Jobs
         elif GpuStyle=='OpenCL':
           try:
             a,m,s=MetropolisOpenCL(circle,Iterations,1,Jobs,ParaStyle,
-                                   Alu,Device)
+                                   Alu,Device,Dense)
           except:
             print "Problem with %i // computations on OpenCL" % Jobs            
         duration=numpy.append(duration,time.time()-start)
@@ -739,12 +925,12 @@
     else:
       if GpuStyle=='CUDA':
         try:
-          avg,med,std=MetropolisCuda(circle,Iterations,Redo,Jobs,ParaStyle)
+          avg,med,std=MetropolisCuda(circle,Iterations,Redo,Jobs,ParaStyle,Dense)
         except:
           print "Problem with %i // computations on Cuda" % Jobs
       elif GpuStyle=='OpenCL':
         try:
-          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device)
+          avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device,Dense)
         except:
           print "Problem with %i // computations on OpenCL" % Jobs            
 
