Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python3

#

# Pi-by-MonteCarlo using PyCUDA/PyOpenCL

#

# CC BY-NC-SA 2011 : Emmanuel QUEMENER <emmanuel.quemener@gmail.com> 

# Cecill v2 : Emmanuel QUEMENER <emmanuel.quemener@gmail.com>

#

# Thanks to Andreas Klockner for PyCUDA:

# http://mathema.tician.de/software/pycuda

# Thanks to Andreas Klockner for PyOpenCL:

# 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

# Common tools

import numpy

from numpy.random import randint as nprnd

import sys

import getopt

import time

import itertools

from socket import gethostname

class PenStacle:

    """Pentacle of Statistics from data"""

    Avg=0

    Med=0

    Std=0

    Min=0

    Max=0

    def __init__(self,Data):

        self.Avg=numpy.average(Data)

        self.Med=numpy.median(Data)

        self.Std=numpy.std(Data)

        self.Max=numpy.max(Data)

        self.Min=numpy.min(Data)

    def display(self):

        print("%s %s %s %s %s" % (self.Avg,self.Med,self.Std,self.Min,self.Max))

class Experience:

    """Metrology for experiences"""

    DeviceStyle=''

    DeviceId=0

    AvgD=0

    MedD=0

    StdD=0

    MinD=0

    MaxD=0

    AvgR=0

    MedR=0

    StdR=0

    MinR=0

    MaxR=0

    def __init__(self,DeviceStyle,DeviceId,Iterations):

        self.DeviceStyle=DeviceStyle

        self.DeviceId=DeviceId

        self.Iterations

    def Metrology(self,Data):

        Duration=PenStacle(Data)

        Rate=PenStacle(Iterations/Data)

        print("Duration %s" % Duration)

        print("Rate %s" % Rate)

def DictionariesAPI():

    Marsaglia={'CONG':0,'SHR3':1,'MWC':2,'KISS':3}

    Computing={'INT32':0,'INT64':1,'FP32':2,'FP64':3}

    return(Marsaglia,Computing)

# 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):

    factors=PrimeFactors(jobs)

    matrix=numpy.append([factors],[factors[::-1]],axis=0)

    threads=1

    for factor in matrix.transpose().ravel():

        threads=threads*factor

        if threads*threads>jobs or threads>512:

            break

    return(long(threads))

# Predicted Amdahl Law (Reduced with s=1-p)  

def AmdahlR(N, T1, p):

    return (T1*(1-p+p/N))

# Predicted Amdahl Law

def Amdahl(N, T1, s, p):

    return (T1*(s+p/N))

# Predicted Mylq Law with first order

def Mylq(N, T1,s,c,p):

    return (T1*(s+p/N)+c*N)

# Predicted Mylq Law with second order

def Mylq2(N, T1,s,c1,c2,p):

    return (T1*(s+p/N)+c1*N+c2*N*N)

