Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python

# Version using Multiprocessing module

#

# Pi-by-MC

#

# CC BY-NC-SA 2013 : <emmanuel.quemener@ens-lyon.fr> 

#

import getopt

import numpy

from random import random

# Multithread library call

from multiprocessing import Pool

# Common tools

import sys

import getopt

import time

import matplotlib.pyplot as plt

import math

from scipy.optimize import curve_fit

from socket import gethostname

# 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+c*N+p/N))

# Predicted Mylq Law with second order

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

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

def MainLoop(iterations):

    total=0

    # Absulute necessary to use xrange instead of range to avoid out of memory !

    for i in xrange(iterations):

        # Random access coordonate

        x,y=random(),random()

        if ((x*x+y*y) < 1.0):

            total+=1

    return(total)

def MetropolisMP(circle,iterations,steps,jobs):

    MyPi=numpy.zeros(steps)

    MyDuration=numpy.array([])

    # Define iterations to send to each node

    if iterations%jobs==0:

        iterationsMP=iterations/jobs

    else:

        iterationsMP=iterations/jobs+1

    print "%i iterations will be send to each core" % iterationsMP

    WorkItems=[]

    for i in xrange(0,jobs):

        WorkItems.append(iterationsMP)

    for i in xrange(steps):

        start=time.time()

        pool=Pool(processes=jobs)

        # Define the number of processes to be launched at a time

        # pool=Pool(processes=CORES)

        print "Start on %i processors..." % jobs

        # MAP: distribution of usecases T to be applied to MetropolisStrip 

        # WorkLaunches=[ pool.apply_async(MainLoop,(wi,)) for wi in WorkItems ]

        # circle=[wl.get() for wl in WorkLaunches]

        circle=pool.map(MainLoop,WorkItems)

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

        MyPi[i]=4.*float(numpy.sum(circle))/float(iterationsMP)/float(jobs)

    return(numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration))

def FitAndPrint(N,D,Curves):

    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[2],

                                                                coeffs_Mylq[3])

        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) 

    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

    # Iterations is integer

    Iterations=1000000

    # JobStart in first number of Jobs to explore

    JobStart=1

    # JobEnd is last number of Jobs to explore

    JobEnd=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

    # Curves is True to print the curves

    Curves=False

    try:

        opts, args = getopt.getopt(sys.argv[1:],"hoci:s:e:r:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","redo=","device="])

    except getopt.GetoptError:

        print '%s -o (Out of Core Metrology) -c (Print Curves) -i <Iterations> -s <JobStart> -e <JobEnd> -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) -i <Iterations> -s <JobStart> -e <JobEnd> -r <RedoToImproveStats>' % sys.argv[0]

            sys.exit()

        elif opt == '-o':

            OutMetrology=True

        elif opt == '-c':

            Curves=True

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

            Iterations = numpy.uint32(arg)

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

            JobStart = int(arg)

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

            JobEnd = int(arg)

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

            Redo = int(arg)

    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 "Metrology done out : %r" % OutMetrology

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

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

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

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

    Jobs=JobStart

    while Jobs <= JobEnd:

        avg,med,std=0,0,0

        ExploredJobs=numpy.append(ExploredJobs,Jobs)

        circle=numpy.zeros(Jobs).astype(numpy.uint32)

        if OutMetrology:

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

            for i in range(Redo):

                start=time.time()

                MetropolisMP(circle,Iterations,Redo,Jobs)

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

            avg=numpy.mean(duration)

            med=numpy.median(duration)

            std=numpy.std(duration)

        else:

            avg,med,std=MetropolisMP(circle,Iterations,Redo,Jobs)

        if (avg,med,std) != (0,0,0):

            print "avg,med,std",avg,med,std

            average=numpy.append(average,avg)

            median=numpy.append(median,med)

            stddev=numpy.append(stddev,std)

        else:

            print "Values seem to be wrong..."

        # THREADS*=2

        numpy.savez("Pi_%s_%i_%.8i_%s" % (JobStart,JobEnd,Iterations,gethostname()),(ExploredJobs,average,median,stddev))

        Jobs+=1

    FitAndPrint(ExploredJobs,median,Curves)

#!/usr/bin/env python

# If numpy is loaded without being used, it does not work !

from math import exp

from random import random

# Multithread library call

#import multiprocessing

from multiprocessing import Pool

import numpy

from numpy.random import rand

def MainLoop(iterations):

    total=numpy.int(0)

    for i in range(0,iterations):

        # Random access coordonate

        #x,y=random(),random()

        x,y=rand(),rand()

        if ((x*x+y*y) < 1.0):

                total+=1

    return(total)

