restrict in generic kernels

[ldobsik]

Hello, have you considered using the restrict qualifier on some of the output parameters in the generic kernels? It seems to unlock auto-vectorization e.g. on gcc for RISC-V /w V-ext.

[giaf]

Hi, I experimented with restrict but admittedly it was quite a few years ago, and back then it didn't help, at least for the cases I looked into.
Could you report here a few examples of the intended changes that allowed to enable auto-vectorization in gcc for RISC-V?

[ldobsik]

E.g. the generic kernel_sgemm_nn_4x4_lib4. You can compare the output with and without restrict.

[giaf]

Strictly speaking, D can (and often does) alias with C, and I could see ways to write a vectorized piece of code that is correct if C and D don't alias, and wrong if they do (but admittedly such code would be a bit artificial as it would make useless operations, so not likely to happen in practice).
But still, in theory restrict should not be used in such case.

I had a look at the assembly code, but I'm not familiar with RISC-V. What is the vector size in this case?
In the vectorized code there seem to be a bunch of gather instructions compared to the amount of fmadd ones, I don't see why this should be the case.
fmadd itself seems to be a fused multiply add variant that overwrited a multiplicand, while I would find much more natural a variant that overwrites an addend.
Did you have a chance to benchmark if the vectorized code (i.e. with restrict) is faster than the scalar one (i.e. without restrict)?

[ldobsik]

Hello,
yes, if anything the situation around parameter aliasing could use some clarification (i.e. what can and what shall not be aliased).

The vector size in the example is 256 bits, i.e. 8 floats. I didn't really analyze the generated code, just noticed it is using the vector instruction after adding restrict.
I'm guessing the the gathers are used to reshuffle the data correctly after using the 8-wide loads. Maybe strided loads could also be used effectively instead.
Regarding actual performance - I did only basic instruction-accurate simulation of a similar code some time ago and I remember the vector code had a lower number of instructions executed. No actual execution time benchmark, though.

[giaf]

In the current implementation of all linear algebra routines in BLASFEO, in the gemm kernels and the like either C and D point to the same location (say in case of destructive operation C = A*B + C), or the pieces of memory accessed don't overlap at all (say in case of non-destructive operation D = A*B + C).

[ldobsik]

Thank you for the clarification! I assume, the A and B can alias between each other while they cannot alias with neither C or D.

[giaf]

Yes A and B often alias (e.g. implementation of symmetric routines).

About other cases like aliasing A or B with C, this never happens in the BLASFEO implementation itself but an user could call D = A*A + A and expect a correct result.

If D is aliased with either factor A or B, it can't be expected a correct result, because the factor would change in implementation-dependent ways as D gets computed.
Anyway, also in the official BLAS routine dgemm there is nothing explicitly preventing such (mis-)use.