def KernelCodeCuda():

    KERNEL_CODE_CUDA="""

#define TCONG 0

#define TSHR3 1

#define TMWC 2

#define TKISS 3

#define TINT32 0

#define TINT64 1

#define TFP32 2

#define TFP64 3

// Marsaglia RNG very simple implementation

#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 MWCfp MWC * 2.328306435454494e-10f

#define KISSfp KISS * 2.328306435454494e-10f

#define SHR3fp SHR3 * 2.328306435454494e-10f

#define CONGfp CONG * 2.328306435454494e-10f

__device__ ulong MainLoop(ulong iterations,uint seed_w,uint seed_z,size_t work)

{

#if TRNG == TCONG

   uint jcong=seed_z+work;

#elif TRNG == TSHR3

   uint jsr=seed_w+work;

#elif TRNG == TMWC

   uint z=seed_z+work;

   uint w=seed_w+work;

#elif TRNG == TKISS

   uint jcong=seed_z+work;

   uint jsr=seed_w+work;

   uint z=seed_z-work;

   uint w=seed_w-work;

#endif

   ulong total=0;

   for (ulong i=0;i<iterations;i++) {

#if TYPE == TINT32

    #define THEONE 1073741824

    #if TRNG == TCONG

        uint x=CONG>>17 ;

        uint y=CONG>>17 ;

    #elif TRNG == TSHR3

        uint x=SHR3>>17 ;

        uint y=SHR3>>17 ;

    #elif TRNG == TMWC

        uint x=MWC>>17 ;

        uint y=MWC>>17 ;

    #elif TRNG == TKISS

        uint x=KISS>>17 ;

        uint y=KISS>>17 ;

    #endif

#elif TYPE == TINT64

    #define THEONE 4611686018427387904

    #if TRNG == TCONG

        ulong x=(ulong)(CONG>>1) ;

        ulong y=(ulong)(CONG>>1) ;

    #elif TRNG == TSHR3

        ulong x=(ulong)(SHR3>>1) ;

        ulong y=(ulong)(SHR3>>1) ;

    #elif TRNG == TMWC

        ulong x=(ulong)(MWC>>1) ;

        ulong y=(ulong)(MWC>>1) ;

    #elif TRNG == TKISS

        ulong x=(ulong)(KISS>>1) ;

        ulong y=(ulong)(KISS>>1) ;

    #endif

#elif TYPE == TFP32

    #define THEONE 1.0f

    #if TRNG == TCONG

        float x=CONGfp ;

        float y=CONGfp ;

    #elif TRNG == TSHR3

        float x=SHR3fp ;

        float y=SHR3fp ;

    #elif TRNG == TMWC

        float x=MWCfp ;

        float y=MWCfp ;

    #elif TRNG == TKISS

      float x=KISSfp ;

      float y=KISSfp ;

    #endif

#elif TYPE == TFP64

    #define THEONE 1.0f

    #if TRNG == TCONG

        double x=(double)CONGfp ;

        double y=(double)CONGfp ;

    #elif TRNG == TSHR3

        double x=(double)SHR3fp ;

        double y=(double)SHR3fp ;

    #elif TRNG == TMWC

        double x=(double)MWCfp ;

        double y=(double)MWCfp ;

    #elif TRNG == TKISS

        double x=(double)KISSfp ;

        double y=(double)KISSfp ;

    #endif

#endif

      ulong inside=((x*x+y*y) <= THEONE) ? 1:0;

      total+=inside;

   }

   return(total);

}

__global__ void MainLoopBlocks(ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,blockIdx.x);

   s[blockIdx.x]=total;

   __syncthreads();

}

__global__ void MainLoopThreads(ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,threadIdx.x);

   s[threadIdx.x]=total;

   __syncthreads();

}

__global__ void MainLoopHybrid(ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,blockDim.x*blockIdx.x+threadIdx.x);

   s[blockDim.x*blockIdx.x+threadIdx.x]=total;

   __syncthreads();

}

"""

    return(KERNEL_CODE_CUDA)

