GCC Optimization.

GCC Optimization.

[Manan]

Hi all,
    I am a Physics student.  We in our research group have been 
developing a molecular simulation software.
So to make our code run fast we are using gcc optimizations but i have 
got a case in which gcc shows different behaviour when run with or 
without optimization.
And i am not very sure how to correct this situation.

Below is the part of code in which i am seeing this behaviour. Red line 
is the exact line which gives different results with and without 
optimization.
The value of  x[i] at that point is 1.5(exact) value of box_xh is  2.1 
and rnx is 0.6 so icelx should evaluate to 6. which it does if i dont 
use optimization but when i switch optimization it gives icelx = 5. 
which is not correct.

So can some please tell me what i am doing wrong in this code.
I will really appreciate your help.
Thanks
-Manan
 


void make_clist(cl *clist,cord position,sim simul,syst sys,out 
*output,files all_files)
{
    //This code makes neighbour list
    real box_x = simul.box_x;
    real box_y = simul.box_y;
    real box_z = simul.box_z;

    real box_xh = simul.box_xh;
    real box_yh = simul.box_yh;
    real box_zh = simul.box_zh;

    int nx,ny,nz;
    nx = clist->nx;
    ny = clist->ny;
    nz = clist->nz;

    real rnx = box_x / nx;//width of cell in x direction
    real rny = box_y / ny;//width of cell in y direction
    real rnz = box_z / nz;//width of cell in z direction

    clist->rnx = rnx;
    clist->rny = rny;
    clist->rnz = rnz;


    int icel,i,icelx,icely,icelz,ncell;
    real *x,*y,*z;

    int *hoc,*fcell_list,*rcell_list;
    x = position.x;
    y = position.y;
    z = position.z;

    hoc = clist->hoc;
    fcell_list = clist->fcell_list;
    rcell_list = clist->rcell_list;
    ncell = clist->ncell;

    for(icel = 0;icel < ncell;icel++) hoc[icel] = -1;//All the cells 
have -1 on there top initially
    for(i = 0;i < simul.nsites;i++) fcell_list[i] = rcell_list[i] = 
-1;//Initialize the neighbour list
    for( i = 0;i < simul.nsites;i++)
    {
        icelx = floor((x[i] + box_xh) / rnx);//This is because the cells 
span fom -box_h to box_h
        icely = floor((y[i] + box_yh) / rny);
        icelz = floor((z[i] + box_zh) / rnz);
        icel = icelx + icely * nx + icelz * nx * ny;
        //printf("%d\n",icelx);

#ifdef DEBUG2
        output->log_ptr += sprintf(output->log_ptr,"%d\t",icel+1);
        ioflush(output,all_files,false);
#endif

        fcell_list[i] = hoc[icel];
        if(hoc[icel] != -1) rcell_list[hoc[icel]] = i;
        hoc[icel] = i;
    }

#ifdef DEBUG2
    output->log_ptr += sprintf(output->log_ptr,"Forward nlist : ");
    ioflush(output,all_files,false);

    for(i = 0;i < simul.nsites;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",fcell_list[i]+1);
    ioflush(output,all_files,false);
   
    output->log_ptr += sprintf(output->log_ptr,"\n");
    ioflush(output,all_files,false);
   
    output->log_ptr += sprintf(output->log_ptr,"Backward nlist : ");
    ioflush(output,all_files,false);
   
    for(i = 0;i < simul.nsites;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",rcell_list[i]+1);
    ioflush(output,all_files,false);
   
    output->log_ptr += sprintf(output->log_ptr,"\n");
    ioflush(output,all_files,false);
   
    output->log_ptr += sprintf(output->log_ptr,"hoc list :");
    ioflush(output,all_files,false);
    for(i = 0;i < ncell;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",hoc[i]+1);
    ioflush(output,all_files,false);
   
    output->log_ptr += sprintf(output->log_ptr,"\n");
    ioflush(output,all_files,false);

#endif           
}

Re: GCC Optimization.

[Nathan Sidwell]

Manan wrote:

> Below is the part of code in which i am seeing this behaviour. Red line 
> is the exact line which gives different results with and without 
> optimization.
What red line?

> The value of  x[i] at that point is 1.5(exact) value of box_xh is  2.1 
> and rnx is 0.6 so icelx should evaluate to 6. which it does if i dont 
> use optimization but when i switch optimization it gives icelx = 5. 
> which is not correct.
> 
> So can some please tell me what i am doing wrong in this code.
You are presuming that the computer represents decimal floating point
values exactly.  Refer to

What Every Computer Scientist Should Know about Floating-Point Arithmetic  by 
David Goldberg, including Doug Priest's supplement (PDF format)

linked to from the gcc 'further readings' page.

nathan
-- 
Nathan Sidwell    ::   http://www.codesourcery.com   ::     CodeSourcery LLC
nathan@codesourcery.com    ::     http://www.planetfall.pwp.blueyonder.co.uk

Re: GCC OPTIMIZATION

[Ian Lance Taylor]

Manan Chopra <mchopra@wisc.edu> writes:

>    I am a Physics student.  We in our research group have been
> developing a molecular simulation software.
> So to make our code run fast we are using gcc optimizations but i have
> got a case in which gcc shows different behaviour when run with or
> without optimization.
> And i am not very sure how to correct this situation.
> 
> Below is the part of code in which i am seeing this behaviour. THIS IS
> BUGGY LINE is the exact line which gives different results with and
> without optimization.
> The value of  x[i] at that point is 1.5(exact) value of box_xh is  2.1
> and rnx is 0.6 so icelx should evaluate to 6. which it does if i dont
> use optimization but when i switch optimization it gives icelx =
> 5. which is not correct. One moe thing i should mention real is double
> in my code.

You neglected to mention what type of machine you are running on.  I
would guess that you are using some sort of x86 architecture.  On the
x86 the type double is 64 bits, but floating point computation is done
by default in 80 bit registers.  The difference in precision can cause
unexpected effects.

In particular neither 0.6 nor 2.1 can be precisely represented in a
binary floating point format.  When compiling without optimization,
the values are forced back to 64-bits on the stack at each step.  When
compiling with optimization this does not happen.  With optimization
the computation using values in 80-bit precision yields a value just
smaller than 6.  When you then call floor, this is reduced to 5.

The quick workaround is to use the -ffloat-store option, which will
probably fix your immediate problem.  Or if you have a sufficiently
new processor and compiler you could try -mfpmath=sse.  Medium term
you should avoid discontinuous functions like floor and you should
avoid comparing floating point values for equality--instead check
whether they differ by some appropriately small epsilon.  Longer term
I would recommend a better understanding of computer approximations of
floating point arithmetic.

Ian

GCC OPTIMIZATION

[Manan Chopra]

Hi all,
   I am a Physics student.  We in our research group have been 
developing a molecular simulation software.
So to make our code run fast we are using gcc optimizations but i have 
got a case in which gcc shows different behaviour when run with or 
without optimization.
And i am not very sure how to correct this situation.

Below is the part of code in which i am seeing this behaviour. THIS IS 
BUGGY LINE is the exact line which gives different results with and 
without optimization.
The value of  x[i] at that point is 1.5(exact) value of box_xh is  2.1 
and rnx is 0.6 so icelx should evaluate to 6. which it does if i dont 
use optimization but when i switch optimization it gives icelx = 5. 
which is not correct. One moe thing i should mention real is double in 
my code.

So can some please tell me what i am doing wrong in this code.
I will really appreciate your help.
Thanks
-Manan



void make_clist(cl *clist,cord position,sim simul,syst sys,out 
*output,files all_files)
{
   //This code makes neighbour list
   real box_x = simul.box_x;
   real box_y = simul.box_y;
   real box_z = simul.box_z;

   real box_xh = simul.box_xh;
   real box_yh = simul.box_yh;
   real box_zh = simul.box_zh;

   int nx,ny,nz;
   nx = clist->nx;
   ny = clist->ny;
   nz = clist->nz;

   real rnx = box_x / nx;//width of cell in x direction
   real rny = box_y / ny;//width of cell in y direction
   real rnz = box_z / nz;//width of cell in z direction

   clist->rnx = rnx;
   clist->rny = rny;
   clist->rnz = rnz;


   int icel,i,icelx,icely,icelz,ncell;
   real *x,*y,*z;

   int *hoc,*fcell_list,*rcell_list;
   x = position.x;
   y = position.y;
   z = position.z;

   hoc = clist->hoc;
   fcell_list = clist->fcell_list;
   rcell_list = clist->rcell_list;
   ncell = clist->ncell;

   for(icel = 0;icel < ncell;icel++) hoc[icel] = -1;//All the cells have 
-1 on there top initially
   for(i = 0;i < simul.nsites;i++) fcell_list[i] = rcell_list[i] = 
-1;//Initialize the neighbour list
   for( i = 0;i < simul.nsites;i++)
   {
       icelx = floor((x[i] + box_xh) / rnx);//This is because the cells 
span fom -box_h to box_h //THIS IS BUGY LINE
       icely = floor((y[i] + box_yh) / rny);
       icelz = floor((z[i] + box_zh) / rnz);
       icel = icelx + icely * nx + icelz * nx * ny;
       //printf("%d\n",icelx);

#ifdef DEBUG2
       output->log_ptr += sprintf(output->log_ptr,"%d\t",icel+1);
       ioflush(output,all_files,false);
#endif

       fcell_list[i] = hoc[icel];
       if(hoc[icel] != -1) rcell_list[hoc[icel]] = i;
       hoc[icel] = i;
   }

#ifdef DEBUG2
   output->log_ptr += sprintf(output->log_ptr,"Forward nlist : ");
   ioflush(output,all_files,false);

   for(i = 0;i < simul.nsites;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",fcell_list[i]+1);
   ioflush(output,all_files,false);
     output->log_ptr += sprintf(output->log_ptr,"\n");
   ioflush(output,all_files,false);
     output->log_ptr += sprintf(output->log_ptr,"Backward nlist : ");
   ioflush(output,all_files,false);
     for(i = 0;i < simul.nsites;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",rcell_list[i]+1);
   ioflush(output,all_files,false);
     output->log_ptr += sprintf(output->log_ptr,"\n");
   ioflush(output,all_files,false);
     output->log_ptr += sprintf(output->log_ptr,"hoc list :");
   ioflush(output,all_files,false);
   for(i = 0;i < ncell;i++) output->log_ptr += 
sprintf(output->log_ptr,"%d  ",hoc[i]+1);
   ioflush(output,all_files,false);
     output->log_ptr += sprintf(output->log_ptr,"\n");
   ioflush(output,all_files,false);

#endif           }