Here is the code:
#include <math.h>
#include <iostream>
#include <omp.h>
#include <fenv.h>
#include <signal.h>
void handler(int sig)
{
printf("Floating Point Exception\n");
exit(0);
}
const float alpha=1.5;
const unsigned int N=2;
struct Particle
{
float x, y, z;
float m;
};
Particle p[N] __attribute__((aligned(64)));
void interaction()
{
double P=0.0;
#pragma omp parallel for reduction(+:P)
for(int i=0; i<N; ++i)
{
float PP=0.0;
#pragma simd reduction(+:PP)
//#pragma novector
for(int j=0; j<N; ++j) if(i!=j)
{
float rdist1=sqrtf((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
PP+=alpha/rdist1;
}
P+=PP;
}
std::cout<<"P="<<P<<std::endl;
}
void randomizeBodies()
{
double pot_en=0.;
const double pi=3.1415926536;
float RADIUS=pow(10.0*N,1.0/3.0);
#pragma omp single
#pragma novector
for(int i=0; i<N; ++i)
{
float DISTANCE=0.0f;
if(i>0)
{
while(DISTANCE<=1.0f)
{
float theta=pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
float phi=2*pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
float rr=RADIUS*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
p[i].x =rr*sin(theta)*cos(phi);
p[i].y =rr*sin(theta)*sin(phi);
p[i].z =rr*cos(theta);
DISTANCE=10000.0f;
#pragma simd reduction(min:DISTANCE)
for(int j=0; j<i; ++j)
{
float dij=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
if(dij<DISTANCE) DISTANCE=dij;
}
}
}
else
{
float theta=pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
float phi=2*pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
float rr=RADIUS*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
p[i].x =rr*sin(theta)*cos(phi);
p[i].y =rr*sin(theta)*sin(phi);
p[i].z =rr*cos(theta);
}
}
#pragma omp parallel for reduction(+:pot_en)
for(int i=0; i<N; ++i)
{
float pp=0.0;
#pragma simd reduction(+:pp)
//#pragma novector
for(int j=0; j<N; ++j) if(i!=j)
{
float rd=sqrtf((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
pp+=alpha/rd;
}
pot_en+=pp;
}
pot_en/=2;
std::cout<<"P="<<pot_en<<std::endl;
}
int main(int argc, char **argv)
{
feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
signal(SIGFPE, handler);
randomizeBodies();
interaction();
}
I compile it using:
-DCMAKE_CXX_COMPILER=icpc -DCMAKE_CXX_FLAGS="-O2"
and with "-g" it gives the following output
Program received signal SIGFPE, Arithmetic exception.
0x00000000004024ab in randomizeBodies ()
at /home/70-gaa/source/GPU/ERROR24/error.cpp:90
90 pp+=alpha/rd;
(gdb)
at N=2,4. But at N>=8 it works. If to comment
#pragma simd reduction
and uncomment
#pragma novector
in both for loops, everything works at any N. If to use "-01", everything also works at any N.
If to to use vectorization and "-O2" or "-O3", the program throws a floating point exception at N=2,4, but at N>=8 works. Why? In principle, i need to use the following compile line:
-DCMAKE_CXX_FLAGS="-march=native -mtune=native
-ipo16 -fp-model fast=2 -O3 -qopt-report=5 -mcmodel=large".
But "-O3" doesn't work.
I work on Intel Core i7-3770 8-core CPU and use Intel compiler icpc version 17.0.1.