def KernelCodeOpenCL():

    KERNEL_CODE_OPENCL="""

#define TCONG 0

#define TSHR3 1

#define TMWC 2

#define TKISS 3

#define TINT32 0

#define TINT64 1

#define TFP32 2

#define TFP64 3

// Marsaglia RNG very simple implementation

#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 MWCfp MWC * 2.328306435454494e-10f

#define KISSfp KISS * 2.328306435454494e-10f

#define CONGfp CONG * 2.328306435454494e-10f

#define SHR3fp SHR3 * 2.328306435454494e-10f

ulong MainLoop(ulong iterations,uint seed_z,uint seed_w,size_t work)

{

#if TRNG == TCONG

   uint jcong=seed_z+work;

#elif TRNG == TSHR3

   uint jsr=seed_w+work;

#elif TRNG == TMWC

   uint z=seed_z+work;

   uint w=seed_w+work;

#elif TRNG == TKISS

   uint jcong=seed_z+work;

   uint jsr=seed_w+work;

   uint z=seed_z-work;

   uint w=seed_w-work;

#endif

   ulong total=0;

   for (ulong i=0;i<iterations;i++) {

#if TYPE == TINT32

    #define THEONE 1073741824

    #if TRNG == TCONG

        uint x=CONG>>17 ;

        uint y=CONG>>17 ;

    #elif TRNG == TSHR3

        uint x=SHR3>>17 ;

        uint y=SHR3>>17 ;

    #elif TRNG == TMWC

        uint x=MWC>>17 ;

        uint y=MWC>>17 ;

    #elif TRNG == TKISS

        uint x=KISS>>17 ;

        uint y=KISS>>17 ;

    #endif

#elif TYPE == TINT64

    #define THEONE 4611686018427387904

    #if TRNG == TCONG

        ulong x=(ulong)(CONG>>1) ;

        ulong y=(ulong)(CONG>>1) ;

    #elif TRNG == TSHR3

        ulong x=(ulong)(SHR3>>1) ;

        ulong y=(ulong)(SHR3>>1) ;

    #elif TRNG == TMWC

        ulong x=(ulong)(MWC>>1) ;

        ulong y=(ulong)(MWC>>1) ;

    #elif TRNG == TKISS

        ulong x=(ulong)(KISS>>1) ;

        ulong y=(ulong)(KISS>>1) ;

    #endif

#elif TYPE == TFP32

    #define THEONE 1.0f

    #if TRNG == TCONG

        float x=CONGfp ;

        float y=CONGfp ;

    #elif TRNG == TSHR3

        float x=SHR3fp ;

        float y=SHR3fp ;

    #elif TRNG == TMWC

        float x=MWCfp ;

        float y=MWCfp ;

    #elif TRNG == TKISS

      float x=KISSfp ;

      float y=KISSfp ;

    #endif

#elif TYPE == TFP64

#pragma OPENCL EXTENSION cl_khr_fp64: enable

    #define THEONE 1.0f

    #if TRNG == TCONG

        double x=(double)CONGfp ;

        double y=(double)CONGfp ;

    #elif TRNG == TSHR3

        double x=(double)SHR3fp ;

        double y=(double)SHR3fp ;

    #elif TRNG == TMWC

        double x=(double)MWCfp ;

        double y=(double)MWCfp ;

    #elif TRNG == TKISS

        double x=(double)KISSfp ;

        double y=(double)KISSfp ;

    #endif

#endif

      ulong inside=((x*x+y*y) <= THEONE) ? 1:0;

      total+=inside;

   }

   return(total);

}

__kernel void MainLoopGlobal(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,get_global_id(0));

   barrier(CLK_GLOBAL_MEM_FENCE);

   s[get_global_id(0)]=total;     

}

__kernel void MainLoopLocal(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,get_local_id(0));

   barrier(CLK_LOCAL_MEM_FENCE);

   s[get_local_id(0)]=total;

}

__kernel void MainLoopHybrid(__global ulong *s,ulong iterations,uint seed_w,uint seed_z)

{

   ulong total=MainLoop(iterations,seed_z,seed_w,get_global_id(0));

   barrier(CLK_GLOBAL_MEM_FENCE || CLK_LOCAL_MEM_FENCE);

   s[get_global_id(0)]=total;

}

"""

    return(KERNEL_CODE_OPENCL)

def MetropolisCuda(InputCU):

    print("Inside ",InputCU)

    iterations=InputCU['Iterations']

    steps=InputCU['Steps']

    blocks=InputCU['Blocks']

    threads=InputCU['Threads']

    Device=InputCU['Device']

    RNG=InputCU['RNG']

    ValueType=InputCU['ValueType']

    Marsaglia,Computing=DictionariesAPI()

    try:

        # For PyCUDA import

        import pycuda.driver as cuda

        from pycuda.compiler import SourceModule

        cuda.init()

        for Id in range(cuda.Device.count()):

            if Id==Device:

                XPU=cuda.Device(Id)

                print("GPU selected %s" % XPU.name())

        print

    except ImportError:

        print("Platform does not seem to support CUDA")

    circle=numpy.zeros(blocks*threads).astype(numpy.uint64)  

    circleCU = cuda.InOut(circle)

    #circleCU = cuda.mem_alloc(circle.size*circle.dtype.itemize)

    #cuda.memcpy_htod(circleCU, circle)

    Context=XPU.make_context()

    try:

        #mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DTRNG=%i -DTYPE=%s' % (Marsaglia[RNG],Computing[ValueType])])

        mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DTRNG=%i -DTYPE=%s' % (Marsaglia[RNG],Computing[ValueType])])

    except:

        print("Compilation seems to broke")

    MetropolisBlocksCU=mod.get_function("MainLoopBlocks")

    MetropolisThreadsCU=mod.get_function("MainLoopThreads")

    MetropolisHybridCU=mod.get_function("MainLoopHybrid")

    MyDuration=numpy.zeros(steps)

    jobs=blocks*threads;

    iterationsCU=numpy.uint64(iterations/jobs)

    if iterations%jobs!=0:

        iterationsCU+=numpy.uint64(1)

    for i in range(steps):

        start_time=time.time()

        try:

            MetropolisHybridCU(circleCU,

                               numpy.uint64(iterationsCU),

                               numpy.uint32(nprnd(2**32)),

                               numpy.uint32(nprnd(2**32)),

                               grid=(blocks,1),block=(threads,1,1))

        except:

            print("Crash during CUDA call")

        elapsed = time.time()-start_time

        print("(Blocks/Threads)=(%i,%i) method done in %.2f s..." % (blocks,threads,elapsed))

        MyDuration[i]=elapsed

    OutputCU={'Inside':sum(circle),'NewIterations':numpy.uint64(iterationsCU*jobs),'Duration':MyDuration}

    print(OutputCU)

    Context.pop()

    #Context.detach()

    return(OutputCU)

