[Lagrangian.Solver] BlockGaussSeidelConstraintSolver: performance improvement to gaussSeidel_increment

[stu-l]

When running fallingBeamLagrangianCollision.scn

The limiting factor seems to be gaussSeidel_increment, one thread is holding up all the others. This example has 'dim' of 291 and 'nbLines' of 3. The big loop (dim) is too short for a separate threads, and 3 is an odd number. I'm hoping these are representative values.

So to code more harmoniously with the pipeline slots and SSE register sizes I've unrolled the loop in 12s (nearest multiple of SSE double word registers size 4 and the nbLines). I tried 15, but there wasn't any improvement (on my hardware).

Ideally padding the data to fit the hardware would be the best solution, even if the values were dummies that don't influence the result. Harmonious memory shaping would more than make up for any superfluous calculations.

On my hardware I get a 25% performance improvement.

[with-all-tests]

---

By submitting this pull request, I acknowledge that
I have read, understand, and agree SOFA Developer Certificate of Origin (DCO).

---

Reviewers will merge this pull-request only if

• it builds with SUCCESS for all platforms on the CI.
• it does not generate new warnings.
• it does not generate new unit test failures.
• it does not generate new scene test failures.
• it does not break API compatibility.
• it is more than 1 week old (or has fast-merge label).

[fredroy]

Not an expert, but it seems this loop is not doing the same as the current version, e.g d[j+0] do += with only k+0, k+3, k+ 6 and k+9. It seems some contributions will be missing

[stu-l]

I'm not saying you're not right, this needs more eyes and some testing. I'll write some confirmation code when I get a minute, it'll be inline rather than formal tests.

In theory though the original inner loopd[j+l] += w[j+l][k] * force[k]; iterated over l [0,2] so these are unrolled and the loop removed leaving this outer dim loop:
(j is fixed by the preceding code)

for(int k=0; k<dim; k++)
        {
//            for(unsigned int l=0; l<nb; l++)    // loop over [0,2]
//            {
                d[j+0 ] += w[j+0 ][k ] * force[k ];
                d[j+1 ] += w[j+1 ][k ] * force[k ];
                d[j+2 ] += w[j+2 ][k ] * force[k ];
                //d[j+l] += w[j+l][k] * force[k];
//           }
        }

dim turns out to be ~261 so is a relatively short, and consistent, amount (although not a power of 2, nor is 3???). I tried offloading to other threads, but this made it slower. Although I might not have got your worker thread system to work properly (why are you limiting it to half the available processors?). So the dim loop is then unrolled and the results profiled (repeatedly) giving an optimal 12 on my hardware.
However if given more obvious memory patterns, eg loop unrolling, hardware pre-fetching will get memory in advance of its use (as it says on the tin). Pre-fetching will always do as much as it can, but the more clues the developer can give it the be the better the result. Repeatedly looping over the same number also helps with branch misprediction expenses, although this is ideally supplied at compile time. Finally filling the available pipelines with as much, non predicted, code that can be executed by the chip in parallel the better. Also the compiler will start to notice and generate more optimised code, I didn't look at the generated assembler.

Changes like this will work on any architecture, it is likely to be even better on newer architecture with even more pipeline stages running even more instructions per cycle.

However; needs proper testing, I'll work on it when I get some time. Also is nb ever not 3, I have no idea what this value is? Or is dim by any chance a constant known at compile time, if it were the compiler would benefit from knowing?

[stu-l]

Not an expert, but it seems this loop is not doing the same as the current version, e.g d[j+0] do += with only k+0, k+3, k+ 6 and k+9. It seems some contributions will be missing

hm, just staring at it some more; you're right. you're saying that
d[j+0 ] += w[j+0 ][k+0 ] * force[k+0 ];
should include:

d[j+0 ] += w[j+0 ][k+0 ] * force[k+0 ];
d[j+0 ] += w[j+0 ][k+1] * force[k+1 ];
d[j+0 ] += w[j+0 ][k+2] * force[k+2 ];

I think you're right, well spotted. I was too focused on the non functional detail ;-)
I'll fix it later today.

[stu-l]

I've pushed a fix, this might make the code more readable, and maybe correct.
The performance is now only %22 better ;-)

ps. still no verification code, I'll work on that.

[fredroy]

[ci-build][with-all-tests]
This will trigger regression tests (on windows)