density, pressure, viscosity calculation method

In SPH method, we basically calculate each particle's pressure and viscosity, then using the gradient of the pressure and the Laplacian of the viscosity to get the force at each particle, once we get the force at each particle, we can then get the acceleration, velocity and ultimately, using the velocity to update each particles' position.

After I studied the paper mentioned in the last post, I believed the CPU serial implementation of SPH method should be quite straight forward.
Either density, pressure or viscosity can be calculated by

and this equation's gradient, Laplacian.

the gradient or Laplacian of this equation will only affect the smoothing kernel function W(r - rj,h).

To implement this equation, a naive method is to use two  for loops.
for example, if we want to calculate the density:

for ( data1 = particle.data; dat1 < dat1_end; dat1 += particle.width ) {

p = data1;

sum = 0.0;

data2_end = particle.data + NUM_PARTICLES*particle.width;

for ( data2 = particle.data; dat2 < dat2_end; dat2 += particle.width ) {

q = data2;

if ( p==q ) continue;

dx = ( p->pos.x - q->pos.x)*d;

dy = ( p->pos.y - q->pos.y)*d;

dz = ( p->pos.z - q->pos.z)*d;

d2 = (dx*dx + dy*dy + dz*dz);

if ( sRadius2> d2 ) {

c =  sRadius2 - d2;

sum += c * c * c;

}

}

p->density = sum * PARTICLE_MASS* kernFun1 ;