def MetropolisOpenCL(InputCL):

    import pyopencl as cl

    print("Inside ",InputCL)

    iterations=InputCL['Iterations']

    steps=InputCL['Steps']

    blocks=InputCL['Blocks']

    threads=InputCL['Threads']

    Device=InputCL['Device']

    RNG=InputCL['RNG']

    ValueType=InputCL['ValueType']

    Marsaglia,Computing=DictionariesAPI()

    # Initialisation des variables en les CASTant correctement

    Id=0

    HasXPU=False

    for platform in cl.get_platforms():

        for device in platform.get_devices():

            if Id==Device:

                XPU=device

                print("CPU/GPU selected: ",device.name.lstrip())

                HasXPU=True

            Id+=1

    if HasXPU==False:

        print("No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1))

        sys.exit()              

    # Je cree le contexte et la queue pour son execution

    try:

        ctx = cl.Context([XPU])

        queue = cl.CommandQueue(ctx,properties=cl.command_queue_properties.PROFILING_ENABLE)

    except:

        print("Crash during context creation")

    # Je recupere les flag possibles pour les buffers

    mf = cl.mem_flags

    circle=numpy.zeros(blocks*threads).astype(numpy.uint64)

    circleCL = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=circle)

    MetropolisCL = cl.Program(ctx,KernelCodeOpenCL()).build( options = "-cl-mad-enable -cl-fast-relaxed-math -DTRNG=%i -DTYPE=%s" % (Marsaglia[RNG],Computing[ValueType]))

    MyDuration=numpy.zeros(steps)

    jobs=blocks*threads;

    iterationsCL=numpy.uint64(iterations/jobs)

    if iterations%jobs!=0:

        iterationsCL+=1

    for i in range(steps):

        start_time=time.time()

        if threads == 1:

            CLLaunch=MetropolisCL.MainLoopGlobal(queue,(blocks,),None,

                                                 circleCL,

                                                 numpy.uint64(iterationsCL),

                                                 numpy.uint32(nprnd(2**32)),

                                                 numpy.uint32(nprnd(2**32)))

        else:

            CLLaunch=MetropolisCL.MainLoopHybrid(queue,(jobs,),(threads,),

                                                 circleCL,

                                                 numpy.uint64(iterationsCL),

                                                 numpy.uint32(nprnd(2**32)),

                                                 numpy.uint32(nprnd(2**32)))

        CLLaunch.wait()

        cl.enqueue_copy(queue, circle, circleCL).wait()

        elapsed = time.time()-start_time

        print("(Blocks/Threads)=(%i,%i) method done in %.2f s..." % (blocks,threads,elapsed))

        # Elapsed method based on CLLaunch doesn't work for Beignet OpenCL

        # elapsed = 1e-9*(CLLaunch.profile.end - CLLaunch.profile.start)

        # print circle,numpy.mean(circle),numpy.median(circle),numpy.std(circle)

        MyDuration[i]=elapsed

        # AllPi=4./numpy.float32(iterationsCL)*circle.astype(numpy.float32)

        # MyPi[i]=numpy.median(AllPi)

        # print MyPi[i],numpy.std(AllPi),MyDuration[i]

    circleCL.release()

    OutputCL={'Inside':sum(circle),'NewIterations':numpy.uint64(iterationsCL*jobs),'Duration':MyDuration}

    print(OutputCL)

    return(OutputCL)