if __name__=='__main__':

    # Machine

    CORES=4

    # A=numpy.zeros(10)

    # Au dessus de 4 10^7, le MPI ne veut plus se lancer...

    Iterations=100000000

    total=0

    # Define iterations to send to each node

    if Iterations%CORES==0:

        iterations=Iterations/CORES

    else:

        iterations=Iterations/CORES+1

        print "%i iterations will be send to each core" % iterations

    IT=[]

    for i in range(0,CORES):

        IT.append(iterations)

    print IT

    # Define the number of processes to be launched at a time

    # POOL: resources to be used ( 

    # pool=Pool(processes=CORES)

    pool=Pool(processes=CORES)

    print "Start on %i processors..." % CORES

    # MAP: distribution of usecases T to be applied to MetropolisStrip 

    Results=pool.map(MainLoop,IT)

    # print sum(Results.get())

    print Results

import numpy

npzfile=numpy.load("Ising2D_MP_256_16777216.npz")

npzfile['arr_0']

T=npzfile['arr_0'].transpose()[0]

E=npzfile['arr_0'].transpose()[1]

M=npzfile['arr_0'].transpose()[2]

print T

print E

print M

import matplotlib.pyplot 

matplotlib.pyplot.plot(T,E)

matplotlib.pyplot.show()

#!/usr/bin/env python

# Version using Multiprocessing module

#

# Pi-by-MC

#

# CC BY-NC-SA 2013 : <emmanuel.quemener@ens-lyon.fr> 

#

import getopt

import numpy

from random import random

# Multithread library call

from multiprocessing import Pool

# Common tools

import sys

import getopt

import time

import matplotlib.pyplot as plt

import math

from scipy.optimize import curve_fit

from socket import gethostname

# 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+c*N+p/N))

# Predicted Mylq Law with second order

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

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

def MainLoop(iterations):

    total=0

    # Absulute necessary to use xrange instead of range to avoid out of memory !

    for i in xrange(iterations):

        # Random access coordonate

        x,y=random(),random()

        if ((x*x+y*y) < 1.0):

            total+=1

    return(total)

def MetropolisMP(circle,iterations,steps,jobs):

    MyPi=numpy.zeros(steps)

    MyDuration=numpy.array([])

    # Define iterations to send to each node

    if iterations%jobs==0:

        iterationsMP=iterations/jobs

    else:

        iterationsMP=iterations/jobs+1

    print "%i iterations will be send to each core" % iterationsMP

    WorkItems=[]

    for i in xrange(0,jobs):

        WorkItems.append(iterationsMP)

    for i in xrange(steps):

        start=time.time()

        pool=Pool(processes=jobs)

        # Define the number of processes to be launched at a time

        # pool=Pool(processes=CORES)

        print "Start on %i processors..." % jobs

        # MAP: distribution of usecases T to be applied to MetropolisStrip 

        # WorkLaunches=[ pool.apply_async(MainLoop,(wi,)) for wi in WorkItems ]

        # circle=[wl.get() for wl in WorkLaunches]

        circle=pool.map(MainLoop,WorkItems)

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

        MyPi[i]=4.*float(numpy.sum(circle))/float(iterationsMP)/float(jobs)

    return(numpy.mean(MyDuration),numpy.median(MyDuration),numpy.std(MyDuration))

def FitAndPrint(N,D,Curves):

    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[2],

                                                                coeffs_Mylq[3])

        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) 

    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

    # Iterations is integer

    Iterations=1000000

    # JobStart in first number of Jobs to explore

    JobStart=1

    # JobEnd is last number of Jobs to explore

    JobEnd=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

    # Curves is True to print the curves

    Curves=False

    try:

        opts, args = getopt.getopt(sys.argv[1:],"hoci:s:e:r:",["alu=","gpustyle=","parastyle=","iterations=","jobstart=","jobend=","redo=","device="])

    except getopt.GetoptError:

        print '%s -o (Out of Core Metrology) -c (Print Curves) -i <Iterations> -s <JobStart> -e <JobEnd> -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) -i <Iterations> -s <JobStart> -e <JobEnd> -r <RedoToImproveStats>' % sys.argv[0]

            sys.exit()

        elif opt == '-o':

            OutMetrology=True

        elif opt == '-c':

            Curves=True

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

            Iterations = numpy.uint32(arg)

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

            JobStart = int(arg)

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

            JobEnd = int(arg)

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

            Redo = int(arg)

    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 "Metrology done out : %r" % OutMetrology

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

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

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

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

    Jobs=JobStart

    while Jobs <= JobEnd:

        avg,med,std=0,0,0

        ExploredJobs=numpy.append(ExploredJobs,Jobs)

        circle=numpy.zeros(Jobs).astype(numpy.uint32)

        if OutMetrology:

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

            for i in range(Redo):

                start=time.time()

                try:

                    MetropolisMP(circle,Iterations,Redo,Jobs)

                except:

                    print "Problem with %i // computations" % Jobs

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

            avg=numpy.mean(duration)

            med=numpy.median(duration)

            std=numpy.std(duration)

        else:

            try:

                avg,med,std=MetropolisMP(circle,Iterations,Redo,Jobs)

            except:

                print "Problem with %i // computations" % Jobs

        if (avg,med,std) != (0,0,0):

            print "avg,med,std",avg,med,std

            average=numpy.append(average,avg)

            median=numpy.append(median,med)

            stddev=numpy.append(stddev,std)

        else:

            print "Values seem to be wrong..."

        # THREADS*=2

        numpy.savez("Pi_%s_%i_%.8i_%s" % (JobStart,JobEnd,Iterations,gethostname()),(ExploredJobs,average,median,stddev))

        Jobs+=1

    FitAndPrint(ExploredJobs,median,Curves)

