Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

"""

Demonstrateur OpenCL d'interaction NCorps

Emmanuel QUEMENER <emmanuel.quemener@ens-lyon.fr> CeCILLv2

"""

import getopt

import sys

import time

import numpy as np

import pyopencl as cl

import pyopencl.array as cl_array

from numpy.random import randint as nprnd

def DictionariesAPI():

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

    Computing={'FP32':0,'FP64':1}

    return(Marsaglia,Computing)

BlobOpenCL= """

#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 TFP32 0

#define TFP64 1

#define LENGTH 1.e0f

#if TYPE == TFP32

#define MYFLOAT4 float4

#define MYFLOAT8 float8

#define MYFLOAT float

#define DISTANCE fast_distance

#else

#if defined(cl_khr_fp64)  // Khronos extension available?

#pragma OPENCL EXTENSION cl_khr_fp64 : enable

#define DOUBLE_SUPPORT_AVAILABLE

#elif defined(cl_amd_fp64)  // AMD extension available?

#pragma OPENCL EXTENSION cl_amd_fp64 : enable

#define DOUBLE_SUPPORT_AVAILABLE

#endif

#define MYFLOAT4 double4

#define MYFLOAT8 double8

#define MYFLOAT double

#define DISTANCE distance

#endif

#define MWCfp (MYFLOAT)(MWC * 2.3283064365386963e-10f)

#define KISSfp (MYFLOAT)(KISS * 2.3283064365386963e-10f)

#define SHR3fp (MYFLOAT)(SHR3 * 2.3283064365386963e-10f)

#define CONGfp (MYFLOAT)(CONG * 2.3283064365386963e-10f)

#define PI (MYFLOAT)3.141592653589793238462643197169399375105820974944592307816406286e0f

#define SMALL_NUM 1.e-9f

MYFLOAT4 Interaction(MYFLOAT4 m,MYFLOAT4 n)

{

    private MYFLOAT r=DISTANCE(n,m);

    return((n-m)/(MYFLOAT)(r*r*r));

}

MYFLOAT4 InteractionCore(MYFLOAT4 m,MYFLOAT4 n)

{

    private MYFLOAT core=(MYFLOAT)1.e5f;

    private MYFLOAT r=DISTANCE(n,m);

    private MYFLOAT d=r*r+core*core;

    return(core*(n-m)/(MYFLOAT)(d*d));

}

MYFLOAT PairPotential(MYFLOAT4 m,MYFLOAT4 n)

{

    return((MYFLOAT)(-1.e0f)/(DISTANCE(n,m)));

}

MYFLOAT AtomicPotential(__global MYFLOAT4* clDataX,int gid)

{

    private MYFLOAT potential=(MYFLOAT)0.e0f;

    private MYFLOAT4 x=clDataX[gid];

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        potential+=PairPotential(x,clDataX[i]);

    }

    barrier(CLK_GLOBAL_MEM_FENCE);

    return(potential);

}

MYFLOAT AtomicPotentialCoM(__global MYFLOAT4* clDataX,__global MYFLOAT4* clCoM,int gid)

{

    return(PairPotential(clDataX[gid],clCoM[0]));

}

MYFLOAT8 AtomicRungeKutta(__global MYFLOAT4* clDataInX,__global MYFLOAT4* clDataInV,int gid,MYFLOAT dt)

{

    private MYFLOAT4 a0,v0,x0,a1,v1,x1,a2,v2,x2,a3,v3,x3,a4,v4,x4,xf,vf;

    a0=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    v0=(MYFLOAT4)clDataInV[gid];

    x0=(MYFLOAT4)clDataInX[gid];

    int N = get_global_size(0);

    for (int i=0;i<N;i++)

    {

        if (gid != i)

        a0+=Interaction(x0,clDataInX[i]);

    }

    a1=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    v1=v0+a0*dt;

    x1=x0+v0*dt;

    for (int i=0;i<N;i++)

    {

        if (gid != i)

        a1+=Interaction(x1,clDataInX[i]);

    }

    a2=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    v2=v0+a1*dt*(MYFLOAT)5.e-1f;

    x2=x0+v1*dt*(MYFLOAT)5.e-1f;

    for (int i=0;i<N;i++)

    {

        if (gid != i)

        a2+=Interaction(x2,clDataInX[i]);

    }

    a3=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    v3=v0+a2*dt*(MYFLOAT)5.e-1f;

    x3=x0+v2*dt*(MYFLOAT)5.e-1f;

    for (int i=0;i<N;i++)

    {

        if (gid != i)

        a3+=Interaction(x3,clDataInX[i]);

    }

    a4=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    v4=v0+a3*dt;

    x4=x0+v3*dt;

    for (int i=0;i<N;i++)

    {

        if (gid != i)

        a4+=Interaction(x4,clDataInX[i]);

    }

    xf=x0+dt*(v1+(MYFLOAT)2.e0f*(v2+v3)+v4)/(MYFLOAT)6.e0f;

    vf=v0+dt*(a1+(MYFLOAT)2.e0f*(a2+a3)+a4)/(MYFLOAT)6.e0f;

    return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,1.e0f));

}

// Elements from : http://doswa.com/2009/01/02/fourth-order-runge-kutta-numerical-integration.html

MYFLOAT8 AtomicHeun(__global MYFLOAT4* clDataInX,__global MYFLOAT4* clDataInV,int gid,MYFLOAT dt)

{

    private MYFLOAT4 x,v,a,xi,vi,ai,xf,vf;

    x=(MYFLOAT4)clDataInX[gid];

    v=(MYFLOAT4)clDataInV[gid];

    a=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        a+=Interaction(x,clDataInX[i]);

    }

    vi=v+dt*a;

    xi=x+dt*vi;

    ai=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        ai+=Interaction(xi,clDataInX[i]);

    }

    vf=v+dt*(a+ai)/(MYFLOAT)2.e0f;

    xf=x+dt*(v+vi)/(MYFLOAT)2.e0f;

    return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,0.e0f));

}

MYFLOAT8 AtomicImplicitEuler(__global MYFLOAT4* clDataInX,__global MYFLOAT4* clDataInV,int gid,MYFLOAT dt)

{

    private MYFLOAT4 x,v,a,xf,vf;

    x=(MYFLOAT4)clDataInX[gid];

    v=(MYFLOAT4)clDataInV[gid];

    a=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        a+=Interaction(x,clDataInX[i]);

    }

    vf=v+dt*a;

    xf=x+dt*vf;

    return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,0.e0f));

}

MYFLOAT8 AtomicExplicitEuler(__global MYFLOAT4* clDataInX,__global MYFLOAT4* clDataInV,int gid,MYFLOAT dt)

{

    MYFLOAT4 x,v,a,xf,vf;

    x=(MYFLOAT4)clDataInX[gid];

    v=(MYFLOAT4)clDataInV[gid];

    a=(MYFLOAT4)(0.e0f,0.e0f,0.e0f,0.e0f);

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        a+=Interaction(x,clDataInX[i]);

    }

    vf=v+dt*a;

    xf=x+dt*v;

    return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,0.e0f));

}

__kernel void SplutterPoints(__global MYFLOAT4* clDataX, MYFLOAT box,

                             uint seed_z,uint seed_w)

{

    int gid = get_global_id(0);

    uint z=seed_z+(uint)gid;

    uint w=seed_w-(uint)gid;

    MYFLOAT x0=box*(MYFLOAT)(MWCfp-(MYFLOAT)5.e-1f);

    MYFLOAT y0=box*(MYFLOAT)(MWCfp-(MYFLOAT)5.e-1f);

    MYFLOAT z0=box*(MYFLOAT)(MWCfp-(MYFLOAT)5.e-1f);

    clDataX[gid].s0123 = (MYFLOAT4) (x0,y0,z0,1.e0f);

}

__kernel void SplutterStress(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,__global MYFLOAT4* clCoM, MYFLOAT velocity,uint seed_z,uint seed_w)

{

    int gid = get_global_id(0);

    MYFLOAT N = (MYFLOAT)get_global_size(0);

    uint z=seed_z+(uint)gid;

    uint w=seed_w-(uint)gid;

    if (velocity<SMALL_NUM) {

       MYFLOAT4 SpeedVector=(MYFLOAT4)normalize(cross(clDataX[gid],clCoM[0]))*sqrt(-AtomicPotential(clDataX,gid)/(MYFLOAT)2.e0f);

       clDataV[gid]=SpeedVector;

    }

    else

    {

       // cast to float for sin,cos are NEEDED by Mesa FP64 implementation!

       MYFLOAT theta=MWCfp*PI;

       MYFLOAT phi=MWCfp*PI*(MYFLOAT)2.e0f;

       MYFLOAT sinTheta=sin((float)theta);

       clDataV[gid].s0=velocity*sinTheta*cos((float)phi);

       clDataV[gid].s1=velocity*sinTheta*sin((float)phi);

       clDataV[gid].s2=velocity*cos((float)theta);

       clDataV[gid].s3=(MYFLOAT)1.e0f;

    }

}

__kernel void RungeKutta(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,MYFLOAT h)

{

    int gid = get_global_id(0);

    MYFLOAT8 clDataGid=AtomicRungeKutta(clDataX,clDataV,gid,h);

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.lo;

    clDataV[gid]=clDataGid.hi;

}

__kernel void ImplicitEuler(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,MYFLOAT h)

{

    int gid = get_global_id(0);

    MYFLOAT8 clDataGid=AtomicImplicitEuler(clDataX,clDataV,gid,h);

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.lo;

    clDataV[gid]=clDataGid.hi;

}

__kernel void Heun(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,MYFLOAT h)

{

    int gid = get_global_id(0);

    MYFLOAT8 clDataGid=AtomicHeun(clDataX,clDataV,gid,h);

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.lo;

    clDataV[gid]=clDataGid.hi;

}

__kernel void ExplicitEuler(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,MYFLOAT h)

{

    int gid = get_global_id(0);

    MYFLOAT8 clDataGid=AtomicExplicitEuler(clDataX,clDataV,gid,h);

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.lo;

    clDataV[gid]=clDataGid.hi;

}

__kernel void CoMPotential(__global MYFLOAT4* clDataX,__global MYFLOAT4* clCoM,__global MYFLOAT* clPotential)

{

    int gid = get_global_id(0);

    clPotential[gid]=PairPotential(clDataX[gid],clCoM[0]);

}

__kernel void Potential(__global MYFLOAT4* clDataX,__global MYFLOAT* clPotential)

{

    int gid = get_global_id(0);

    MYFLOAT potential=(MYFLOAT)0.e0f;

    MYFLOAT4 x=clDataX[gid];

    for (int i=0;i<get_global_size(0);i++)

    {

        if (gid != i)

        potential+=PairPotential(x,clDataX[i]);

    }

    barrier(CLK_GLOBAL_MEM_FENCE);

    clPotential[gid]=potential*(MYFLOAT)5.e-1f;

}

__kernel void CenterOfMass(__global MYFLOAT4* clDataX,__global MYFLOAT4* clCoM,int Size)

{

    MYFLOAT4 CoM=clDataX[0];

    for (int i=1;i<Size;i++)

    {

        CoM+=clDataX[i];

    }

    barrier(CLK_GLOBAL_MEM_FENCE);

    clCoM[0]=(MYFLOAT4)(CoM.s0,CoM.s1,CoM.s2,1.e0f)/(MYFLOAT)Size;

}

__kernel void Kinetic(__global MYFLOAT4* clDataV,__global MYFLOAT* clKinetic)

{

    int gid = get_global_id(0);

    barrier(CLK_GLOBAL_MEM_FENCE);

    MYFLOAT d=(MYFLOAT)length(clDataV[gid]);

    clKinetic[gid]=(MYFLOAT)5.e-1f*(MYFLOAT)(d*d);

}

"""