def FitAndPrint(N,D,Curves):

    from scipy.optimize import curve_fit

    import matplotlib.pyplot as plt

    try:

        coeffs_Amdahl, matcov_Amdahl = curve_fit(Amdahl, N, D)

        D_Amdahl=Amdahl(N,coeffs_Amdahl[0],coeffs_Amdahl[1],coeffs_Amdahl[2])

        coeffs_Amdahl[1]=coeffs_Amdahl[1]*coeffs_Amdahl[0]/D[0]

        coeffs_Amdahl[2]=coeffs_Amdahl[2]*coeffs_Amdahl[0]/D[0]

        coeffs_Amdahl[0]=D[0]

        print("Amdahl Normalized: T=%.2f(%.6f+%.6f/N)" % (coeffs_Amdahl[0],coeffs_Amdahl[1],coeffs_Amdahl[2]))

    except:

        print("Impossible to fit for Amdahl law : only %i elements" % len(D))

    try:

        coeffs_AmdahlR, matcov_AmdahlR = curve_fit(AmdahlR, N, D)

        D_AmdahlR=AmdahlR(N,coeffs_AmdahlR[0],coeffs_AmdahlR[1])

        coeffs_AmdahlR[1]=coeffs_AmdahlR[1]*coeffs_AmdahlR[0]/D[0]

        coeffs_AmdahlR[0]=D[0]

        print("Amdahl Reduced Normalized: T=%.2f(%.6f+%.6f/N)" % (coeffs_AmdahlR[0],1-coeffs_AmdahlR[1],coeffs_AmdahlR[1]))

    except:

        print("Impossible to fit for Reduced Amdahl law : only %i elements" % len(D))

    try:

        coeffs_Mylq, matcov_Mylq = curve_fit(Mylq, N, D)

        coeffs_Mylq[1]=coeffs_Mylq[1]*coeffs_Mylq[0]/D[0]

        # coeffs_Mylq[2]=coeffs_Mylq[2]*coeffs_Mylq[0]/D[0]

        coeffs_Mylq[3]=coeffs_Mylq[3]*coeffs_Mylq[0]/D[0]

        coeffs_Mylq[0]=D[0]

        print("Mylq Normalized : T=%.2f(%.6f+%.6f/N)+%.6f*N" % (coeffs_Mylq[0],

                                                                coeffs_Mylq[1],

                                                                coeffs_Mylq[3],

                                                                coeffs_Mylq[2]))

        D_Mylq=Mylq(N,coeffs_Mylq[0],coeffs_Mylq[1],coeffs_Mylq[2],

                    coeffs_Mylq[3])

    except:

        print("Impossible to fit for Mylq law : only %i elements" % len(D))

    try:

        coeffs_Mylq2, matcov_Mylq2 = curve_fit(Mylq2, N, D)

        coeffs_Mylq2[1]=coeffs_Mylq2[1]*coeffs_Mylq2[0]/D[0]

        # coeffs_Mylq2[2]=coeffs_Mylq2[2]*coeffs_Mylq2[0]/D[0]

        # coeffs_Mylq2[3]=coeffs_Mylq2[3]*coeffs_Mylq2[0]/D[0]

        coeffs_Mylq2[4]=coeffs_Mylq2[4]*coeffs_Mylq2[0]/D[0]

        coeffs_Mylq2[0]=D[0]

        print("Mylq 2nd order Normalized: T=%.2f(%.6f+%.6f/N)+%.6f*N+%.6f*N^2" % (coeffs_Mylq2[0],coeffs_Mylq2[1],coeffs_Mylq2[4],coeffs_Mylq2[2],coeffs_Mylq2[3]))

    except:

        print("Impossible to fit for 2nd order Mylq law : only %i elements" % len(D))

    if Curves:

        plt.xlabel("Number of Threads/work Items")

        plt.ylabel("Total Elapsed Time")

        Experience,=plt.plot(N,D,'ro') 

    try:

        pAmdahl,=plt.plot(N,D_Amdahl,label="Loi de Amdahl")    

        pMylq,=plt.plot(N,D_Mylq,label="Loi de Mylq")

    except:

        print("Fit curves seem not to be available")

    plt.legend()

    plt.show()

