Index: Pi/GPU/Pi-GPU.py
===================================================================
--- Pi/GPU/Pi-GPU.py	(revision 49)
+++ Pi/GPU/Pi-GPU.py	(revision 50)
@@ -148,9 +148,79 @@
    __syncthreads();
 
 }
+
+__global__ void MainLoopBlocks64(ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(blockIdx.x+1);
+   uint w=seed_w/(blockIdx.x+1);
+
+   ulong total=0;
+
+   for (ulong i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+
+   }
+
+   s[blockIdx.x]=total;
+   __syncthreads();
+
+}
+
+__global__ void MainLoopThreads64(ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(threadIdx.x+1);
+   uint w=seed_w/(threadIdx.x+1);
+
+   ulong total=0;
+
+   for (ulong i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+
+   }
+
+   s[threadIdx.x]=total;
+   __syncthreads();
+
+}
+
+__global__ void MainLoopHybrid64(ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(blockDim.x*blockIdx.x+threadIdx.x+1);
+   uint w=seed_w/(blockDim.x*blockIdx.x+threadIdx.x+1);
+
+   ulong total=0;
+
+   for (ulong i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+
+   }
+
+   s[blockDim.x*blockIdx.x+threadIdx.x]=total;
+   __syncthreads();
+
+}
 """
 
 KERNEL_CODE_OPENCL="""
+#pragma OPENCL EXTENSION cl_khr_fp64: enable
 
 // Marsaglia RNG very simple implementation
 #define znew  ((z=36969*(z&65535)+(z>>16))<<16)
@@ -225,9 +295,72 @@
    s[get_group_id(0)*get_num_groups(0)+get_local_id(0)]=total;
       
 }