#!/usr/bin/env python

# If numpy is loaded without being used, it does not work !

from math import exp

from random import random

# Multithread library call

from multiprocessing import Pool

def MainLoop(iterations):

    total=0

    for i in xrange(iterations):

        # Random access coordonate

        x,y=random(),random()

        if ((x*x+y*y) < 1.0):

                total+=1

    return(total)

if __name__=='__main__':

   # Machine

   CORES=4

   # Au dessus de 4 10^7, le MPI ne veut plus se lancer...

   Iterations=1000000000

   total=0

   # Define iterations to send to each node

   if Iterations%CORES==0:

	   iterations=Iterations/CORES

   else:

	   iterations=Iterations/CORES+1

   print "%i iterations will be send to each core" % iterations

   IT=[]

   for i in range(0,CORES):

	   IT.append(iterations)

   print IT

   # Define the number of processes to be launched at a time

   # POOL: resources to be used ( 

   #pool=Pool(processes=CORES)

   pool=Pool(CORES)

   print "Start on %i processors..." % CORES

   # MAP: distribution of usecases T to be applied to MetropolisStrip 

   Results=[ pool.apply_async(MainLoop,(i,)) for i in IT]

   results=[r.get() for r in Results]

   print results,sum(results)

import numpy

import sys

npzfile=numpy.load(sys.argv[1])

T=npzfile['arr_0'][0]

E=npzfile['arr_0'][1]

ED=npzfile['arr_0'][2]

import matplotlib.pyplot

matplotlib.pyplot.plot(T,E)

#matplotlib.pyplot.plot(T,E,T,ED)

matplotlib.pyplot.show()

import numpy

import sys

npzSerial=numpy.load("Ising2D_Serial_256_16777216.npz")

npzOCLG=numpy.load("Ising2D_OCLG_256_16777216.npz")

npzOCLL=numpy.load("Ising2D_OCLL_256_16777216.npz")

npzCUDAG=numpy.load("Ising2D_CUDAG_256_16777216.npz")

npzCUDAL=numpy.load("Ising2D_CUDAL_256_16777216.npz")

T=npzSerial['arr_0'][0]

ES=npzSerial['arr_0'][1]

EOCLG=npzOCLG['arr_0'][1]

EOCLL=npzOCLL['arr_0'][1]

ECUDAG=npzCUDAG['arr_0'][1]

ECUDAL=npzCUDAL['arr_0'][1]

import matplotlib.pyplot

matplotlib.pyplot.plot(T,ES,T,EOCLG,T,EOCLL,T,ECUDAG,T,ECUDAL)

matplotlib.pyplot.show()

import numpy

import sys

npzSerial=numpy.load("Ising2D_Serial_256_16777216.npz")

npzMP=numpy.load("Ising2D_MP_256_16777216.npz")

npzMPI=numpy.load("Ising2D_MPI_256_16777216.npz")

T=npzSerial['arr_0'][0]

ES=npzSerial['arr_0'][1]

EMP=npzMP['arr_0'].transpose()[1]

EMPI=npzMPI['arr_0'].transpose()[1]

import matplotlib.pyplot

matplotlib.pyplot.plot(T,ES,T,EMP,T,EMPI)

matplotlib.pyplot.show()

import numpy

import sys

npzOCLLG=numpy.load("Ising2D_OCLL_256_16777216_GPU_atlas.npz")

npzOCLLC=numpy.load("Ising2D_OCLL_256_16777216_CPU_atlas.npz")

npzCUDAL=numpy.load("Ising2D_CUDAL_256_16777216_GPU_atlas.npz")

npzOCLLG2=numpy.load("Ising2D_OCLL_256_16777216_GPU_o745-2.npz")

npzCUDAL2=numpy.load("Ising2D_CUDAL_256_16777216_GPU_o745-2.npz")

T=npzOCLLG['arr_0'][0]

EOCLLG=npzOCLLG['arr_0'][1]

POCLLG=npzOCLLG['arr_0'][2]

EOCLLC=npzOCLLC['arr_0'][1]

POCLLC=npzOCLLC['arr_0'][2]

ECUDAL=npzCUDAL['arr_0'][1]

PCUDAL=npzCUDAL['arr_0'][2]

EOCLLG2=npzOCLLG2['arr_0'][1]

POCLLG2=npzOCLLG2['arr_0'][2]

ECUDAL2=npzCUDAL2['arr_0'][1]

PCUDAL2=npzCUDAL2['arr_0'][2]

dPOCLLG=numpy.diff(POCLLG)

dPOCLLC=numpy.diff(POCLLC)

dPCUDAL=numpy.diff(PCUDAL)

dPOCLLG2=numpy.diff(POCLLG2)

