Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python3

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

"""

NBody Demonstrator implemented in OpenCL, rendering OpenGL

By default, rendering in OpenGL is disabled. Add -g option to activate.

Part of matrix programs from: https://forge.cbp.ens-lyon.fr/svn/bench4gpu/

CC BY-NC-SA 2011 : Emmanuel QUEMENER <emmanuel.quemener@ens-lyon.fr>

Thanks to Andreas Klockner for PyOpenCL:

http://mathema.tician.de/software/pyopencl

"""

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

import string, sys

from OpenGL.GL import *

from OpenGL.GLUT import *

def DictionariesAPI():

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

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

    Interaction={'Force':0,'Potential':1}

    Artevasion={'None':0,'NegExp':1,'CorRad':2}

    return(Marsaglia,Computing,Interaction,Artevasion)

BlobOpenCL= """

#define TFP32 0

#define TFP64 1

#define TFORCE 0

#define TPOTENTIAL 1

#define NONE 0

#define NEGEXP 1

#define CORRAD 2

#if TYPE == TFP32

#define MYFLOAT4 float4

#define MYFLOAT8 float8

#define MYFLOAT float

#define DISTANCE fast_distance

#else

#define MYFLOAT4 double4

#define MYFLOAT8 double8

#define MYFLOAT double

#define DISTANCE distance

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

#pragma OPENCL EXTENSION cl_khr_fp64 : enable

#endif

#endif

#define znew  ((zmwc=36969*(zmwc&65535)+(zmwc>>16))<<16)

#define wnew  ((wmwc=18000*(wmwc&65535)+(wmwc>>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 (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.141592653589793238e0f

#define SMALL_NUM (MYFLOAT)1.e-9f

#define CoreRadius (MYFLOAT)(1.e0f)

// Create my own Distance implementation: distance buggy on Oland AMD chipset

MYFLOAT MyDistance(MYFLOAT4 n,MYFLOAT4 m)

{

    private MYFLOAT x2,y2,z2;

    x2=n.s0-m.s0;

    x2*=x2;

    y2=n.s1-m.s1;

    y2*=y2;

    z2=n.s2-m.s2;

    z2*=z2;

    return(sqrt(x2+y2+z2));

}

// Potential between 2 m,n bodies

MYFLOAT PairPotential(MYFLOAT4 m,MYFLOAT4 n)

#if ARTEVASION == NEGEXP

// Add exp(-r) to numerator to avoid divergence for low distances

{

    MYFLOAT r=DISTANCE(n,m);

    return((-1.e0f+exp(-r))/r);

}

#elif ARTEVASION == CORRAD

// Add Core Radius to avoid divergence for low distances

{

    MYFLOAT r=DISTANCE(n,m);

    return(-1.e0f/sqrt(r*r+CoreRadius*CoreRadius));

}

#else

// Classical potential in 1/r

{

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

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

}

#endif

// Interaction based of Force as gradient of Potential

MYFLOAT4 Interaction(MYFLOAT4 m,MYFLOAT4 n)

#if INTERACTION == TFORCE

#if ARTEVASION == NEGEXP

// Force gradient of potential, set as (1-exp(-r))/r

{

    private MYFLOAT r=MyDistance(n,m);

    private MYFLOAT num=1.e0f+exp(-r)*(r-1.e0f);

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

}

#elif ARTEVASION == CORRAD

// Force gradient of potential, (Core Radius) set as 1/sqrt(r**2+CoreRadius**2)

{

    private MYFLOAT r=MyDistance(n,m);

    private MYFLOAT den=sqrt(r*r+CoreRadius*CoreRadius);

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

}

#else

// Simplest implementation of force (equals to acceleration)

// seems to bo bad (numerous artevasions)

// MYFLOAT4 InteractionForce(MYFLOAT4 m,MYFLOAT4 n)

{

    private MYFLOAT r=MyDistance(n,m);

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

}

#endif

#else

// Force definited as gradient of potential

// Estimate potential and proximate potential to estimate force

{

    // 1/1024 seems to be a good factor: larger one provides bad results

    private MYFLOAT epsilon=(MYFLOAT)(1.e0f/1024);

    private MYFLOAT4 er=normalize(n-m);

    private MYFLOAT4 dr=er*(MYFLOAT)epsilon;

    return(er/epsilon*(PairPotential(m,n)-PairPotential(m+dr,n)));

}

#endif

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]));

}

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

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;

    MYFLOAT4 DT=dt*(MYFLOAT4)(1.e0f,1.e0f,1.e0f,1.e0f);

    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 (private 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=a0*dt+v0;

    x1=v0*dt+x0;

    for (private int j=0;j<N;j++)

    {

        if (gid != j)

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

    }

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

    v2=a1*(MYFLOAT)(dt/2.e0f)+v0;

    x2=v1*(MYFLOAT)(dt/2.e0f)+x0;

    for (private int k=0;k<N;k++)

    {

        if (gid != k)

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

    }

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

    v3=a2*(MYFLOAT)(dt/2.e0f)+v0;

    x3=v2*(MYFLOAT)(dt/2.e0f)+x0;

    for (private int l=0;l<N;l++)

    {

        if (gid != l)

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

    }

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

    v4=a3*dt+v0;

    x4=v3*dt+x0;

    for (private int m=0;m<N;m++)

    {

        if (gid != m)

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

    }

    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,0.e0f,vf.s0,vf.s1,vf.s2,0.e0f));

}

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

{

    private MYFLOAT4 x0,v0,a0,x1,v1,a1,xf,vf;

    MYFLOAT4 Dt=dt*(MYFLOAT4)(1.e0f,1.e0f,1.e0f,1.e0f);

    x0=(MYFLOAT4)clDataInX[gid];

    v0=(MYFLOAT4)clDataInV[gid];

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

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

    {

        if (gid != i)

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

    }

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

    //v1=v0+dt*a0;

    //x1=x0+dt*v0;

    v1=dt*a0+v0;

    x1=dt*v0+x0;

    for (private int j=0;j<get_global_size(0);j++)

    {

        if (gid != j)

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

    }

    vf=v0+dt*(a0+a1)/(MYFLOAT)2.e0f;

    xf=x0+dt*(v0+v1)/(MYFLOAT)2.e0f;

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

}

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

{

    MYFLOAT4 x0,v0,a,xf,vf;

    x0=(MYFLOAT4)clDataInX[gid];

    v0=(MYFLOAT4)clDataInV[gid];

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

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

    {

        if (gid != i)

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

    }

    vf=v0+dt*a;

    xf=x0+dt*vf;

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

}

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

{

    MYFLOAT4 x0,v0,a,xf,vf;

    x0=(MYFLOAT4)clDataInX[gid];

    v0=(MYFLOAT4)clDataInV[gid];

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

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

    {

        if (gid != i)

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

    }

    vf=v0+dt*a;

    xf=x0+dt*v0;

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

}

__kernel void InBallSplutterPoints(__global MYFLOAT4* clDataX,

                                   MYFLOAT diameter,uint seed_z,uint seed_w)

{

    private int gid=get_global_id(0);

    private uint zmwc=seed_z+gid;

    private uint wmwc=seed_w+(gid+1)%2;

    private MYFLOAT Heat;

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

    {

        Heat=MWCfp;

    }

// More accurate distribution based on spherical coordonates

// Disactivated because of AMD Oland GPU crash on launch

//     private MYFLOAT Radius,Theta,Phi,PosX,PosY,PosZ,SinTheta;

//     Radius=MWCfp*diameter/2.e0f;

//     Theta=(MYFLOAT)acos((float)(-2.e0f*MWCfp+1.0e0f));

//     Phi=(MYFLOAT)(2.e0f*PI*MWCfp);

//     SinTheta=sin((float)Theta);

//     PosX=cos((float)Phi)*Radius*SinTheta;

//     PosY=sin((float)Phi)*Radius*SinTheta;

//     PosZ=cos((float)Theta)*Radius;

//     clDataX[gid]=(MYFLOAT4)(PosX,PosY,PosZ,0.e0f);

    private MYFLOAT Radius=diameter/2.e0f;

    private MYFLOAT Length=diameter;

    private MYFLOAT4 Position;

    while (Length>Radius) {

       Position=(MYFLOAT4)((MWCfp-0.5e0f)*diameter,(MWCfp-0.5e0f)*diameter,(MWCfp-0.5e0f)*diameter,0.e0f);

       Length=(MYFLOAT)length((MYFLOAT4)Position);

    }

    clDataX[gid]=Position;

    barrier(CLK_GLOBAL_MEM_FENCE);

}

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

                             uint seed_z,uint seed_w)

{

    int gid=get_global_id(0);

    uint zmwc=seed_z+gid;

    uint wmwc=seed_w-gid;

    private MYFLOAT Heat;

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

    {

        Heat=MWCfp;

    }

    clDataX[gid]=(MYFLOAT4)((MWCfp-0.5e0f)*box,(MWCfp-0.5e0f)*box,(MWCfp-0.5e0f)*box,0.e0f);

    barrier(CLK_GLOBAL_MEM_FENCE);

}

__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 zmwc=seed_z+(uint)gid;

    uint wmwc=seed_w-(uint)gid;

    MYFLOAT4 CrossVector,SpeedVector,FromCoM;

    MYFLOAT Heat,ThetaA,PhiA,ThetaB,PhiB,Length,tA,tB,Polar;

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

    {

        Heat=MWCfp;

    }

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

    // Implemention on AMD Oland are probably broken in float

    FromCoM=(MYFLOAT4)(clDataX[gid]-clCoM[0]);

    Length=length(FromCoM);

    //Theta=acos(FromCoM.z/Length);

    //Phi=atan(FromCoM.y/FromCoM.x);

    // First tangential vector to sphere of length radius

    ThetaA=acos(FromCoM.x/Length)+5.e-1f*PI;

    PhiA=atan(FromCoM.y/FromCoM.z);

    // Second tangential vector to sphere of length radius

    ThetaB=acos((float)(FromCoM.x/Length));

    PhiB=atan((float)(FromCoM.y/FromCoM.z))+5.e-1f*PI;

    // (x,y) random coordonates to plane tangential to sphere

    Polar=MWCfp*2.e0f*PI;

    tA=cos((float)Polar);

    tB=sin((float)Polar);

    // Exception for 2 particules to ovoid shifting

    if (get_global_size(0)==2) {

       CrossVector=(MYFLOAT4)(1.e0f,1.e0f,1.e0f,0.e0f);

    } else {

       CrossVector.s0=tA*cos((float)ThetaA)+tB*cos((float)ThetaB);

       CrossVector.s1=tA*sin((float)ThetaA)*sin((float)PhiA)+tB*sin((float)ThetaB)*sin((float)PhiB);

       CrossVector.s2=tA*sin((float)ThetaA)*cos((float)PhiA)+tB*sin((float)ThetaB)*cos((float)PhiB);

       CrossVector.s3=0.e0f;

    }

    if (velocity<SMALL_NUM) {

       SpeedVector=(MYFLOAT4)normalize(cross(FromCoM,CrossVector))*sqrt((-AtomicPotential(clDataX,gid)/(MYFLOAT)2.e0f));

    }

    else

    {

       SpeedVector=(MYFLOAT4)((MWCfp-5e-1f)*velocity,(MWCfp-5e-1f)*velocity,

                              (MWCfp-5e-1f)*velocity,0.e0f);

    }

    clDataV[gid]=SpeedVector;

    barrier(CLK_GLOBAL_MEM_FENCE);

}

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

{

    private int gid = get_global_id(0);

    private MYFLOAT8 clDataGid;

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

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.s0123;

    clDataV[gid]=clDataGid.s4567;

}

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

{

    private int gid = get_global_id(0);

    private MYFLOAT8 clDataGid;

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

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.s0123;

    clDataV[gid]=clDataGid.s4567;

}

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

{

    private int gid = get_global_id(0);

    private MYFLOAT8 clDataGid;

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

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.s0123;

    clDataV[gid]=clDataGid.s4567;

}

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

{

    private int gid = get_global_id(0);

    private MYFLOAT8 clDataGid;

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

    barrier(CLK_GLOBAL_MEM_FENCE);

    clDataX[gid]=clDataGid.s0123;

    clDataV[gid]=clDataGid.s4567;

}

__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,0.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 MainOpenCL(clDataX,clDataV,Step,Method):

    time_start=time.time()

    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()

    Elapsed=time.time()-time_start

    return(Elapsed)

def display(*args):

    global MyDataX,MyDataV,clDataX,clDataV,Step,Method,Number,Iterations,Durations,Verbose,SpeedRendering

    glClearColor(0.0, 0.0, 0.0, 0.0)

    glClear(GL_COLOR_BUFFER_BIT)

    glColor3f(1.0,1.0,1.0)

    Elapsed=MainOpenCL(clDataX,clDataV,Step,Method)

    if SpeedRendering:

        cl.enqueue_copy(queue, MyDataV, clDataV)

        MyDataV.reshape(Number,4)[:,3]=1

        glVertexPointerf(MyDataV.reshape(Number,4))

    else:

        cl.enqueue_copy(queue, MyDataX, clDataX)

        MyDataX.reshape(Number,4)[:,3]=1

        glVertexPointerf(MyDataX.reshape(Number,4))

    if Verbose:

        print("Positions for #%s iteration: %s" % (Iterations,MyDataX))

    else:

        sys.stdout.write('.')

        sys.stdout.flush()

    Durations=np.append(Durations,MainOpenCL(clDataX,clDataV,Step,Method))

    glEnableClientState(GL_VERTEX_ARRAY)

    glDrawArrays(GL_POINTS, 0, Number)

    glDisableClientState(GL_VERTEX_ARRAY)

    glFlush()

    Iterations+=1

    glutSwapBuffers()

def halt():

    pass

def keyboard(k,x,y):

    global ViewRZ,SpeedRendering

    LC_Z = as_8_bit( 'z' )

    UC_Z = as_8_bit( 'Z' )

    Plus = as_8_bit( '+' )

    Minus = as_8_bit( '-' )

    Switch = as_8_bit( 's' )

    Zoom=1

    if k == LC_Z:

        ViewRZ += 1.0

    elif k == UC_Z:

        ViewRZ -= 1.0

    elif k == Plus:

        Zoom *= 2.0

    elif k == Minus:

        Zoom /= 2.0

    elif k == Switch:

        if SpeedRendering:

            SpeedRendering=False

        else:

            SpeedRendering=True

    elif ord(k) == 27: # Escape

        glutLeaveMainLoop()

        return(False)

    else:

        return

    glRotatef(ViewRZ, 0.0, 0.0, 1.0)

    glScalef(Zoom,Zoom,Zoom)

    glutPostRedisplay()

def special(k,x,y):

    global ViewRX, ViewRY

    Step=1.

    if k == GLUT_KEY_UP:

        ViewRX += Step

    elif k == GLUT_KEY_DOWN:

        ViewRX -= Step

    elif k == GLUT_KEY_LEFT:

        ViewRY += Step

    elif k == GLUT_KEY_RIGHT:

        ViewRY -= Step

    else:

        return

    glRotatef(ViewRX, 1.0, 0.0, 0.0)

    glRotatef(ViewRY, 0.0, 1.0, 0.0)

    glutPostRedisplay()

def setup_viewport():

    global SizeOfBox

    glMatrixMode(GL_PROJECTION)

    glLoadIdentity()

    glOrtho(-SizeOfBox, SizeOfBox, -SizeOfBox, SizeOfBox, -SizeOfBox, SizeOfBox)

    glutPostRedisplay()

def reshape(w, h):

    glViewport(0, 0, w, h)

    setup_viewport()

if __name__=='__main__':

    global Number,Step,clDataX,clDataV,MyDataX,MyDataV,Method,SizeOfBox,Iterations,Verbose,Durations

    # ValueType

    ValueType='FP32'

    class MyFloat(np.float32):pass

    #    clType8=cl_array.vec.float8

    # Set defaults values

    np.set_printoptions(precision=2)

    # Id of Device : 1 is for first find !

    Device=0

    # Number of bodies is integer

    Number=2

    # Number of iterations (for standalone execution)

    Iterations=10

    # Size of shape

    SizeOfShape=MyFloat(1.)

    # Initial velocity of particules

    Velocity=MyFloat(1.)

    # Step

    Step=MyFloat(1./32)

    # Method of integration

    Method='ImplicitEuler'

    # InitialRandom

    InitialRandom=False

    # RNG Marsaglia Method

    RNG='MWC'

    # Viriel Distribution of stress

    VirielStress=True

    # Verbose

    Verbose=False

    # OpenGL real time rendering

    OpenGL=False

    # Speed rendering

    SpeedRendering=False

    # Counter ArtEvasions Measures (artefact evasion)

    CoArEv='None'

    # Shape to distribute

    Shape='Ball'

    # Type of Interaction

    InterType='Force'

    HowToUse='%s -h [Help] -r [InitialRandom] -g [OpenGL] -e [VirielStress] -o [Verbose] -p [Potential] -x <None|NegExp|CorRad> -d <DeviceId> -n <NumberOfParticules> -i <Iterations> -z <SizeOfBoxOrBall> -v <Velocity> -s <Step> -b <Ball|Box> -m <ImplicitEuler|RungeKutta|ExplicitEuler|Heun> -t <FP32|FP64>'

    try:

        opts, args = getopt.getopt(sys.argv[1:],"rpgehod:n:i:z:v:s:m:t:b:x:",["random","potential","coarev","opengl","viriel","verbose","device=","number=","iterations=","size=","velocity=","step=","method=","valuetype=","shape="])

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

                        # Failed now because of POCL implementation

                        #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

            else:

                class MyFloat(np.float32):pass

            ValueType = arg

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

            Device=int(arg)

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

            Method=arg

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

            Shape=arg

            if Shape!='Ball' or Shape!='Box':

                print('Wrong argument: set to Ball')

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