Asked by Francesco
on 31 Jul 2013

Hello,

I have a Matlab version and MEX C version of a function and it works fine. I tried implementing a CUDA C version and compiling it using nvcc into PTX (no problems here) and using the GPU kernel Matlab function in the Parallel Computing Toolbox, but I get NaN values as outputs for my inputs (for the Matlab native and MEX C version I get the correct outputs).

I have tested some example CUDA codes adding vectors etc. and they seem to work OK.

Here is my procedure:

kernel = parallel.gpu.CUDAKernel('StressDueToSeg.ptx','StressDueToSeg.cu');

kernel.ThreadBlockSize = 1024;

kernel.GridSize = ceil(N/1024);

[~,~,~,~,~,~,~,~,~,~,~,~,...

s11, s22, s33, s12, s13, s23] = feval(kern,N,S,px,py,pz,p1x,p1y,p1z,p2x,p2y,p2z,bx,by,bz,a,NU,MU,...

sxx,syy,szz,sxy,syz,sxz);

The inputs are scalars 'N,S,a,MU,NU', 1D vectors of length N 'px,py,pz', and 1D vectors of length S 'p1x,p1y,p1z,p2x,p2y,p2z,bx,by,bz'.

My .cu code is as follows:

#include <math.h>

__global__ void StressDueToSeg(int N, int S, double *px, double *py, double *pz,

double *p1x, double *p1y, double *p1z,

double *p2x, double *p2y, double *p2z,

double *bx, double *by, double *bz,

const double a, const double MU, const double NU,

double *sxx, double *syy, double *szz,

double *sxy, double *syz, double *sxz)