dPCUDAL2=numpy.diff(PCUDAL2)

dPOCLLG=numpy.insert(dPOCLLG,0,0)

dPOCLLC=numpy.insert(dPOCLLC,0,0)

dPCUDAL=numpy.insert(dPCUDAL,0,0)

dPOCLLG2=numpy.insert(dPOCLLG2,0,0)

dPCUDAL2=numpy.insert(dPCUDAL2,0,0)

from matplotlib.pyplot import *

plot(T,EOCLLG,T,EOCLLC,T,ECUDAL,T,dPOCLLG,T,dPOCLLC,T,dPCUDAL,

     T,EOCLLG2,T,ECUDAL2,T,dPOCLLG2,T,dPCUDAL2)

legend((r'OCL GPU GTX560',r'OCL CPU Xeon4c',r'Cuda GTX560',

        r'OCL GPU GT8400',r'Cuda GT8400'))

xlabel('Temperature',{'fontsize':20})

ylabel('Energy',{'fontsize':20})

matplotlib.pyplot.show()

import numpy

import sys

npzSerial=numpy.load("Ising2D_Serial_256_16777216.npz")

npzMP=numpy.load("Ising2D_MP_256_16777216.npz")

npzMPI=numpy.load("Ising2D_MPI_256_16777216.npz")

npzOCLG=numpy.load("Ising2D_OCLG_256_16777216.npz")

npzOCLL=numpy.load("Ising2D_OCLL_256_16777216.npz")

npzCUDAG=numpy.load("Ising2D_CUDAG_256_16777216.npz")

npzCUDAL=numpy.load("Ising2D_CUDAL_256_16777216.npz")

T=npzSerial['arr_0'][0]

ES=npzSerial['arr_0'][1]

EMP=npzMP['arr_0'].transpose()[1]

EMPI=npzMPI['arr_0'].transpose()[1]

EOCLG=npzOCLG['arr_0'][1]

EOCLL=npzOCLL['arr_0'][1]

ECUDAG=npzCUDAG['arr_0'][1]

ECUDAL=npzCUDAL['arr_0'][1]

import matplotlib.pyplot

matplotlib.pyplot.plot(T,ES,T,EMP,T,EMPI,T,EOCLG,T,EOCLL,T,ECUDAG,T,ECUDAL)

matplotlib.pyplot.show()

#include <Python.h>

#include <math.h>

// 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

int _InsideCircle(unsigned int iterations,

		  unsigned int seed_w,unsigned int seed_z)

{

   unsigned int z=seed_z;

   unsigned int w=seed_w;

   unsigned int i;

   int total=0;

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

      float x=MWCfp ;

      float y=MWCfp ;

      // Matching test

      int inside=((x*x+y*y) < 1.0f) ? 1:0;

      total+=inside;

   }

   return(total);

}

static PyObject* InsideCircle(PyObject* self, PyObject* args)

{

  int i,z,w;

  if (!PyArg_ParseTuple(args, "iii", &i,&z,&w))

    return NULL;

  return Py_BuildValue("i", _InsideCircle(i,z,w));

}

static PyMethodDef MonteCarloMethods[] = {

  {"InsideCircle", InsideCircle, METH_VARARGS, "Inside Circle Randomly Check"},

  {NULL, NULL, 0, NULL}

};

PyMODINIT_FUNC

initMonteCarlo(void)

{

  (void) Py_InitModule("MonteCarlo", MonteCarloMethods);

}

#!/usr/bin/env python

import MonteCarlo

print MonteCarlo.InsideCircle(1000000000,1,2)

from distutils.core import setup, Extension

module1 = Extension('MonteCarlo', sources = ['pimodule.c'])

setup (name = 'PackageName',

       version = '1.0',

       description = 'This is a demo package',

       ext_modules = [module1])

rm -r build MonteCarlo.so 

python setup.py build

cp build/lib.linux-x86_64-2.7/MonteCarlo.so .

python MC.py

/*hello.i*/

%module hello

extern void say_hello(const char* name);

# This file was automatically generated by SWIG (http://www.swig.org).

# Version 2.0.7

#

# Do not make changes to this file unless you know what you are doing--modify

# the SWIG interface file instead.

from sys import version_info

if version_info >= (2,6,0):

    def swig_import_helper():

        from os.path import dirname

        import imp

        fp = None

        try:

            fp, pathname, description = imp.find_module('_hello', [dirname(__file__)])

        except ImportError:

            import _hello

            return _hello

        if fp is not None:

            try:

                _mod = imp.load_module('_hello', fp, pathname, description)

            finally:

                fp.close()

            return _mod

    _hello = swig_import_helper()

    del swig_import_helper

else:

    import _hello

del version_info

try:

    _swig_property = property

except NameError:

    pass # Python < 2.2 doesn't have 'property'.

def _swig_setattr_nondynamic(self,class_type,name,value,static=1):

    if (name == "thisown"): return self.this.own(value)

    if (name == "this"):

        if type(value).__name__ == 'SwigPyObject':

            self.__dict__[name] = value

            return

    method = class_type.__swig_setmethods__.get(name,None)

    if method: return method(self,value)

    if (not static):

        self.__dict__[name] = value

    else:

        raise AttributeError("You cannot add attributes to %s" % self)

