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

            Number=int(arg)

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

            Iterations=int(arg)

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

            SizeOfShape=MyFloat(arg)

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

            Velocity=MyFloat(arg)

            VirielStress=False

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

            Step=MyFloat(arg)

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

            InitialRandom=True

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

            CheckEnergies=True

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

            VirielStress=True

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

            OpenGL=True

        elif opt in ("-p", "--potential"):

            InterType='Potential'

        elif opt in ("-x", "--coarev"):

            CoArEv=arg

        elif opt in ("-o", "--verbose"):

            Verbose=True

    SizeOfShape=np.sqrt(MyFloat(SizeOfShape*Number))

    Velocity=MyFloat(Velocity)

    Step=MyFloat(Step)

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

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

    print("Size of Shape : %s" % SizeOfShape)

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

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

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

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

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

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

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

    print("OpenGL real time rendering : %s" % OpenGL)

    print("Speed rendering : %s" % SpeedRendering)

    print("Interaction type : %s" % InterType)

    print("Counter Artevasion type : %s" % CoArEv)

    # Create Numpy array of CL vector with 8 FP32    

    MyCoM = np.zeros(4,dtype=MyFloat)

    MyDataX = np.zeros(Number*4, dtype=MyFloat)

    MyDataV = np.zeros(Number*4, dtype=MyFloat)

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

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

    Marsaglia,Computing,Interaction,Artevasion=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")

    # Build all routines used for the computing

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

    BuildOptions="-cl-mad-enable -cl-fast-relaxed-math -DTRNG=%i -DTYPE=%i -DINTERACTION=%i -DARTEVASION=%i" % (Marsaglia[RNG],Computing[ValueType],Interaction[InterType],Artevasion[CoArEv])

    if 'Intel' in PlatForm.name or 'Experimental' 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

    # Read/Write approach for buffering

    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)

    # Write/HostPointer approach for buffering

    # 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.')

    # Set particles to RNG points

    if InitialRandom:

        seed_w=np.uint32(nprnd(2**32))

        seed_z=np.uint32(nprnd(2**32))

    else:

        seed_w=np.uint32(19710211)

        seed_z=np.uint32(20081010)

    if Shape=='Ball':

        MyRoutines.InBallSplutterPoints(queue,(Number,1),None,clDataX,SizeOfShape,seed_w,seed_z)

    else:

        MyRoutines.InBoxSplutterPoints(queue,(Number,1),None,clDataX,SizeOfShape,seed_w,seed_z)

    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 estimated: (%s,%s,%s)' % (MyCoM[0],MyCoM[1],MyCoM[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()

    print('All particules stressed')

    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('Energy estimated: Viriel=%s Potential=%s Kinetic=%s\n'% (np.sum(MyPotential)+2*np.sum(MyKinetic),np.sum(MyPotential),np.sum(MyKinetic)))

    if SpeedRendering:

        SizeOfBox=max(2*MyKinetic)

    else:

        SizeOfBox=SizeOfShape        

    if OpenGL:

        print('\tTiny documentation to interact OpenGL rendering:\n')

        print('\t<Left|Right> Rotate around X axis')

        print('\t  <Up|Down>  Rotate around Y axis')

        print('\t   <z|Z>     Rotate around Z axis')

        print('\t   <-|+>     Unzoom/Zoom')

        print('\t    <s>      Toggle to display Positions or Velocities')

        print('\t   <Esc>     Quit\n')

    wall_time_start=time.time()

    Durations=np.array([],dtype=MyFloat)

    print('Starting!')

    if OpenGL:

        global ViewRX,ViewRY,ViewRZ

        Iterations=0

        ViewRX,ViewRY,ViewRZ = 0.,0.,0.

        # Launch OpenGL Loop

        glutInit(sys.argv)

        glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB)

        glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE,GLUT_ACTION_CONTINUE_EXECUTION)

        glutInitWindowSize(512,512)

        glutCreateWindow(b'NBodyGL')

        setup_viewport()

        glutReshapeFunc(reshape)

        glutDisplayFunc(display)

        glutIdleFunc(display)

        #   glutMouseFunc(mouse)

        glutSpecialFunc(special)

        glutKeyboardFunc(keyboard)

        glutMainLoop()

    else:

        for iteration in range(Iterations):

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

            if Verbose:

                # print("Duration of #%s iteration: %s" % (iteration,Elapsed))

                cl.enqueue_copy(queue, MyDataX, clDataX)

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

            else:

                sys.stdout.write('.')

                sys.stdout.flush()

            Durations=np.append(Durations,Elapsed)

    print('\nEnding!')

    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('\nCenter Of Mass estimated: (%s,%s,%s)' % (MyCoM[0],MyCoM[1],MyCoM[2]))

    print('Energy estimated: Viriel=%s Potential=%s Kinetic=%s\n'% (np.sum(MyPotential)+2.*np.sum(MyKinetic),np.sum(MyPotential),np.sum(MyKinetic)))

    print("Duration stats on device %s with %s iterations :\n\tMean:\t%s\n\tMedian:\t%s\n\tStddev:\t%s\n\tMin:\t%s\n\tMax:\t%s\n\n\tVariability:\t%s\n" % (Device,Iterations,np.mean(Durations),np.median(Durations),np.std(Durations),np.min(Durations),np.max(Durations),np.std(Durations)/np.median(Durations)))

    # FPS: 1/Elapsed

    FPS=np.ones(len(Durations))

    FPS/=Durations

    print("FPS stats on device %s with %s iterations :\n\tMean:\t%s\n\tMedian:\t%s\n\tStddev:\t%s\n\tMin:\t%s\n\tMax:\t%s\n" % (Device,Iterations,np.mean(FPS),np.median(FPS),np.std(FPS),np.min(FPS),np.max(FPS)))

    # Contraction of Square*Size*Hertz: Size*Size/Elapsed

    Squertz=np.ones(len(Durations))

    Squertz*=Number*Number

    Squertz/=Durations

    print("Squertz in log10 & complete stats on device %s with %s iterations :\n\tMean:\t%s\t%s\n\tMedian:\t%s\t%s\n\tStddev:\t%s\t%s\n\tMin:\t%s\t%s\n\tMax:\t%s\t%s\n" % (Device,Iterations,np.log10(np.mean(Squertz)),np.mean(Squertz),np.log10(np.median(Squertz)),np.median(Squertz),np.log10(np.std(Squertz)),np.std(Squertz),np.log10(np.min(Squertz)),np.min(Squertz),np.log10(np.max(Squertz)),np.max(Squertz)))

    clDataX.release()

    clDataV.release()

    clKinetic.release()

    clPotential.release()