def Energy(MyData):

    return(sum(MyData*MyData))

if __name__=='__main__':

    # ValueType

    ValueType='FP32'

    class MyFloat(np.float32):pass

    #    clType8=cl_array.vec.float8

    clType4=cl_array.vec.float4

    # Set defaults values

    np.set_printoptions(precision=2)

    # Id of Device : 1 is for first find !

    Device=0

    # Iterations is integer

    Number=2

    # Size of box

    SizeOfBox=MyFloat(1.)

    # Initial velocity of particules

    Velocity=MyFloat(1.)

    # Redo the last process

    Iterations=int(np.pi*1024)

    # Step

    Step=MyFloat(1./1024)

    # Method of integration

    Method='ImplicitEuler'

    # InitialRandom

    InitialRandom=False

    # RNG Marsaglia Method

    RNG='MWC'

    # CheckEnergies

    CheckEnergies=False

    # Display samples in 3D

    GraphSamples=False

    # Viriel Distribution of stress

    VirielStress=True

    HowToUse='%s -h [Help] -r [InitialRandom] -e [VirielStress] -g [GraphSamples] -c [CheckEnergies] -d <DeviceId> -n <NumberOfParticules> -z <SizeOfBox> -v <Velocity> -s <Step> -i <Iterations> -m <ImplicitEuler|RungeKutta|ExplicitEuler|Heun> -t <FP32|FP64>'

    try:

        opts, args = getopt.getopt(sys.argv[1:],"rehgcd:n:z:v:i:s:m:t:",["random","viriel","graph","check","device=","number=","size=","velocity=","iterations=","step=","method=","valuetype="])

    except getopt.GetoptError:

        print(HowToUse % sys.argv[0])

        sys.exit(2)

    for opt, arg in opts:

        if opt == '-h':

            print(HowToUse % sys.argv[0])

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

            try:

                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(),deviceType,device.name.lstrip()))

                        Id=Id+1

                sys.exit()

            except ImportError:

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

                sys.exit()

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

            if arg=='FP64':

                class MyFloat(np.float64): pass

                clType4=cl_array.vec.double4

            else:

                class MyFloat(np.float32):pass

                clType4=cl_array.vec.float4

            ValueType = arg

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

            Device=int(arg)

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

            Method=arg

        elif opt in ("-n", "--number"):

            Number=int(arg)

        elif opt in ("-z", "--size"):

            SizeOfBox=MyFloat(arg)

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

            Velocity=MyFloat(arg)

            VirielStress=False

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

            Step=MyFloat(arg)

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

            Iterations=int(arg)

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

            InitialRandom=True

        elif opt in ("-c", "--check"):

            CheckEnergies=True

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

            GraphSamples=True

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

            VirielStress=True

    SizeOfBox=MyFloat(Number*SizeOfBox)

    Velocity=MyFloat(Velocity)

    Step=MyFloat(Step)

    print("Device choosed : %s" % Device)

    print("Number of particules : %s" % Number)

    print("Size of Box : %s" % SizeOfBox)

    print("Initial velocity : %s" % Velocity)

    print("Number of iterations : %s" % Iterations)

    print("Step of iteration : %s" % Step)

    print("Method of resolution : %s" % Method)

    print("Initial Random for RNG Seed : %s" % InitialRandom)

    print("Check for Energies : %s" % CheckEnergies)

    print("Graph for Samples : %s" % GraphSamples)

    print("ValueType is : %s" % ValueType)

    print("Viriel distribution of stress %s" % VirielStress)

    # Create Numpy array of CL vector with 8 FP32

    MyCoM = np.zeros(1,dtype=clType4)

    MyDataX = np.zeros(Number, dtype=clType4)

    MyDataV = np.zeros(Number, dtype=clType4)

    MyPotential = np.zeros(Number, dtype=MyFloat)

    MyKinetic = np.zeros(Number, dtype=MyFloat)

    Marsaglia,Computing=DictionariesAPI()

    # Scan the OpenCL arrays

    Id=0

    HasXPU=False

    for platform in cl.get_platforms():

        for device in platform.get_devices():

            if Id==Device:

                PlatForm=platform

                XPU=device

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

                print("Platform selected: ",platform.name)

                HasXPU=True

            Id+=1

    if HasXPU==False:

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

        sys.exit()

    # Create Context

    try:

        ctx = cl.Context([XPU])

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

    except:

        print("Crash during context creation")

    print(Marsaglia[RNG],Computing[ValueType])

    # Build all routines used for the computing

    #BuildOptions="-DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType])

    #BuildOptions="-cl-mad-enable -cl-fast-relaxed-math -DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType])

    BuildOptions="-cl-mad-enable -cl-kernel-arg-info -cl-fast-relaxed-math -cl-std=CL1.2 -DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType])

    if 'Intel' in PlatForm.name or 'Clover' in PlatForm.name or 'Portable' in PlatForm.name :

        MyRoutines = cl.Program(ctx, BlobOpenCL).build(options = BuildOptions)

    else:

        MyRoutines = cl.Program(ctx, BlobOpenCL).build(options = BuildOptions+" -cl-strict-aliasing")

    mf = cl.mem_flags

    # clDataX = cl.Buffer(ctx, mf.READ_WRITE, MyDataX.nbytes)

    # clDataV = cl.Buffer(ctx, mf.READ_WRITE, MyDataV.nbytes)

    # clPotential = cl.Buffer(ctx, mf.READ_WRITE, MyPotential.nbytes)

    # clKinetic = cl.Buffer(ctx, mf.READ_WRITE, MyKinetic.nbytes)

    # clCoM = cl.Buffer(ctx, mf.READ_WRITE, MyCoM.nbytes)

    clDataX = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyDataX)

    clDataV = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyDataV)

    clPotential = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyPotential)

    clKinetic = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyKinetic)

    clCoM = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyCoM)

    print('All particles superimposed.')

    print(SizeOfBox.dtype)

    # Set particles to RNG points

    if InitialRandom:

        MyRoutines.SplutterPoints(queue,(Number,1),None,clDataX,SizeOfBox,np.uint32(nprnd(2**32)),np.uint32(nprnd(2**32)))

    else:

        MyRoutines.SplutterPoints(queue,(Number,1),None,clDataX,SizeOfBox,np.uint32(110271),np.uint32(250173))

    print('All particules distributed')

    CLLaunch=MyRoutines.CenterOfMass(queue,(1,1),None,clDataX,clCoM,np.int32(Number))

    CLLaunch.wait()

    cl.enqueue_copy(queue,MyCoM,clCoM)

    print('Center Of Mass: (%s,%s,%s)' % (MyCoM[0][0],MyCoM[0][1],MyCoM[0][2]))

    if VirielStress:

        CLLaunch=MyRoutines.SplutterStress(queue,(Number,1),None,clDataX,clDataV,clCoM,MyFloat(0.),np.uint32(110271),np.uint32(250173))

    else:

        CLLaunch=MyRoutines.SplutterStress(queue,(Number,1),None,clDataX,clDataV,clCoM,Velocity,np.uint32(110271),np.uint32(250173))

    CLLaunch.wait()

    if GraphSamples:

        cl.enqueue_copy(queue, MyDataX, clDataX)

        t0=np.array([[MyDataX[0][0],MyDataX[0][1],MyDataX[0][2]]])

        t1=np.array([[MyDataX[1][0],MyDataX[1][1],MyDataX[1][2]]])

        tL=np.array([[MyDataX[-1][0],MyDataX[-1][1],MyDataX[-1][2]]])

    CLLaunch=MyRoutines.Potential(queue,(Number,1),None,clDataX,clPotential)

    CLLaunch=MyRoutines.Kinetic(queue,(Number,1),None,clDataV,clKinetic)

    CLLaunch.wait()

    cl.enqueue_copy(queue,MyPotential,clPotential)

    cl.enqueue_copy(queue,MyKinetic,clKinetic)

    print('Viriel=%s Potential=%s Kinetic=%s'% (np.sum(MyPotential)+2*np.sum(MyKinetic),np.sum(MyPotential),np.sum(MyKinetic)))

    if GraphSamples:

        cl.enqueue_copy(queue, MyDataX, clDataX)

        t0=np.array([[MyDataX[0][0],MyDataX[0][1],MyDataX[0][2]]])

        t1=np.array([[MyDataX[1][0],MyDataX[1][1],MyDataX[1][2]]])

        tL=np.array([[MyDataX[-1][0],MyDataX[-1][1],MyDataX[-1][2]]])

    time_start=time.time()

    for i in range(Iterations):

        if Method=="RungeKutta":

            CLLaunch=MyRoutines.RungeKutta(queue,(Number,1),None,clDataX,clDataV,Step)

        elif Method=="ExplicitEuler":

            CLLaunch=MyRoutines.ExplicitEuler(queue,(Number,1),None,clDataX,clDataV,Step)

        elif Method=="Heun":

            CLLaunch=MyRoutines.Heun(queue,(Number,1),None,clDataX,clDataV,Step)

        else:

            CLLaunch=MyRoutines.ImplicitEuler(queue,(Number,1),None,clDataX,clDataV,Step)

        CLLaunch.wait()

        if CheckEnergies:

            CLLaunch=MyRoutines.Potential(queue,(Number,1),None,clDataX,clPotential)

            CLLaunch=MyRoutines.Kinetic(queue,(Number,1),None,clDataV,clKinetic)

            CLLaunch.wait()

            cl.enqueue_copy(queue,MyPotential,clPotential)

            cl.enqueue_copy(queue,MyKinetic,clKinetic)

            print(np.sum(MyPotential)+2.*np.sum(MyKinetic),np.sum(MyPotential),np.sum(MyKinetic))

            print(MyPotential,MyKinetic)

        if GraphSamples:

            cl.enqueue_copy(queue, MyDataX, clDataX)

            t0=np.append(t0,[MyDataX[0][0],MyDataX[0][1],MyDataX[0][2]])

            t1=np.append(t1,[MyDataX[1][0],MyDataX[1][1],MyDataX[1][2]])

            tL=np.append(tL,[MyDataX[-1][0],MyDataX[-1][1],MyDataX[-1][2]])

    print("\nDuration on %s for each %s\n" % (Device,(time.time()-time_start)/Iterations))

    MyRoutines.CenterOfMass(queue,(1,1),None,clDataX,clCoM,np.int32(Number))

    CLLaunch=MyRoutines.Potential(queue,(Number,1),None,clDataX,clPotential)

    CLLaunch=MyRoutines.Kinetic(queue,(Number,1),None,clDataV,clKinetic)

    CLLaunch.wait()

    cl.enqueue_copy(queue,MyCoM,clCoM)

    cl.enqueue_copy(queue,MyPotential,clPotential)

    cl.enqueue_copy(queue,MyKinetic,clKinetic)

    print('Center Of Mass: (%s,%s,%s)' % (MyCoM[0][0],MyCoM[0][1],MyCoM[0][2]))

    print('Viriel=%s Potential=%s Kinetic=%s'% (np.sum(MyPotential)+2.*np.sum(MyKinetic),np.sum(MyPotential),np.sum(MyKinetic)))

    if GraphSamples:

        t0=np.transpose(np.reshape(t0,(Iterations+1,3)))

        t1=np.transpose(np.reshape(t1,(Iterations+1,3)))

        tL=np.transpose(np.reshape(tL,(Iterations+1,3)))

        import matplotlib.pyplot as plt

        from mpl_toolkits.mplot3d import Axes3D

        fig = plt.figure()

        ax = fig.gca(projection='3d')

        ax.scatter(t0[0],t0[1],t0[2], marker='^',color='blue')

        ax.scatter(t1[0],t1[1],t1[2], marker='o',color='red')

        ax.scatter(tL[0],tL[1],tL[2], marker='D',color='green')

        ax.set_xlabel('X Label')

        ax.set_ylabel('Y Label')

        ax.set_zlabel('Z Label')

        plt.show()

    clDataX.release()

    clDataV.release()

    clKinetic.release()

    clPotential.release()