def _swig_setattr(self,class_type,name,value):

    return _swig_setattr_nondynamic(self,class_type,name,value,0)

def _swig_getattr(self,class_type,name):

    if (name == "thisown"): return self.this.own()

    method = class_type.__swig_getmethods__.get(name,None)

    if method: return method(self)

    raise AttributeError(name)

def _swig_repr(self):

    try: strthis = "proxy of " + self.this.__repr__()

    except: strthis = ""

    return "<%s.%s; %s >" % (self.__class__.__module__, self.__class__.__name__, strthis,)

try:

    _object = object

    _newclass = 1

except AttributeError:

    class _object : pass

    _newclass = 0

def say_hello(*args):

  return _hello.say_hello(*args)

say_hello = _hello.say_hello

# This file is compatible with both classic and new-style classes.

/*hellomodule.c*/

#include <stdio.h>

void say_hello(const char* name) {

    printf("Hello %s!\n", name);

}

/* ----------------------------------------------------------------------------

 * This file was automatically generated by SWIG (http://www.swig.org).

 * Version 2.0.7

 * 

 * This file is not intended to be easily readable and contains a number of 

 * coding conventions designed to improve portability and efficiency. Do not make

 * changes to this file unless you know what you are doing--modify the SWIG 

 * interface file instead. 

 * ----------------------------------------------------------------------------- */

#define SWIGPYTHON

#define SWIG_PYTHON_DIRECTOR_NO_VTABLE

/* -----------------------------------------------------------------------------

 *  This section contains generic SWIG labels for method/variable

 *  declarations/attributes, and other compiler dependent labels.

 * ----------------------------------------------------------------------------- */

/* template workaround for compilers that cannot correctly implement the C++ standard */

#ifndef SWIGTEMPLATEDISAMBIGUATOR

# if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x560)

#  define SWIGTEMPLATEDISAMBIGUATOR template

# elif defined(__HP_aCC)

/* Needed even with `aCC -AA' when `aCC -V' reports HP ANSI C++ B3910B A.03.55 */

/* If we find a maximum version that requires this, the test would be __HP_aCC <= 35500 for A.03.55 */

#  define SWIGTEMPLATEDISAMBIGUATOR template

# else

#  define SWIGTEMPLATEDISAMBIGUATOR

# endif

#endif

/* inline attribute */

#ifndef SWIGINLINE

# if defined(__cplusplus) || (defined(__GNUC__) && !defined(__STRICT_ANSI__))

#   define SWIGINLINE inline

# else

#   define SWIGINLINE

# endif

#endif

/* attribute recognised by some compilers to avoid 'unused' warnings */

#ifndef SWIGUNUSED

# if defined(__GNUC__)

#   if !(defined(__cplusplus)) || (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4))

#     define SWIGUNUSED __attribute__ ((__unused__)) 

#   else

#     define SWIGUNUSED

#   endif

# elif defined(__ICC)

#   define SWIGUNUSED __attribute__ ((__unused__)) 

# else

#   define SWIGUNUSED 

# endif

#endif

#ifndef SWIG_MSC_UNSUPPRESS_4505

# if defined(_MSC_VER)

#   pragma warning(disable : 4505) /* unreferenced local function has been removed */

# endif 

#endif

#ifndef SWIGUNUSEDPARM

# ifdef __cplusplus

#   define SWIGUNUSEDPARM(p)

# else

#   define SWIGUNUSEDPARM(p) p SWIGUNUSED 

# endif

#endif

/* internal SWIG method */

#ifndef SWIGINTERN

# define SWIGINTERN static SWIGUNUSED

#endif

/* internal inline SWIG method */

#ifndef SWIGINTERNINLINE

# define SWIGINTERNINLINE SWIGINTERN SWIGINLINE

#endif

/* exporting methods */

#if (__GNUC__ >= 4) || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)

#  ifndef GCC_HASCLASSVISIBILITY

#    define GCC_HASCLASSVISIBILITY

#  endif

#endif

#ifndef SWIGEXPORT

# if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__)

#   if defined(STATIC_LINKED)

#     define SWIGEXPORT

#   else

#     define SWIGEXPORT __declspec(dllexport)

#   endif

# else

#   if defined(__GNUC__) && defined(GCC_HASCLASSVISIBILITY)

#     define SWIGEXPORT __attribute__ ((visibility("default")))

#   else

#     define SWIGEXPORT

#   endif

# endif

#endif

/* calling conventions for Windows */

#ifndef SWIGSTDCALL

# if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__)

#   define SWIGSTDCALL __stdcall

# else

#   define SWIGSTDCALL

# endif 

#endif

/* Deal with Microsoft's attempt at deprecating C standard runtime functions */

#if !defined(SWIG_NO_CRT_SECURE_NO_DEPRECATE) && defined(_MSC_VER) && !defined(_CRT_SECURE_NO_DEPRECATE)

# define _CRT_SECURE_NO_DEPRECATE

#endif

/* Deal with Microsoft's attempt at deprecating methods in the standard C++ library */

#if !defined(SWIG_NO_SCL_SECURE_NO_DEPRECATE) && defined(_MSC_VER) && !defined(_SCL_SECURE_NO_DEPRECATE)

# define _SCL_SECURE_NO_DEPRECATE

#endif

/* Python.h has to appear first */

#include <Python.h>

/* -----------------------------------------------------------------------------

 * swigrun.swg

 *

 * This file contains generic C API SWIG runtime support for pointer

 * type checking.

 * ----------------------------------------------------------------------------- */

/* This should only be incremented when either the layout of swig_type_info changes,

   or for whatever reason, the runtime changes incompatibly */

#define SWIG_RUNTIME_VERSION "4"

/* define SWIG_TYPE_TABLE_NAME as "SWIG_TYPE_TABLE" */

#ifdef SWIG_TYPE_TABLE

# define SWIG_QUOTE_STRING(x) #x

# define SWIG_EXPAND_AND_QUOTE_STRING(x) SWIG_QUOTE_STRING(x)

# define SWIG_TYPE_TABLE_NAME SWIG_EXPAND_AND_QUOTE_STRING(SWIG_TYPE_TABLE)

#else

# define SWIG_TYPE_TABLE_NAME

#endif

/*

  You can use the SWIGRUNTIME and SWIGRUNTIMEINLINE macros for

  creating a static or dynamic library from the SWIG runtime code.

  In 99.9% of the cases, SWIG just needs to declare them as 'static'.

  But only do this if strictly necessary, ie, if you have problems

  with your compiler or suchlike.

*/

#ifndef SWIGRUNTIME

# define SWIGRUNTIME SWIGINTERN

#endif

#ifndef SWIGRUNTIMEINLINE

# define SWIGRUNTIMEINLINE SWIGRUNTIME SWIGINLINE

#endif

/*  Generic buffer size */

#ifndef SWIG_BUFFER_SIZE

# define SWIG_BUFFER_SIZE 1024

#endif

/* Flags for pointer conversions */

#define SWIG_POINTER_DISOWN        0x1

#define SWIG_CAST_NEW_MEMORY       0x2

/* Flags for new pointer objects */

#define SWIG_POINTER_OWN           0x1

/* 

   Flags/methods for returning states.

   The SWIG conversion methods, as ConvertPtr, return an integer 

   that tells if the conversion was successful or not. And if not,

   an error code can be returned (see swigerrors.swg for the codes).

   Use the following macros/flags to set or process the returning

   states.

   In old versions of SWIG, code such as the following was usually written:

     if (SWIG_ConvertPtr(obj,vptr,ty.flags) != -1) {

       // success code

     } else {

       //fail code

     }

   Now you can be more explicit:

    int res = SWIG_ConvertPtr(obj,vptr,ty.flags);

    if (SWIG_IsOK(res)) {

      // success code

    } else {

      // fail code

    }

   which is the same really, but now you can also do

    Type *ptr;

    int res = SWIG_ConvertPtr(obj,(void **)(&ptr),ty.flags);

    if (SWIG_IsOK(res)) {

      // success code

      if (SWIG_IsNewObj(res) {

        ...

	delete *ptr;

      } else {

        ...

      }

    } else {

      // fail code

    }

   I.e., now SWIG_ConvertPtr can return new objects and you can

   identify the case and take care of the deallocation. Of course that

   also requires SWIG_ConvertPtr to return new result values, such as

      int SWIG_ConvertPtr(obj, ptr,...) {         

        if (<obj is ok>) {			       

          if (<need new object>) {		       

            *ptr = <ptr to new allocated object>; 

            return SWIG_NEWOBJ;		       

          } else {				       

            *ptr = <ptr to old object>;	       

            return SWIG_OLDOBJ;		       

          } 				       

        } else {				       

          return SWIG_BADOBJ;		       

        }					       

      }

   Of course, returning the plain '0(success)/-1(fail)' still works, but you can be

   more explicit by returning SWIG_BADOBJ, SWIG_ERROR or any of the

   SWIG errors code.

   Finally, if the SWIG_CASTRANK_MODE is enabled, the result code

   allows to return the 'cast rank', for example, if you have this

       int food(double)

       int fooi(int);

   and you call

      food(1)   // cast rank '1'  (1 -> 1.0)

      fooi(1)   // cast rank '0'

   just use the SWIG_AddCast()/SWIG_CheckState()

*/

#define SWIG_OK                    (0) 

#define SWIG_ERROR                 (-1)

#define SWIG_IsOK(r)               (r >= 0)

#define SWIG_ArgError(r)           ((r != SWIG_ERROR) ? r : SWIG_TypeError)  

/* The CastRankLimit says how many bits are used for the cast rank */

#define SWIG_CASTRANKLIMIT         (1 << 8)

/* The NewMask denotes the object was created (using new/malloc) */

#define SWIG_NEWOBJMASK            (SWIG_CASTRANKLIMIT  << 1)

/* The TmpMask is for in/out typemaps that use temporal objects */

#define SWIG_TMPOBJMASK            (SWIG_NEWOBJMASK << 1)

/* Simple returning values */

#define SWIG_BADOBJ                (SWIG_ERROR)

#define SWIG_OLDOBJ                (SWIG_OK)

#define SWIG_NEWOBJ                (SWIG_OK | SWIG_NEWOBJMASK)

#define SWIG_TMPOBJ                (SWIG_OK | SWIG_TMPOBJMASK)

/* Check, add and del mask methods */

#define SWIG_AddNewMask(r)         (SWIG_IsOK(r) ? (r | SWIG_NEWOBJMASK) : r)

#define SWIG_DelNewMask(r)         (SWIG_IsOK(r) ? (r & ~SWIG_NEWOBJMASK) : r)

#define SWIG_IsNewObj(r)           (SWIG_IsOK(r) && (r & SWIG_NEWOBJMASK))

#define SWIG_AddTmpMask(r)         (SWIG_IsOK(r) ? (r | SWIG_TMPOBJMASK) : r)

#define SWIG_DelTmpMask(r)         (SWIG_IsOK(r) ? (r & ~SWIG_TMPOBJMASK) : r)

#define SWIG_IsTmpObj(r)           (SWIG_IsOK(r) && (r & SWIG_TMPOBJMASK))

/* Cast-Rank Mode */

#if defined(SWIG_CASTRANK_MODE)

#  ifndef SWIG_TypeRank

#    define SWIG_TypeRank             unsigned long

#  endif

#  ifndef SWIG_MAXCASTRANK            /* Default cast allowed */

#    define SWIG_MAXCASTRANK          (2)

#  endif

#  define SWIG_CASTRANKMASK          ((SWIG_CASTRANKLIMIT) -1)

#  define SWIG_CastRank(r)           (r & SWIG_CASTRANKMASK)

SWIGINTERNINLINE int SWIG_AddCast(int r) { 

  return SWIG_IsOK(r) ? ((SWIG_CastRank(r) < SWIG_MAXCASTRANK) ? (r + 1) : SWIG_ERROR) : r;

}

SWIGINTERNINLINE int SWIG_CheckState(int r) { 

  return SWIG_IsOK(r) ? SWIG_CastRank(r) + 1 : 0; 

}

#else /* no cast-rank mode */

#  define SWIG_AddCast

#  define SWIG_CheckState(r) (SWIG_IsOK(r) ? 1 : 0)

#endif

#include <string.h>

#ifdef __cplusplus

