Note that there are some explanatory texts on larger screens.

plurals
  1. POAdvice on optimization floating point calculations: solving a quadratic equation
    text
    copied!<p>I'm working on a granular dynamics problem. The computationally expensive part is the function below that solve a quadratic equation to detect the collision of two particles. </p> <p>I was wondering if this can be easily optimized or if I'm doing something noticeably stupid? For example, is it a good idea to use these <code>const double x1 = p1-&gt;x;</code> constructions to give the compiler a hint? Looking at the assembler code, the compiler uses SSE instructions, but I have no idea if they are optimal in any way (probably not). According to the profiler, most of the time is spend in calculating the expressions <code>a</code>, <code>b</code> and <code>c</code>. What do you (in general) do when you are trying to optimize some kernel function like the one below? </p> <pre><code>void detect_collision_of_pair(struct particle* p1, struct particle* p2){ const double x1 = p1-&gt;x; const double y1 = p1-&gt;y; const double z1 = p1-&gt;z; const double x2 = p2-&gt;x; const double y2 = p2-&gt;y; const double z2 = p2-&gt;z; const double vx1 = p1-&gt;vx; const double vy1 = p1-&gt;vy; const double vz1 = p1-&gt;vz; const double vx2 = p2-&gt;vx; const double vy2 = p2-&gt;vy; const double vz2 = p2-&gt;vz; const double a = vx1*vx1 - 2.*vx1*vx2 + vx2*vx2 + vy1*vy1 - 2.*vy1*vy2 + vy2*vy2 + vz1*vz1 - 2.*vz1*vz2 + vz2*vz2; const double b = 2.*vx1*x1 - 2.*vx2*x1 - 2.*vx1*x2 + 2.*vx2*x2 + 2.*vy1*y1 - 2.*vy2*y1 - 2.*vy1*y2 + 2.*vy2*y2 + 2.*vz1*z1 - 2.*vz2*z1 - 2.*vz1*z2 + 2.*vz2*z2; const double c = -4.*particle_radius*particle_radius + x1*x1 - 2.*x1*x2 + x2*x2 + y1*y1 - 2.*y1*y2 + y2*y2 + z1*z1 - 2.*z1*z2 + z2*z2; double root = b*b-4.*a*c; if (root&gt;=0.){ root = sqrt(root); double time1 = (-b-root)/(2.*a); double time2 = (-b+root)/(2.*a); if ( (time1&gt;-dt &amp;&amp; time1&lt;0.) || (time1&lt;-dt &amp;&amp; time2&gt;0) ){ double times = -dt; if (time1&gt;-dt || time2&lt;0){ if (time1&gt;-dt){ times = time1; }else{ times = time2; } } resolve_collision(p1,p2,times); } } } </code></pre>
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload