Checking for NaN and -1.#IND

Accelerated Computing CUDA CUDA Programming and Performance
1

Hey guys,

I have a program im working with and currently have up and running.
However the results I am getting back are either incredibly small, NaN, or -1.#IND.

I know the results I SHOULD be getting since we have a VB version of the program running with the same data.

I can’t seem to track down where the error is coming from exactly so I want to put in checks to catch where the NaN/-1.#IND is coming from.

I’ve seen suggestions to use isnan() but that is not recognized within a kernel. Could someone suggest a good way to do this within a kernel?

Thanks a lot!

2

isnan should work in kernels and is what you want, can you show your code?

3

Hmm that’s odd, I should only need to #include math.h right?

Maybe I did something wrong then, I will check again.

My code is pretty large out of where it could be going wrong. But So far I have tracked the error down to one function.

__device__ double getInteraction(int isoN, int run, double xtemp, double ytemp, double ztemp, double thxt, double thyt, double thzt, 

									double Isomer1_arr[100][17], double Isomer2_arr[100][17], double SP_arr[100][17], double Shift_arr[100][17]){

//rotation variables

double cx, cy, cz;									//cos xyz

double sx, sy, sz;									//sin xyz

double xt, yt, zt;									//temp/prime

double R11, R12, R13, R21, R22, R23, R31, R32, R33;	//temp values

double thx, thy, thz;								//shift angle values

double xshift, yshift, zshift;						//x,y,z shift values

int iMax, jMax;

double x, y, z;										//holds values from isomer 1 or 2

double temp_arr[100][17];							//temp array (replaces sheet 5) 100 rows of 17 columns

//energy variables

double xs, ys, zs;									//coordinates of SP, reuse x,y,z from before for coord of isomers

double r, r2;										//Inter atome (mol-SP) distance and square of same

double Eq, Eqt;										//electrostatic energy

double Enb, Enbt;									//non-bonded energy			

double A, B, rStar, epsilon;						//non-bonded parameters

double rStarS, epsilonS;							//non-bonded parameters

int line;											//non-bonded parameters

double Et;											//total energy

double Ett;	

double q, qs;										//charges for isomer and SP

		

//////////////////////ROTATION/////////////////////////

//check if first run or not, if first use shift_arr[][] else use new xyz, txyz from input

//these first values should NOT be set, they should be obtained somehow **fix**

if(run == 1){

thx = -66.25976621 * pi / 180.000;

thy = 120.893644 * pi / 180.000;

thz = 261.891781 * pi / 180.000;

xshift = 0.18333688;

yshift = 9.258853604;

zshift = -0.041139504;

}

else{

	thx = thxt * pi / 180.000;

	thy = thyt * pi / 180.000;

	thz = thzt * pi / 180.000;

	xshift = xtemp;

	yshift = ytemp;

	zshift = ztemp;

}

xt = 0;

yt = 0;

zt = 0;

//do trig

cx = cos(thx);

sx = sin(thx);

cy = cos(thy);

sy = sin(thy);

cz = cos(thz);

sz = sin(thz);

//get matrix elements

R11 = cy * cz;

R12 = cx * sz + sx * sy * cz;

R13 = sx * sz - cx * sy * cz;

R21 = -cy * sz;

R22 = cx * cz - sx * sy * sz;

R23 = sx * cz + cx * sy * sz;

R31 = sy;

R32 = -sx * cy;

R33 = cx * cy;

//get number of atoms in Isomer1 (iMax) and stationary phase (jMax)

if(isoN == 1)								

	iMax = Isomer1_arr[0][2];			//num rows [0,2] Isomer1_arr ----- iMax = Isomer1_arr_1D[0*100+2];	

else

	iMax = Isomer2_arr[0][2];			//num rows [0,2] Isomer2_arr ----- iMax = Isomer2_arr_1D[0*100+2];

jMax = 17;									//num cols

//copy isomer1 or 2 to temp_arr (sheet 5)

for(int i = 0; i < iMax; i++){				//Y

	for(int j = 0; j < jMax; j++){			//X

		if(isoN == 1)

			temp_arr[i][j] = Isomer1_arr[i][j];	//[i,j]

		else

			temp_arr[i][j] = Isomer2_arr[i][j];	//[i,j]

	}

}

//make transformation

for(int i = 2; i < iMax + 2; i++){

	if(isoN == 1){

		x = Isomer1_arr[i][2];						//[i,2]

		y = Isomer1_arr[i][3];						//[i,3]

		z = Isomer1_arr[i][4];						//[i,4]

	}

	else{

		x = Isomer2_arr[i][2];						//[i,2]

		y = Isomer2_arr[i][3];						//[i,3]

		z = Isomer2_arr[i][4];						//[i,4]

	}

	xt = x * R11 + y * R12 + z * R13;

	yt = x * R21 + y * R22 + z * R23;

	zt = x * R31 + y * R32 + z * R33;

	temp_arr[i][2] = xt + xshift;						//[i,2]

	temp_arr[i][3] = yt + yshift;						//[i,3]

	temp_arr[i][4] = zt + zshift;						//[i,4]

}

////////////////////END ROTATION --> ENERGY///////////////////////

//reset x,y,z values to 0.0 for reuse

x = 0.0;

y = 0.0;

z = 0.0;

//get number of atoms in Isomer 1 and SP

iMax = Isomer1_arr[0][2];								//[0,2]

jMax = SP_arr[0][2];									//[0,2]

//initialize energies

Eqt = 0.0;											//sum of the electrostatic interactions to date					

Enbt = 0.0;											//sum of the non-bonded (6-12) interactions to date

Ett = 0.0;											//sum of the total energies to date

line = 2;

for(int i = 2; i < iMax + 2; i++){

	x = temp_arr[i][2];								//x for isomer		//[i + 2, 2]

	y = temp_arr[i][3];								//y for isomer		//[i + 2, 3]

	z = temp_arr[i][4];								//z for isomer		//[i + 2, 4]

	rStar = temp_arr[i][5];							//rStar for isomer	//[i + 2, 5]

	epsilon = temp_arr[i][6];						//epsilon for isomer//[i + 2, 6]

	q = temp_arr[i][7];								//charge for isomer	//[i + 2, 7]

	for(int j = 2; j < jMax + 2; j++){

		xs = SP_arr[j][2];							//x for SP			//[j + 2, 2]

		ys = SP_arr[j][3];							//y for SP			//[j + 2, 3]

		zs = SP_arr[j][4];							//z for SP			//[j + 2, 4]

		rStarS = SP_arr[j][5];						//rStar for SP		//[j + 2, 5]

		epsilonS = SP_arr[j][6];					//epsilon for SP	//[j + 2, 6]

		qs = SP_arr[j][7];							//charge for SP		//[j + 2, 7]

		//calculate distance from atoms on SP

		r2 = ((x - xs) * (x - xs)) + ((y - ys) * (y - ys)) + ((z - zs) * (z - zs));

		r = sqrt(r2);								//distance between atome j on stationary phase and atom i on isomer

		//printf(i);

		//printf(j);

		//printf(r);

		if(r2 == 0.0){

			Eq = 1000000.0;

		}

		else{

			Eq = q * qs / r * 0.529177249 * 627.509;	//A to Bohr gives hartrees which then needs conversion to kcal/mol

		}

		Eqt = Eqt + Eq;

		//printf(Eq);

		//calculate non bonded (6-12) potentioals and ouput

		if(r == 0.0){

			Enb = 1000000.0;

		}

		else{

			A = pow((rStar + rStarS), 12) * sqrt(epsilon * epsilonS);

			B = 2 * pow((rStar + rStarS), 6) * sqrt(epsilon * epsilonS);

			Enb =  (A / pow(r, 12)) - (B / pow(r, 6));

		}

		

		Enbt = Enbt + Enb;

		//printf(Enb);

		//calculate total interaction and output

		Et = Enb + Eq;

		Ett = Et + Ett;

		//printf(Et);

		line = line + 1;

	}

}

//printout total energy associated with interaction of isomer with stationary phase

//printf("total");

//printf(Eqt);

//printf(Enbt);

//printf(Ett);							//this is the total interaction energy

return Ett;

}

The arrays read in data in .csv files to be used. I have checked all the input and it is correct. But somewhere in here there is an error causing a NaN in the end.

4

I didn’t see an “isnan” in your code, it should be a valid function in device code.

It seems your question is really, “why is my code generating NAN values?”

One reason could be intermediate overflow. Intermediate values can overflow to +infinity or -infinity. Then when you use these values you can easily end up with infinity-infinity, which returns a NAN. In particular, in your code near the end I see some calls to “pow” calculating the 6th and 12th powers of some numbers then a few additions and subtractions. High powers with the pow function can easily overflow. I would start by looking there.

5

To slightly elaborate on the above comment: your code checks for the case r == 0, but if r is really small, then pow(r, 12) or pow(r, 6) might both still be 0 in floating point. Which would result in the subtraction of two infinities.