extern "C" {

#endif

typedef void *(*swig_converter_func)(void *, int *);

typedef struct swig_type_info *(*swig_dycast_func)(void **);

/* Structure to store information on one type */

typedef struct swig_type_info {

  const char             *name;			/* mangled name of this type */

  const char             *str;			/* human readable name of this type */

  swig_dycast_func        dcast;		/* dynamic cast function down a hierarchy */

  struct swig_cast_info  *cast;			/* linked list of types that can cast into this type */

  void                   *clientdata;		/* language specific type data */

  int                    owndata;		/* flag if the structure owns the clientdata */

} swig_type_info;

/* Structure to store a type and conversion function used for casting */

typedef struct swig_cast_info {

  swig_type_info         *type;			/* pointer to type that is equivalent to this type */

  swig_converter_func     converter;		/* function to cast the void pointers */

  struct swig_cast_info  *next;			/* pointer to next cast in linked list */

  struct swig_cast_info  *prev;			/* pointer to the previous cast */

} swig_cast_info;

/* Structure used to store module information

 * Each module generates one structure like this, and the runtime collects

 * all of these structures and stores them in a circularly linked list.*/

typedef struct swig_module_info {

  swig_type_info         **types;		/* Array of pointers to swig_type_info structures that are in this module */

  size_t                 size;		        /* Number of types in this module */

  struct swig_module_info *next;		/* Pointer to next element in circularly linked list */

  swig_type_info         **type_initial;	/* Array of initially generated type structures */

  swig_cast_info         **cast_initial;	/* Array of initially generated casting structures */

  void                    *clientdata;		/* Language specific module data */

} swig_module_info;

/* 

  Compare two type names skipping the space characters, therefore

  "char*" == "char *" and "Class<int>" == "Class<int >", etc.

  Return 0 when the two name types are equivalent, as in

  strncmp, but skipping ' '.

*/

SWIGRUNTIME int

SWIG_TypeNameComp(const char *f1, const char *l1,

		  const char *f2, const char *l2) {

  for (;(f1 != l1) && (f2 != l2); ++f1, ++f2) {

    while ((*f1 == ' ') && (f1 != l1)) ++f1;

    while ((*f2 == ' ') && (f2 != l2)) ++f2;

    if (*f1 != *f2) return (*f1 > *f2) ? 1 : -1;

  }

  return (int)((l1 - f1) - (l2 - f2));

}

/*

  Check type equivalence in a name list like <name1>|<name2>|...

  Return 0 if not equal, 1 if equal

*/

SWIGRUNTIME int

SWIG_TypeEquiv(const char *nb, const char *tb) {

  int equiv = 0;

  const char* te = tb + strlen(tb);

  const char* ne = nb;

  while (!equiv && *ne) {

    for (nb = ne; *ne; ++ne) {

      if (*ne == '|') break;

    }

    equiv = (SWIG_TypeNameComp(nb, ne, tb, te) == 0) ? 1 : 0;

    if (*ne) ++ne;

  }

  return equiv;

}

/*

  Check type equivalence in a name list like <name1>|<name2>|...

  Return 0 if equal, -1 if nb < tb, 1 if nb > tb

*/

SWIGRUNTIME int

SWIG_TypeCompare(const char *nb, const char *tb) {

  int equiv = 0;

  const char* te = tb + strlen(tb);

  const char* ne = nb;

  while (!equiv && *ne) {

    for (nb = ne; *ne; ++ne) {

      if (*ne == '|') break;

    }

    equiv = (SWIG_TypeNameComp(nb, ne, tb, te) == 0) ? 1 : 0;

    if (*ne) ++ne;

  }

  return equiv;

}

/*

  Check the typename

*/

SWIGRUNTIME swig_cast_info *

SWIG_TypeCheck(const char *c, swig_type_info *ty) {

  if (ty) {

    swig_cast_info *iter = ty->cast;

    while (iter) {

      if (strcmp(iter->type->name, c) == 0) {

        if (iter == ty->cast)

          return iter;

        /* Move iter to the top of the linked list */

        iter->prev->next = iter->next;

        if (iter->next)

          iter->next->prev = iter->prev;

        iter->next = ty->cast;

        iter->prev = 0;

        if (ty->cast) ty->cast->prev = iter;

        ty->cast = iter;

        return iter;

      }

      iter = iter->next;

    }

  }

  return 0;

}

/* 

  Identical to SWIG_TypeCheck, except strcmp is replaced with a pointer comparison

*/

SWIGRUNTIME swig_cast_info *

SWIG_TypeCheckStruct(swig_type_info *from, swig_type_info *ty) {

  if (ty) {

    swig_cast_info *iter = ty->cast;

    while (iter) {

      if (iter->type == from) {

        if (iter == ty->cast)

          return iter;

        /* Move iter to the top of the linked list */

        iter->prev->next = iter->next;

        if (iter->next)

          iter->next->prev = iter->prev;

        iter->next = ty->cast;

        iter->prev = 0;

        if (ty->cast) ty->cast->prev = iter;

        ty->cast = iter;

        return iter;

      }

      iter = iter->next;

    }

  }

  return 0;

}

/*

  Cast a pointer up an inheritance hierarchy

*/

SWIGRUNTIMEINLINE void *

SWIG_TypeCast(swig_cast_info *ty, void *ptr, int *newmemory) {

  return ((!ty) || (!ty->converter)) ? ptr : (*ty->converter)(ptr, newmemory);

}

/* 

   Dynamic pointer casting. Down an inheritance hierarchy

*/

SWIGRUNTIME swig_type_info *

SWIG_TypeDynamicCast(swig_type_info *ty, void **ptr) {

  swig_type_info *lastty = ty;

  if (!ty || !ty->dcast) return ty;

  while (ty && (ty->dcast)) {

    ty = (*ty->dcast)(ptr);

    if (ty) lastty = ty;

  }

  return lastty;

}

/*

  Return the name associated with this type

*/

SWIGRUNTIMEINLINE const char *

SWIG_TypeName(const swig_type_info *ty) {

  return ty->name;

}

/*

  Return the pretty name associated with this type,

  that is an unmangled type name in a form presentable to the user.

*/

SWIGRUNTIME const char *

SWIG_TypePrettyName(const swig_type_info *type) {

  /* The "str" field contains the equivalent pretty names of the

     type, separated by vertical-bar characters.  We choose

     to print the last name, as it is often (?) the most

     specific. */

  if (!type) return NULL;

  if (type->str != NULL) {

    const char *last_name = type->str;

    const char *s;

    for (s = type->str; *s; s++)

      if (*s == '|') last_name = s+1;

    return last_name;

  }

  else

    return type->name;

}

/* 

   Set the clientdata field for a type

*/

SWIGRUNTIME void

SWIG_TypeClientData(swig_type_info *ti, void *clientdata) {

  swig_cast_info *cast = ti->cast;

  /* if (ti->clientdata == clientdata) return; */

  ti->clientdata = clientdata;

  while (cast) {

    if (!cast->converter) {

      swig_type_info *tc = cast->type;

      if (!tc->clientdata) {

	SWIG_TypeClientData(tc, clientdata);

      }

    }    

    cast = cast->next;

  }

}

SWIGRUNTIME void

SWIG_TypeNewClientData(swig_type_info *ti, void *clientdata) {

  SWIG_TypeClientData(ti, clientdata);