+
+__kernel void MainLoopGlobal64(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_global_id(0)+1);
+   uint w=seed_w/(get_global_id(0)+1);
+
+   ulong total=0;
+
+   for (ulong i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+   }
+   s[get_global_id(0)]=total;
+   barrier(CLK_GLOBAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopLocal64(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_local_id(0)+1);
+   uint w=seed_w/(get_local_id(0)+1);
+
+   ulong total=0;
+
+   for (ulong i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+   }
+   s[get_local_id(0)]=total;
+   barrier(CLK_LOCAL_MEM_FENCE);
+      
+}
+
+__kernel void MainLoopHybrid64(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)
+{
+   uint z=seed_z/(get_group_id(0)*get_num_groups(0)+get_local_id(0)+1);
+   uint w=seed_w/(get_group_id(0)*get_num_groups(0)+get_local_id(0)+1);
+
+   ulong total=0;
+
+   for (uint i=0;i<iterations;i++) {
+
+      double x=(double)MWCfp ;
+      double y=(double)MWCfp ;
+
+      // Matching test
+      ulong inside=((x*x+y*y) < 1.0f) ? 1:0;
+      total+=inside;
+   }
+   barrier(CLK_LOCAL_MEM_FENCE);
+   s[get_group_id(0)*get_num_groups(0)+get_local_id(0)]=total;
+      
+}
 """
 
-def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle):
+def MetropolisCuda(circle,iterations,steps,jobs,ParaStyle,DoublePrecision):
 
   # Avec PyCUDA autoinit, rien a faire !
   
@@ -238,13 +371,16 @@
   MetropolisBlocksCU=mod.get_function("MainLoopBlocks")
   MetropolisJobsCU=mod.get_function("MainLoopThreads")
   MetropolisHybridCU=mod.get_function("MainLoopHybrid")
+  MetropolisBlocks64CU=mod.get_function("MainLoopBlocks64")
+  MetropolisJobs64CU=mod.get_function("MainLoopThreads64")
+  MetropolisHybrid64CU=mod.get_function("MainLoopHybrid64")
   
   start = pycuda.driver.Event()
   stop = pycuda.driver.Event()
   
   MyPi=numpy.zeros(steps)
   MyDuration=numpy.zeros(steps)
-  
+
   if iterations%jobs==0:
     iterationsCL=numpy.uint64(iterations/jobs)
     iterationsNew=iterationsCL*jobs
@@ -256,51 +392,76 @@
     start.record()
     start.synchronize()
     if ParaStyle=='Blocks':
-      MetropolisBlocksCU(circleCU,
-                         numpy.uint64(iterationsCL),
-                         numpy.uint32(nprnd(2**30/jobs)),
-                         numpy.uint32(nprnd(2**30/jobs)),
-                         grid=(jobs,1),
-                         block=(1,1,1))
+      if DoublePrecision:
+        MetropolisBlocksCU(circleCU,
+                           numpy.uint64(iterationsCL),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           grid=(jobs,1),
+                           block=(1,1,1))
+      else:
+        MetropolisBlocks64CU(circleCU,
+                             numpy.uint64(iterationsCL),
+                             numpy.uint32(nprnd(2**30/jobs)),
+                             numpy.uint32(nprnd(2**30/jobs)),
+                             grid=(jobs,1),
+                             block=(1,1,1))
+        
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs,1,ParaStyle)      
     elif ParaStyle=='Hybrid':
       threads=BestThreadsNumber(jobs)
-      MetropolisHybridCU(circleCU,
-                          numpy.uint64(iterationsCL),
-                          numpy.uint32(nprnd(2**30/jobs)),
-                          numpy.uint32(nprnd(2**30/jobs)),
-                          grid=(jobs,1),
-                          block=(threads,1,1))
+      if DoublePrecision:
+        MetropolisHybrid64CU(circleCU,
+                             numpy.uint64(iterationsCL),
+                             numpy.uint32(nprnd(2**30/jobs)),
+                             numpy.uint32(nprnd(2**30/jobs)),
+                             grid=(jobs,1),
+                             block=(threads,1,1))
+      else:
+        MetropolisHybridCU(circleCU,
+                           numpy.uint64(iterationsCL),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           grid=(jobs,1),
+                           block=(threads,1,1))
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs/threads,threads,ParaStyle)
     else:
-      MetropolisJobsCU(circleCU,
-                          numpy.uint64(iterationsCL),
-                          numpy.uint32(nprnd(2**30/jobs)),
-                          numpy.uint32(nprnd(2**30/jobs)),
-                          grid=(1,1),
-                          block=(jobs,1,1))
+      if DoublePrecision:
+        MetropolisJobs64CU(circleCU,
+                           numpy.uint64(iterationsCL),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           numpy.uint32(nprnd(2**30/jobs)),
+                           grid=(1,1),
+                           block=(jobs,1,1))
+      else:
+        MetropolisJobsCU(circleCU,
+                         numpy.uint64(iterationsCL),
+                         numpy.uint32(nprnd(2**30/jobs)),
+                         numpy.uint32(nprnd(2**30/jobs)),
+                         grid=(1,1),
+                         block=(jobs,1,1))
       print "%s with (WorkItems/Threads)=(%i,%i) %s method done" % \
             (Alu,jobs,1,ParaStyle)
     stop.record()
     stop.synchronize()
                 
-    #elapsed = stop.time_since(start)*1e-3
     elapsed = start.time_till(stop)*1e-3
 
-    #print circle,float(numpy.sum(circle))
-    MyPi[i]=4.*float(numpy.sum(circle))/float(iterationsCL)
     MyDuration[i]=elapsed
-    #print MyPi[i],MyDuration[i]
-    #time.sleep(1)
+    AllPi=4./numpy.float32(iterationsCL)*circle.astype(numpy.float32)
+    MyPi[i]=numpy.median(AllPi)
+    print MyPi[i],numpy.std(AllPi),MyDuration[i]
 
+
   print jobs,numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration)
 
   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,
+                     DoublePrecision):
 	
   # Initialisation des variables en les CASTant correctement
     
@@ -338,7 +499,6 @@
       return(0,0,0)
 				
   # Je cree le contexte et la queue pour son execution
-  #ctx = cl.create_some_context()
   ctx = cl.Context([XPU])
   queue = cl.CommandQueue(ctx,
                           properties=cl.command_queue_properties.PROFILING_ENABLE)
@@ -351,8 +511,6 @@
   MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build( \
     options = "-cl-mad-enable -cl-fast-relaxed-math")
 
-  #MetropolisCL = cl.Program(ctx,KERNEL_CODE_OPENCL).build()
-
   i=0
 
   MyPi=numpy.zeros(steps)
@@ -375,30 +533,52 @@
       # SeedZCL is lattice translated in CL format
       # SeedWCL is lattice translated in CL format
       # step is number of iterations
-      CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
-                                           circleCL,
-                                           numpy.uint64(iterationsCL),
-                                           numpy.uint32(nprnd(2**30/jobs)),
-                                           numpy.uint32(nprnd(2**30/jobs)))
+      if DoublePrecision:
+        CLLaunch=MetropolisCL.MainLoopGlobal64(queue,(jobs,),None,
+                                               circleCL,
+                                               numpy.uint64(iterationsCL),
+                                               numpy.uint32(nprnd(2**30/jobs)),
+                                               numpy.uint32(nprnd(2**30/jobs)))
+      else:
+        CLLaunch=MetropolisCL.MainLoopGlobal(queue,(jobs,),None,
+                                             circleCL,
+                                             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,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,),
-                                          circleCL,
-                                          numpy.uint64(iterationsCL),
-                                          numpy.uint32(nprnd(2**30/jobs)),
-                                          numpy.uint32(nprnd(2**30/jobs)))
+      if DoublePrecision:
+        CLLaunch=MetropolisCL.MainLoopHybrid64(queue,(jobs,),(threads,),
+                                               circleCL,
+                                               numpy.uint64(iterationsCL),
+                                               numpy.uint32(nprnd(2**30/jobs)),
+                                               numpy.uint32(nprnd(2**30/jobs)))
+      else:
+        CLLaunch=MetropolisCL.MainLoopHybrid(queue,(jobs,),(threads,),
+                                             circleCL,
+                                             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,),
-                                          circleCL,
-                                          numpy.uint64(iterationsCL),
-                                          numpy.uint32(nprnd(2**30/jobs)),
-                                          numpy.uint32(nprnd(2**30/jobs)))
+      if DoublePrecision:
+        CLLaunch=MetropolisCL.MainLoopLocal64(queue,(jobs,),(jobs,),
+                                              circleCL,
+                                              numpy.uint64(iterationsCL),
+                                              numpy.uint32(nprnd(2**30/jobs)),
+                                              numpy.uint32(nprnd(2**30/jobs)))
+      else:
+        CLLaunch=MetropolisCL.MainLoopLocal(queue,(jobs,),(jobs,),
+                                            circleCL,
+                                            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()
@@ -515,16 +695,18 @@
   Metrology='InMetro'
   # Curves is True to print the curves
   Curves=False
+  # DoublePrecision on FP calculus
+  DoublePrecision=False
 
   try:
-    opts, args = getopt.getopt(sys.argv[1:],"hoca:g:p:i:s:e:t:r:d:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","jobstep=","redo=","device="])
+    opts, args = getopt.getopt(sys.argv[1:],"hocla: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) -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) -l (Double Precision) -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) -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) -l (Double Precision) -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
@@ -540,6 +722,8 @@
     elif opt == '-o':
       OutMetrology=True
       Metrology='OutMetro'
+    elif opt == '-l':
+      DoublePrecision=True
     elif opt == '-c':
       Curves=True
     elif opt in ("-a", "--alu"):
@@ -578,6 +762,7 @@
   print "Number of threads on end : %s" % JobEnd
   print "Number of redo : %s" % Redo
   print "Metrology done out of CPU/GPU : %r" % OutMetrology
+  print "Double Precision in Kernels : %r" % DoublePrecision
 
   if GpuStyle=='CUDA':
     # For PyCUDA import
@@ -594,6 +779,9 @@
       for device in platform.get_devices():
         deviceType=cl.device_type.to_string(device.type)
         print "Device #%i of type %s : %s" % (Id,deviceType,device.name)
+        if Id == Device:
+          # Set the Alu as detected Device Type
+          Alu=deviceType
         Id=Id+1
 
   average=numpy.array([]).astype(numpy.float32)
@@ -615,12 +803,14 @@
         start=time.time()
         if GpuStyle=='CUDA':
           try:
-            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle)
+            a,m,s=MetropolisCuda(circle,Iterations,1,Jobs,ParaStyle,
+                                 DoublePrecision)
           except:
             print "Problem with %i // computations on Cuda" % Jobs
         elif GpuStyle=='OpenCL':
           try:
-            a,m,s=MetropolisOpenCL(circle,Iterations,1,Jobs,ParaStyle,Alu,Device)
+            a,m,s=MetropolisOpenCL(circle,Iterations,1,Jobs,ParaStyle,
+                                   Alu,Device,DoublePrecision)
           except:
             print "Problem with %i // computations on OpenCL" % Jobs            
         duration=numpy.append(duration,time.time()-start)
@@ -633,7 +823,8 @@
     else:
       if GpuStyle=='CUDA':
         try:
-          avg,med,std=MetropolisCuda(circle,Iterations,Redo,Jobs,ParaStyle)
+          avg,med,std=MetropolisCuda(circle,Iterations,Redo,Jobs,ParaStyle,
+                                     DoublePrecision)
         except:
           print "Problem with %i // computations on Cuda" % Jobs
       elif GpuStyle=='OpenCL':
@@ -641,7 +832,7 @@
         #   avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device)
         # except:
         #   print "Problem with %i // computations on OpenCL" % Jobs            
-        avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device)
+        avg,med,std=MetropolisOpenCL(circle,Iterations,Redo,Jobs,ParaStyle,Alu,Device,DoublePrecision)
 
     if (avg,med,std) != (0,0,0):
       print "jobs,avg,med,std",Jobs,avg,med,std
@@ -651,9 +842,14 @@
     else:
       print "Values seem to be wrong..."
     #THREADS*=2
+    if DoublePrecision:
+      Precision='DP'
+    else:
+      Precision='SP'
     if len(average)!=0:
-      numpy.savez("Pi_%s_%s_%s_%s_%i_%.8i_Device%i_%s_%s" % (Alu,GpuStyle,ParaStyle,JobStart,JobEnd,Iterations,Device,Metrology,gethostname()),(ExploredJobs,average,median,stddev))
-      numpy.savetxt("Pi_%s_%s_%s_%s_%i_%.8i_Device%i_%s_%s" % (Alu,GpuStyle,ParaStyle,JobStart,JobEnd,Iterations,Device,Metrology,gethostname()),(ExploredJobs,average,median,stddev))
+      numpy.savez("Pi%s_%s_%s_%s_%s_%i_%.8i_Device%i_%s_%s" % (Precision,Alu,GpuStyle,ParaStyle,JobStart,JobEnd,Iterations,Device,Metrology,gethostname()),(ExploredJobs,average,median,stddev))
+      ToSave=[ ExploredJobs,average,median,stddev ]
+      numpy.savetxt("Pi%s_%s_%s_%s_%s_%i_%.8i_Device%i_%s_%s" % (Precision,Alu,GpuStyle,ParaStyle,JobStart,JobEnd,Iterations,Device,Metrology,gethostname()),numpy.transpose(ToSave))
     Jobs+=JobStep
 
   FitAndPrint(ExploredJobs,median,Curves)