{

double oneoverLp, common;

double vec1x, vec1y, vec1z;

double tpx, tpy, tpz;

double Rx, Ry, Rz, Rdt;

double ndx, ndy, ndz;

double d2, s1, s2, a2, a2_d2, a2d2inv;

double Ra, Rainv, Ra3inv, sRa3inv;

double s_03a, s_13a, s_05a, s_15a, s_25a;

double s_03b, s_13b, s_05b, s_15b, s_25b;

double s_03, s_13, s_05, s_15, s_25;

double m4p, m8p, m4pn, mn4pn, a2m8p;

double txbx, txby, txbz;

double dxbx, dxby, dxbz;

double dxbdt, dmdxx, dmdyy, dmdzz, dmdxy, dmdyz, dmdxz;

double tmtxx, tmtyy, tmtzz, tmtxy, tmtyz, tmtxz;

double tmdxx, tmdyy, tmdzz, tmdxy, tmdyz, tmdxz;

double tmtxbxx, tmtxbyy, tmtxbzz, tmtxbxy, tmtxbyz, tmtxbxz;

double dmtxbxx, dmtxbyy, dmtxbzz, dmtxbxy, dmtxbyz, dmtxbxz;

double tmdxbxx, tmdxbyy, tmdxbzz, tmdxbxy, tmdxbyz, tmdxbxz;

double I_03xx, I_03yy, I_03zz, I_03xy, I_03yz, I_03xz;

double I_13xx, I_13yy, I_13zz, I_13xy, I_13yz, I_13xz;

double I_05xx, I_05yy, I_05zz, I_05xy, I_05yz, I_05xz;

double I_15xx, I_15yy, I_15zz, I_15xy, I_15yz, I_15xz;

double I_25xx, I_25yy, I_25zz, I_25xy, I_25yz, I_25xz;

int seg;

int coord = threadIdx.x + blockIdx.x * blockDim.x ;

if (coord<N) {

//precompute some constants

m4p = 0.25 * MU / M_PI;

m8p = 0.5 * m4p;

m4pn = m4p / (1 - NU);

mn4pn = m4pn * NU;

a2 = a * a;

a2m8p = a2 * m8p;

for (seg=0; seg<S; seg++) { //loop over the segments

vec1x = p2x[seg] - p1x[seg];

vec1y = p2y[seg] - p1y[seg];

vec1z = p2z[seg] - p1z[seg];

oneoverLp = 1 / sqrt(vec1x*vec1x + vec1y*vec1y + vec1z*vec1z);

tpx = vec1x * oneoverLp;

tpy = vec1y * oneoverLp;

tpz = vec1z * oneoverLp;

Rx = px[coord] - p1x[seg];

Ry = py[coord] - p1y[seg];

Rz = pz[coord] - p1z[seg];

Rdt = Rx*tpx + Ry*tpy + Rz*tpz;

ndx = Rx - Rdt*tpx;

ndy = Ry - Rdt*tpy;

ndz = Rz - Rdt*tpz;

d2 = ndx*ndx + ndy*ndy + ndz*ndz;

s1 = -Rdt;

s2 = -((px[coord]-p2x[seg])*tpx + (py[coord]-p2y[seg])*tpy + (pz[coord]-p2z[seg])*tpz);

a2_d2 = a2 + d2;

a2d2inv = 1 / a2_d2;

Ra = sqrt(a2_d2 + s1*s1);

Rainv = 1 / Ra;

Ra3inv = Rainv * Rainv * Rainv;

sRa3inv = s1 * Ra3inv;

s_03a = s1 * Rainv * a2d2inv;

s_13a = -Rainv;

s_05a = (2*s_03a + sRa3inv) * a2d2inv;

s_15a = -Ra3inv;

s_25a = s_03a - sRa3inv;

Ra = sqrt(a2_d2 + s2*s2);

Rainv = 1 / Ra;

Ra3inv = Rainv * Rainv * Rainv;

sRa3inv = s2 * Ra3inv;

s_03b = s2 * Rainv * a2d2inv;

s_13b = -Rainv;

s_05b = (2*s_03b + sRa3inv) * a2d2inv;

s_15b = -Ra3inv;

s_25b = s_03b - sRa3inv;

s_03 = s_03b - s_03a;

s_13 = s_13b - s_13a;

s_05 = s_05b - s_05a;

s_15 = s_15b - s_15a;

s_25 = s_25b - s_25a;

txbx = tpy*bz[seg] - tpz*by[seg];

txby = tpz*bx[seg] - tpx*bz[seg];

txbz = tpx*by[seg] - tpy*bx[seg];

dxbx = ndy*bz[seg] - ndz*by[seg];

dxby = ndz*bx[seg] - ndx*bz[seg];

dxbz = ndx*by[seg] - ndy*bx[seg];

dxbdt = dxbx*tpx + dxby*tpy + dxbz*tpz;

dmdxx = ndx * ndx;

dmdyy = ndy * ndy;

dmdzz = ndz * ndz;

dmdxy = ndx * ndy;

dmdyz = ndy * ndz;

dmdxz = ndx * ndz;

tmtxx = tpx * tpx;

tmtyy = tpy * tpy;

tmtzz = tpz * tpz;

tmtxy = tpx * tpy;

tmtyz = tpy * tpz;

tmtxz = tpx * tpz;

tmdxx = 2 * tpx * ndx;

tmdyy = 2 * tpy * ndy;

tmdzz = 2 * tpz * ndz;

tmdxy = tpx*ndy + tpy*ndx;

tmdyz = tpy*ndz + tpz*ndy;

tmdxz = tpx*ndz + tpz*ndx;

tmtxbxx = 2 * tpx * txbx;

tmtxbyy = 2 * tpy * txby;

tmtxbzz = 2 * tpz * txbz;

tmtxbxy = tpx*txby + tpy*txbx;

tmtxbyz = tpy*txbz + tpz*txby;

tmtxbxz = tpx*txbz + tpz*txbx;

dmtxbxx = 2 * ndx * txbx;

dmtxbyy = 2 * ndy * txby;

dmtxbzz = 2 * ndz * txbz;

dmtxbxy = ndx*txby + ndy*txbx;

dmtxbyz = ndy*txbz + ndz*txby;

dmtxbxz = ndx*txbz + ndz*txbx;

tmdxbxx = 2 * tpx * dxbx;

tmdxbyy = 2 * tpy * dxby;

tmdxbzz = 2 * tpz * dxbz;

tmdxbxy = tpx*dxby + tpy*dxbx;

tmdxbyz = tpy*dxbz + tpz*dxby;

tmdxbxz = tpx*dxbz + tpz*dxbx;

common = m4pn * dxbdt;

I_03xx = common + m4pn*dmtxbxx - m4p*tmdxbxx;

I_03yy = common + m4pn*dmtxbyy - m4p*tmdxbyy;

I_03zz = common + m4pn*dmtxbzz - m4p*tmdxbzz;

I_03xy = m4pn*dmtxbxy - m4p*tmdxbxy;

I_03yz = m4pn*dmtxbyz - m4p*tmdxbyz;

I_03xz = m4pn*dmtxbxz - m4p*tmdxbxz;

I_13xx = -mn4pn * tmtxbxx;

I_13yy = -mn4pn * tmtxbyy;

I_13zz = -mn4pn * tmtxbzz;

I_13xy = -mn4pn * tmtxbxy;

I_13yz = -mn4pn * tmtxbyz;

I_13xz = -mn4pn * tmtxbxz;

I_05xx = common*(a2+dmdxx) - a2m8p*tmdxbxx;

I_05yy = common*(a2+dmdyy) - a2m8p*tmdxbyy;

I_05zz = common*(a2+dmdzz) - a2m8p*tmdxbzz;

I_05xy = common*dmdxy - a2m8p*tmdxbxy;

I_05yz = common*dmdyz - a2m8p*tmdxbyz;

I_05xz = common*dmdxz - a2m8p*tmdxbxz;

I_15xx = a2m8p*tmtxbxx - common*tmdxx;

I_15yy = a2m8p*tmtxbyy - common*tmdyy;

I_15zz = a2m8p*tmtxbzz - common*tmdzz;

I_15xy = a2m8p*tmtxbxy - common*tmdxy;

I_15yz = a2m8p*tmtxbyz - common*tmdyz;

I_15xz = a2m8p*tmtxbxz - common*tmdxz;

I_25xx = common * tmtxx;

I_25yy = common * tmtyy;

I_25zz = common * tmtzz;

I_25xy = common * tmtxy;

I_25yz = common * tmtyz;

I_25xz = common * tmtxz;

sxx[coord] += I_03xx*s_03 + I_13xx*s_13 + I_05xx*s_05 +

I_15xx*s_15 + I_25xx*s_25;

syy[coord] += I_03yy*s_03 + I_13yy*s_13 + I_05yy*s_05 +

I_15yy*s_15 + I_25yy*s_25;

szz[coord] += I_03zz*s_03 + I_13zz*s_13 + I_05zz*s_05 +

I_15zz*s_15 + I_25zz*s_25;

sxy[coord] += I_03xy*s_03 + I_13xy*s_13 + I_05xy*s_05 +

I_15xy*s_15 + I_25xy*s_25;

syz[coord] += I_03yz*s_03 + I_13yz*s_13 + I_05yz*s_05 +

I_15yz*s_15 + I_25yz*s_25;

sxz[coord] += I_03xz*s_03 + I_13xz*s_13 + I_05xz*s_05 +

I_15xz*s_15 + I_25xz*s_25;

}

}

return;

}

