Towards fearless SIMD

I find intrinsics to be a useful level of abstraction because it allows me to control and automate the code generation in C++ without having to know the details of register scheduling, etc. It also means I don't have to hand write intrinsic algorithms for every architecture -- I have template libraries that can take a rough specification of the architecture characteristics and generate approximately appropriate intrinsic code. It is neither hand writing all the architecture-specific intrinsic code nor trusting the compiler to do a good job. I started working this way when I used to target code at a really weird mix of exotic vector/parallel computing hardware, and it worked nicely.

The real problem is that, in my experience, compilers universally do a poor job of auto-vectorization/parallelization unless the code pattern is blindingly obvious. All of the "higher level abstractions" I've seen are no better at finding vectorizable patterns than the compiler does -- which makes sense since they are both just software. It ends up being syntactic sugar. The amount of effort required to trick compilers into doing the right thing, plus the things they won't do at all, is such that the complexity is worse than just writing the intrinsic code directly. Writing a lot of intrinsic code directly far from optimal but fortunately those intrinsics are available in languages like C++ that can abstract much of that.

Compilers will not consistently generate competent auto-vectorized for the foreseeable future. I've seen very little forward progress in the decade I've been exposed to it. There are still some relatively simple non-vectorized code patterns that compilers currently have a difficult time detecting and doing optimized code generation for.

One problem is that although OpenCL can run on both CPU and GPU, OpenCL code that performs well on CPU and OpenCL code that performs well on GPU can be different.

_This_ is why I linked to my undefined behavior post. Your comment is about as clear an example of being in the semi-portable camp as any I've seen. And I'm not blaming you, because in C you _can't_ do SIMD in the standard camp, so it's actually one of the more compelling reasons to remain in semi-portable. Rust is different though.

Also: the linked crate does use the type system in pretty much this way so that code that clients can be safe. However, there are limitations; it's not just whether a particular instruction can be used, which remains immutable once it's detected, but also which _registers_ (and, by extension, calling convention) can be used. That varies from function to function, and requires the `#[target_feature(enable)]` annotation to control, so just having an `Avx` type in hand is not quite enough to ensure that you're in a context where using the ymm registers is ok, and the intrinsic will be inlined to a "V" variant asm instruction. This is discussed in some detail in the "caveats" section.

Ya too hard for the programmer. What it comes down to is that there are a few main categories of multiprocessing from lowest to highest level:

1. SIMD - manually deal with packed elements (like in SSE, MMX etc)

2. DSP - Maybe someone knows a better term for this, but treating each slice of the data array as an independent serial stream (shaders, OpenCL, CUDA)

3. MIMD - freeform vector/matrix operations that get compiled down to the first two categories (MATLAB, GNU Octave, Scilab)

I'm not sure if TensorFlow fits best in 2 or 3 but my gut feeling is that it's closest to 2. The problem with the lower level abstractions is that it's more work to format the data for the problem space.

So with MATLAB, everything is a vector and operations applied to each vector happen in parallel across all elements. Then if you need to, you can drop down to less efficient code and operate on elements of the vector manually with C-like code.

Unfortunately that becomes more difficult in shaders, because you can't just magically access a neighboring element. And getting general computation to work with SIMD is often infeasible because you have to rewrite your code and potentially alter the layout of the data in memory to achieve better performance.

I also agree that currently nothing really exists to give us general-purpose vector math akin to MATLAB within a language like C.

In rust this looks rather painful. Might I suggest trying this in Julia, might be a good comparison of performance, ease of use, and readability. Julia does a very nice job of compiling directly to SIMD instruction and lets you inspect the low level code generated.

inline function sin9_shaper(x) c0 = 6.28308759 c1 = -41.33318707 c2 = 81.39900205 c3 = -74.66884436 c4 = 33.15324345

    a = abs(x - round(x)) - 0.25
    a2 = a * a
    ((((a2 * c4 + c3) * a2 + c2) * a2 + c1) * a2 + c0) * a

end

function gen_sinwave(freq, init=0.0, step=0.1) wave = [sin9_shaper(x) for x = init:step:freq] end

julia> @code_native gen_sinwave(1113.0); .section __TEXT,__text,regular,pure_instructions ; Function gen_sinwave { ; Location: REPL[60]:2 pushl %ebx decl %eax subl $48, %esp vmovaps %xmm0, %xmm2 ; Function gen_sinwave; { ; Location: REPL[60]:2 decl %eax movl $769501344, %eax ## imm = 0x2DDDA8A0 addl %eax, (%eax) addb %al, (%eax) decl %eax movl $773805120, %ecx ## imm = 0x2E1F5440 addl %eax, (%eax) addb %al, (%eax) vmovsd (%ecx), %xmm1 ## xmm1 = mem[0],zero decl %eax movl %esp, %ebx vxorps %xmm0, %xmm0, %xmm0 decl %eax movl %ebx, %edi calll %eax decl %eax movl $769544608, %eax ## imm = 0x2DDE51A0 addl %eax, (%eax) addb %al, (%eax) decl %eax movl %ebx, %edi calll %eax ;} decl %eax addl $48, %esp popl %ebx retl nopw %cs:(%eax,%eax) ;}

end # module