if __name__=='__main__':

    # Set defaults values

    # Id of Device : 1 is for first find !

    Device=1

    # GPU style can be Cuda (Nvidia implementation) or OpenCL

    GpuStyle='OpenCL'

    # Iterations is integer

    Iterations=10000000

    # BlocksBlocks in first number of Blocks to explore

    BlocksBegin=1

    # BlocksEnd is last number of Blocks to explore

    BlocksEnd=1

    # BlocksStep is the step of Blocks to explore

    BlocksStep=1

    # ThreadsBlocks in first number of Blocks to explore

    ThreadsBegin=1

    # ThreadsEnd is last number of Blocks to explore

    ThreadsEnd=1

    # ThreadsStep is the step of Blocks to explore

    ThreadsStep=1

    # Redo is the times to redo the test to improve metrology

    Redo=1

    # OutMetrology is method for duration estimation : False is GPU inside

    OutMetrology=False

    Metrology='InMetro'

    # Curves is True to print the curves

    Curves=False

    # Fit is True to print the curves

    Fit=False

    # Marsaglia RNG

    RNG='MWC'

    # Value type : INT32, INT64, FP32, FP64

    ValueType='FP32'

    HowToUse='%s -o (Out of Core Metrology) -c (Print Curves) -d <DeviceId> -g <CUDA/OpenCL> -i <Iterations> -b <BlocksBegin> -e <BlocksEnd> -s <BlocksStep> -f <ThreadsFirst> -l <ThreadsLast> -t <ThreadssTep> -r <RedoToImproveStats> -m <SHR3/CONG/MWC/KISS> -v <INT32/INT64/FP32/FP64>'

    try:

        opts, args = getopt.getopt(sys.argv[1:],"hocg:i:b:e:s:f:l:t:r:d:m:v:",["gpustyle=","iterations=","blocksBegin=","blocksEnd=","blocksStep=","threadsFirst=","threadsLast=","threadssTep=","redo=","device=","marsaglia=","valuetype="])

    except getopt.GetoptError:

        print(HowToUse % sys.argv[0])

        sys.exit(2)

    # List of Devices

    Devices=[]

    Alu={}

    for opt, arg in opts:

        if opt == '-h':

            print(HowToUse % sys.argv[0])

            print("\nInformations about devices detected under OpenCL API:")

            # For PyOpenCL import

            try:

                import pyopencl as cl

                Id=0

                for platform in cl.get_platforms():

                    for device in platform.get_devices():

                        #deviceType=cl.device_type.to_string(device.type)

                        deviceType="*PU"

                        print("Device #%i from %s of type %s : %s" % (Id,platform.vendor.lstrip(),deviceType,device.name.lstrip()))

                        Id=Id+1

            except:

                print("Your platform does not seem to support OpenCL")

            print("\nInformations about devices detected under CUDA API:")

            # For PyCUDA import

            try:

                import pycuda.driver as cuda

                cuda.init()

                for Id in range(cuda.Device.count()):

                    device=cuda.Device(Id)

                    print("Device #%i of type GPU : %s" % (Id,device.name()))

                print

            except:

                print("Your platform does not seem to support CUDA")

            sys.exit()

        elif opt == '-o':

            OutMetrology=True

            Metrology='OutMetro'

        elif opt == '-c':

            Curves=True

        elif opt in ("-d", "--device"):

            Devices.append(int(arg))

        elif opt in ("-g", "--gpustyle"):

            GpuStyle = arg

        elif opt in ("-m", "--marsaglia"):

            RNG = arg

        elif opt in ("-v", "--valuetype"):

            ValueType = arg

        elif opt in ("-i", "--iterations"):

            Iterations = numpy.uint64(arg)

        elif opt in ("-b", "--blocksbegin"):

            BlocksBegin = int(arg)

        elif opt in ("-e", "--blocksend"):

            BlocksEnd = int(arg)

        elif opt in ("-s", "--blocksstep"):

            BlocksStep = int(arg)

        elif opt in ("-f", "--threadsfirst"):

            ThreadsBegin = int(arg)

        elif opt in ("-l", "--threadslast"):

            ThreadsEnd = int(arg)

        elif opt in ("-t", "--threadsstep"):

            ThreadsStep = int(arg)

        elif opt in ("-r", "--redo"):

            Redo = int(arg)

    print("Devices Identification : %s" % Devices)

    print("GpuStyle used : %s" % GpuStyle)

    print("Iterations : %s" % Iterations)

    print("Number of Blocks on begin : %s" % BlocksBegin)

    print("Number of Blocks on end : %s" % BlocksEnd)

    print("Step on Blocks : %s" % BlocksStep)

    print("Number of Threads on begin : %s" % ThreadsBegin)

    print("Number of Threads on end : %s" % ThreadsEnd)

    print("Step on Threads : %s" % ThreadsStep)

    print("Number of redo : %s" % Redo)

    print("Metrology done out of XPU : %r" % OutMetrology)

    print("Type of Marsaglia RNG used : %s" % RNG)

    print("Type of variable : %s" % ValueType)

    if GpuStyle=='CUDA':

        try:

            # For PyCUDA import

            import pycuda.driver as cuda

            cuda.init()

            for Id in range(cuda.Device.count()):

                device=cuda.Device(Id)

                print("Device #%i of type GPU : %s" % (Id,device.name()))

                if Id in Devices:

                    Alu[Id]='GPU'

        except ImportError:

            print("Platform does not seem to support CUDA")

    if GpuStyle=='OpenCL':

        try:

            # For PyOpenCL import

            import pyopencl as cl

            Id=0

            for platform in cl.get_platforms():

                for device in platform.get_devices():

                    #deviceType=cl.device_type.to_string(device.type)

                    deviceType="xPU"

                    print("Device #%i from %s of type %s : %s" % (Id,platform.vendor.lstrip().rstrip(),deviceType,device.name.lstrip().rstrip()))

                    if Id in Devices:

                    # Set the Alu as detected Device Type

                        Alu[Id]=deviceType

                    Id=Id+1

        except ImportError:

            print("Platform does not seem to support OpenCL")

    print(Devices,Alu)

    BlocksList=range(BlocksBegin,BlocksEnd+BlocksStep,BlocksStep)

    ThreadsList=range(ThreadsBegin,ThreadsEnd+ThreadsStep,ThreadsStep)

    ExploredJobs=numpy.array([]).astype(numpy.uint32)

    ExploredBlocks=numpy.array([]).astype(numpy.uint32)

    ExploredThreads=numpy.array([]).astype(numpy.uint32)

    avgD=numpy.array([]).astype(numpy.float32)

    medD=numpy.array([]).astype(numpy.float32)

    stdD=numpy.array([]).astype(numpy.float32)

    minD=numpy.array([]).astype(numpy.float32)

    maxD=numpy.array([]).astype(numpy.float32)

    avgR=numpy.array([]).astype(numpy.float32)

    medR=numpy.array([]).astype(numpy.float32)

    stdR=numpy.array([]).astype(numpy.float32)

    minR=numpy.array([]).astype(numpy.float32)

    maxR=numpy.array([]).astype(numpy.float32)

    for Blocks,Threads in itertools.product(BlocksList,ThreadsList):

        # print Blocks,Threads

        circle=numpy.zeros(Blocks*Threads).astype(numpy.uint64)

        ExploredJobs=numpy.append(ExploredJobs,Blocks*Threads)

        ExploredBlocks=numpy.append(ExploredBlocks,Blocks)

        ExploredThreads=numpy.append(ExploredThreads,Threads)

        if OutMetrology: 

            DurationItem=numpy.array([]).astype(numpy.float32)

            Duration=numpy.array([]).astype(numpy.float32)

            Rate=numpy.array([]).astype(numpy.float32)

            for i in range(Redo):

                start=time.time()

                if GpuStyle=='CUDA':

                    try:

                        InputCU={}

                        InputCU['Iterations']=Iterations

                        InputCU['Steps']=1

                        InputCU['Blocks']=Blocks

                        InputCU['Threads']=Threads

                        InputCU['Device']=Devices[0]

                        InputCU['RNG']=RNG

                        InputCU['ValueType']=ValueType                    

                        OutputCU=MetropolisCuda(InputCU)

                        Inside=OutputCU['Circle']

                        NewIterations=OutputCU['NewIterations']

                        Duration=OutputCU['Duration']

                    except:

                        print("Problem with (%i,%i) // computations on Cuda" % (Blocks,Threads))

                elif GpuStyle=='OpenCL':

                    try:

                        InputCL={}

                        InputCL['Iterations']=Iterations

                        InputCL['Steps']=1

                        InputCL['Blocks']=Blocks

                        InputCL['Threads']=Threads

                        InputCL['Device']=Devices[0]

                        InputCL['RNG']=RNG

                        InputCL['ValueType']=ValueType                    

                        OutputCL=MetropolisOpenCL(InputCL)

                        Inside=OutputCL['Circle']

                        NewIterations=OutputCL['NewIterations']

                        Duration=OutputCL['Duration']

                    except:

                        print("Problem with (%i,%i) // computations on OpenCL" % (Blocks,Threads))

                Duration=numpy.append(Duration,time.time()-start)

                Rate=numpy.append(Rate,NewIterations/Duration[-1])

        else:

            if GpuStyle=='CUDA':

                try:

                    InputCU={}

                    InputCU['Iterations']=Iterations

                    InputCU['Steps']=Redo

                    InputCU['Blocks']=Blocks

                    InputCU['Threads']=Threads

                    InputCU['Device']=Devices[0]

                    InputCU['RNG']=RNG

                    InputCU['ValueType']=ValueType

                    OutputCU=MetropolisCuda(InputCU)

                    Inside=OutputCU['Inside']

                    NewIterations=OutputCU['NewIterations']

                    Duration=OutputCU['Duration']

                except:

                    print("Problem with (%i,%i) // computations on Cuda" % (Blocks,Threads))

                try:

                    pycuda.context.pop()

                except:

                    pass

            elif GpuStyle=='OpenCL':

                try:

                    InputCL={}

                    InputCL['Iterations']=Iterations

                    InputCL['Steps']=Redo

                    InputCL['Blocks']=Blocks

                    InputCL['Threads']=Threads

                    InputCL['Device']=Devices[0]

                    InputCL['RNG']=RNG

                    InputCL['ValueType']=ValueType

                    OutputCL=MetropolisOpenCL(InputCL)

                    Inside=OutputCL['Inside']

                    NewIterations=OutputCL['NewIterations']

                    Duration=OutputCL['Duration']

                except:

                    print("Problem with (%i,%i) // computations on OpenCL" % (Blocks,Threads))

            Rate=NewIterations/Duration

            print("Pi estimation %.8f" % (4./NewIterations*Inside))

        avgD=numpy.append(avgD,numpy.average(Duration))

        medD=numpy.append(medD,numpy.median(Duration))

        stdD=numpy.append(stdD,numpy.std(Duration))

        minD=numpy.append(minD,numpy.min(Duration))

        maxD=numpy.append(maxD,numpy.max(Duration))

        avgR=numpy.append(avgR,numpy.average(Rate))

        medR=numpy.append(medR,numpy.median(Rate))

        stdR=numpy.append(stdR,numpy.std(Rate))

        minR=numpy.append(minR,numpy.min(Rate))

        maxR=numpy.append(maxR,numpy.max(Rate))

        print("%.2f %.2f %.2f %.2f %.2f %i %i %i %i %i" % (avgD[-1],medD[-1],stdD[-1],minD[-1],maxD[-1],avgR[-1],medR[-1],stdR[-1],minR[-1],maxR[-1]))

        numpy.savez("Pi_%s_%s_%s_%s_%s_%s_%s_%s_%.8i_Device%i_%s_%s" % (ValueType,RNG,Alu[Devices[0]],GpuStyle,BlocksBegin,BlocksEnd,ThreadsBegin,ThreadsEnd,Iterations,Devices[0],Metrology,gethostname()),(ExploredBlocks,ExploredThreads,avgD,medD,stdD,minD,maxD,avgR,medR,stdR,minR,maxR))

        ToSave=[ ExploredBlocks,ExploredThreads,avgD,medD,stdD,minD,maxD,avgR,medR,stdR,minR,maxR ]

        numpy.savetxt("Pi_%s_%s_%s_%s_%s_%s_%s_%i_%.8i_Device%i_%s_%s" % (ValueType,RNG,Alu[Devices[0]],GpuStyle,BlocksBegin,BlocksEnd,ThreadsBegin,ThreadsEnd,Iterations,Devices[0],Metrology,gethostname()),numpy.transpose(ToSave),fmt='%i %i %e %e %e %e %e %i %i %i %i %i')

    if Fit:

        FitAndPrint(ExploredJobs,median,Curves)