I am new to CUDA so perhaps it's an implementation problem; however, the MEX C code is essentially very similar (it has a 2nd loop with counter 'coord', rather than coord = threadIdx.x + blockIdx.x * blockDim.x) and that works OK.

Please enlighten me.

Thank you

F

Answer by Ben Tordoff
on 1 Aug 2013

Accepted Answer

First thing to do to track this down is to clearly specify your inputs and outputs. Inputs should be const. Based on your description, your function should probably look like:

__global__ void StressDueToSeg(int N, int S,

double const *px, double const *py, double const *pz,

double const *p1x, double const *p1y, double const *p1z,

double const *p2x, double const *p2y, double const *p2z,

double const *bx, double const *by, double const *bz,

const double a, const double MU, const double NU,

double *sxx, double *syy, double *szz,

double *sxy, double *syz, double *sxz)

with that done, you only need to capture the outputs (final 6 arrays). This is both easier to read and will probably be significantly faster:

[s11, s22, s33, s12, s13, s23] = feval( ...

kern, N, S, ...

px, py, pz, ...

p1x, p1y, p1z, ...

p2x, p2y, p2z, ...

bx, by, bz, ...

a, NU, MU, ...

sxx, syy, szz, sxy, syz, sxz);

(NB: it kind-of looks like you have the last two outputs switched). I have tried this with some dummy values of MU, NU, A etc. and don't get NaNs. There are two important considerations though:

- You need to get your input arguments in the right order (you appear to have put NU first in the call site, but MU first in the kernel)
- You need to intialize the output arrays to zero before calling the kernel since the kernel uses += without initializing the array. (Even better would be for the kernel to zero the values at the start.)
- Passing NU=1.0 will guarantee inf/nan everywhere since you have the line m4p / (1-NU).

You can debug this kind of thing by putting the intermediate values into the output arguments (possibly additional ones added for the purpose) - that's what I did to realise that MU and NU were swapped.

I hope that helps.

Ben

Francesco
on 1 Aug 2013

Good morning.

Thank you so much for you input. It was indeed that I had switched MU and NU around when calling the kernel.

I guess this leads to the question: is there a debugger that can go inside the kernel functions that Matlab is calling? I am not aware that the Matlab debugger can do this, right?

Thank you,

You helped a lot! :)

F

One more question: when initializing sxx,syy..etc, I use parallel.gpu.GPUArray.zeros(N,1) to create them on the GPU and store them in device memory. However, I could, as you said, initialize 'on the fly' in the kernel, i.e. before the S++ loop have:

sxx[coord] = 0;

If I do this, what should I pass as sxx,syy etc. in the kernel function in Matlab?

Thanks again!

Ben Tordoff
on 1 Aug 2013

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.