Fast Neon 3-Term Cross Product

JesseT · Jul 22, 2010

So, I'm trying to come up with a fast 3D cross product using NEON instructions, and this is currently what I have:

Code:

// cross3(v1, v2) = [ y1*z2 - y2*z1 , z1*x2 - z2*x1 , x1*y2 - x2*y1 ]
float32x4_t cross3(float32x4_t v1, float32x4_t v2) {
    static uint32x4_t const cross3_mask = { 0x00000000, 0x80000000, 0x00000000, 0x00000000 };
    float32x4_t result, temp1, temp2;
    
    // swizzle v1 and v2 into correct terms
    float32x2_t xy1 = vget_low_f32(v1);
    float32x2_t zz1 = vget_high_f32(v1);
    float32x2_t xy2 = vget_low_f32(v2);
    float32x2_t zz2 = vget_high_f32(v2);
    float32x2_t yx1 = vrev64_f32(xy1);
    float32x2_t yx2 = vrev64_f32(xy2);
    zz1 = vdup_lane_f32(zz1, 0);
    zz2 = vdup_lane_f32(zz2, 0);
    
    // result = v1.yxx_ * v2.zzy_
    temp1 = vcombine_f32(yx1, xy1);
    temp2 = vcombine_f32(zz2, yx2);
    result = vmulq_f32(temp1, temp2); 
    
    // result = result - (v1.zzy_ * v2.yxx_)
    temp1 = vcombine_f32(zz1, yx1);
    temp2 = vcombine_f32(yx2, xy2);
    result = vmlsq_f32(result, temp1, temp2); 
    
    // result = result * [ 1, -1, 1 ]
    result = vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(result), cross3_mask));
    return result;
}

And here's the generated assembly output:

Code:

__Z6cross319__simd128_float32_t19__simd128_float32_t:
	@ args = 16, pretend = 0, frame = 0
	@ frame_needed = 0, uses_anonymous_args = 0
	@ link register save eliminated.
	vldmia	sp, {d6-d7}
	vmov	d26, r0, r1  @ v4sf
	vmov	d27, r2, r3
	vmov	d20, d6
	vmov	d2, d26
	vmov	d3, d7
	vrev64.32	d5, d2
	vrev64.32	d1, d20
	vmov	d17, d27
	vdup.32	d6, d3[0]
	ldr	r3, L3
	vmov	d24, d6
	vmov	d25, d1
	vdup.32	d4, d17[0]
	vmov	d22, d5
	vmov	d23, d2
	vmov	d18, d1
	vmov	d19, d20
	vmul.f32	q1, q11, q12
	vmov	d16, d4
	vmov	d17, d5
	vldmia	r3, {d0-d1}
	vmls.f32	q1, q8, q9
	veor	q3, q1, q0
	vmov	r0, r1, d6  @ v4sf
	vmov	r2, r3, d7
	@ lr needed for prologue
	bx	lr

Looks like I'll have to rewrite it using inline assembly to remove the redundant vmov instructions. But other than that, I'm wondering if the general algorithm can at all be improved. There's some ARM marketing material floating around saying that a 3-term cross product can be done in 3 cycles. I'm not sure how that's possible with 2x vrev64 + 2x vdups + vmul + vldmia + vmls + veor, unless the data is already pre-formatted for a vmul + vmls + veor/vmul. Am I missing something obvious here that can make this faster?

JesseT · Jul 23, 2010

Well, I think I've answered my own questions through some hard work. But if anyone is interested, here's what I arrived at.

I rewrote my original version (using vrev/vdup to order terms) in inline assembly to test it, and it was a bit faster than the C version. I also came up a with a better cross product using vzip/vuzp to interleave the elements across multiple vectors, and despite it using up more multiplies, it performs a bit better due to pipelining. However, I then realized I could perform up to 4 cross-products in one shot if I use the interleaved load and store instructions to load in 4 sets of vectors.

Unfortunately, I've discovered the hard way that GCC really sucks at optimization when using NEON intrinsics in some instances. It's often slower than just coding in straight C/C++. I also benchmarked what happens if you've got two vectors in NEON registers, and attempt perform the C-based cross product on them, causing them to be stored into memory, followed by loading the result from memory. It makes a lot of sense to keep things in NEON registers as much as possible.

Here's the benchmarks, run on a 600MHz Cortex A8.

Code:

c 				iterations: 50000000 time: 4.82s
c w/ neon vld/vst1		iterations: 50000000 time: 8.25s
neon intrin vrev/vdup		iterations: 50000000 time: 10.48s
neon asm vrev/vdup		iterations: 50000000 time: 2.79s
neon intrin vzip/vuzp		iterations: 50000000 time: 10.92s
neon asm vzip/vuzp		iterations: 50000000 time: 2.32s
neon asm vld3/vld3 4x		iterations: 50000000 time: 1.53s

And here's the source code for the various routines.

Code:

inline __attribute__((always_inline)) float32x4_t cross3_neon_rev_dup_asm(float32x4_t a, float32x4_t B) {
    register float32x4_t v1 asm("q0") = a;
    register float32x4_t v2 asm("q2") = b;
    register float32x4_t result asm("q9");
    
    asm volatile(
        "vrev64.32 d3, d0             \n\t" // d0,d3 = xy1, yx1
        "vrev64.32 d7, d4             \n\t" // d4,d7 = xy2, yx2
        "vdup.32 d2, d1[0]            \n\t" // q1 = v1.zzy
        "vdup.32 d6, d5[0]            \n\t" // q3 = v2.zzy
        "vmov    d1, d0               \n\t"
        "vmov    d0, d3               \n\t" // q0 = v1.yxx
        "vmov    d5, d4               \n\t"
        "vmov    d4, d7               \n\t" // q2 = v2.yxx
        "vmul.f32 q0, q0, q3          \n\t" // q0 = v1.yxx * v2.zzy
        "vldmia  %[mask], {d16-d17}   \n\t"
        "vmls.f32 q0, q1, q2          \n\t" // q0 = q0 - v1.zzy * v2.xxy
        "veor    q9, q0, q8"                // result = q0 * { 1, -1, 1 }
        : "=w" (v1), "=w" (v2), "=w" (result)
        : "0" (v1), "1" (v2), [mask] "r" (&cross3_mask)
        : "q1", "q3", "q8"
    );

    return result;
}

inline __attribute__((always_inline)) float32x4_t cross3_neon_zip_unzip(float32x4_t v1, float32x4_t v2) {
    float32x4x2_t xxyyzz1 = vzipq_f32(v1, v1);
    float32x4x2_t xxyyzz2 = vzipq_f32(v2, v2);    
    float32x2_t xx1 = vget_low_f32(xxyyzz1.val[0]);
    float32x2_t yy1 = vget_high_f32(xxyyzz1.val[0]);    
    float32x2_t zz1 = vget_low_f32(xxyyzz1.val[1]);
    float32x2_t xx2 = vget_low_f32(xxyyzz2.val[0]);
    float32x2_t yy2 = vget_high_f32(xxyyzz2.val[0]);    
    float32x2_t zz2 = vget_low_f32(xxyyzz2.val[1]);

    float32x2_t x = vmul_f32(yy1, zz2); // x' = v1.y * v2.z
    float32x2_t y = vmul_f32(zz1, xx2); // x' = v1.z * v2.x  
    float32x2_t z = vmul_f32(xx1, yy2); // x' = v1.x * v2.y
    x = vmls_f32(x, zz1, yy2);          // x' = x' - v1.z * v2.y
    y = vmls_f32(y, xx1, zz2);          // y' = y' - v1.x * v2.z
    z = vmls_f32(z, yy1, xx2);          // z' = z' - v1.y * v2.x

    float32x2x2_t result = vuzp_f32(x, y);
    return vcombine_f32(result.val[0], z);
}

inline __attribute__((always_inline)) float32x4_t cross3_neon_zip_unzip_asm(float32x4_t a, float32x4_t B) {
    register float32x4_t v1 asm("q0") = a;
    register float32x4_t v2 asm("q2") = b;
    register float32x4_t result asm("q8");
    
    asm volatile(
        "vmov    q1, q0              \n\t"
        "vmov    q3, q2              \n\t"
        "vzip.32 q0, q1              \n\t" // d0,d1,d3 = xx1, yy1, zz1
        "vzip.32 q2, q3              \n\t" // d4,d5,d6 = xx2, yy2, zz2
        "vmul.f32 d16, d1, d6        \n\t" // x' = yy1 * zz2
        "vmul.f32 d18, d3, d4        \n\t" // y' = zz1 * xx2
        "vmul.f32 d17, d0, d5        \n\t" // z' = xx1 * yy2
        "vmls.f32 d16, d3, d5        \n\t" // x' = x' - zz1 * yy2
        "vmls.f32 d18, d0, d6        \n\t" // y' = y' - xx1 * zz2
        "vmls.f32 d17, d1, d4        \n\t" // z' = z' - yy1 * xx2
        "vuzp.32 d16, d18"                 // result.xyzw = xyzz'
        : "=w" (v1), "=w" (v2), "=w" (result)
        : "0" (v1), "1" (v2)
        : "q1", "q3", "q9"
    );

    return result;
}

inline __attribute__((always_inline)) void cross3_neon_vld3_vst3_asm_4x(float const* v1, float const* v2, float* result) {
    unsigned int const offset = 16; // 128-bit vectors

    asm volatile(
        "vld3.f32 {d0[0], d2[0], d4[0]}, [%[v1]], %[offset]         \n\t"
        "vld3.f32 {d0[1], d2[1], d4[1]}, [%[v1]], %[offset]         \n\t"
        "vld3.f32 {d1[0], d3[0], d5[0]}, [%[v1]], %[offset]         \n\t"
        "vld3.f32 {d1[1], d3[1], d5[1]}, [%[v1]], %[offset]         \n\t"
        "vld3.f32 {d16[0], d18[0], d20[0]}, [%[v2]], %[offset]      \n\t"
        "vld3.f32 {d16[1], d18[1], d20[1]}, [%[v2]], %[offset]      \n\t"
        "vld3.f32 {d17[0], d19[0], d21[0]}, [%[v2]], %[offset]      \n\t"
        "vld3.f32 {d17[1], d19[1], d21[1]}, [%[v2]], %[offset]      \n\t"
        "vmul.f32 q12, q1, q10                                      \n\t"
        "vmul.f32 q13, q2, q8                                       \n\t"
        "vmul.f32 q14, q0, q9                                       \n\t"        
        "vmls.f32 q12, q2, q9                                       \n\t"
        "vmls.f32 q13, q0, q10                                      \n\t"
        "vmls.f32 q14, q1, q8                                       \n\t"        
        "vst3.f32 {d24[0], d26[0], d28[0]}, [%[result]], %[offset]  \n\t"
        "vst3.f32 {d24[1], d26[1], d28[1]}, [%[result]], %[offset]  \n\t"
        "vst3.f32 {d25[0], d27[0], d29[0]}, [%[result]], %[offset]  \n\t"
        "vst3.f32 {d25[1], d27[1], d29[1]}, [%[result]], %[offset]"        
        : [v1] "=r" (v1), [v2] "=r" (v2), [result] "=r" (result)
        : "0" (v1), "1" (v2), "2" (result), [offset] "r" (offset)
        : "memory", "q0", "q1", "q2", "q8", "q9", "q10", "q12", "q13", "q14"
    );
}

JesseT · Jul 23, 2010

I just noticed the cross product implementation in the math-neon library here, and I wanted to run it through the same benchmark as my other routines.

For reference, here's the link: math-neon homepage.

I tried out two different versions of it, the first without inlining as is the case in the math-neon library, and the second version with it force-inlined like my other routines. Call overhead between the two is negligible.

Code:

pure c time			iterations: 50000000 time: 4.79s
c w/ neon vld/vst		iterations: 50000000 time: 8.35s
neon intrin rev/dup		iterations: 50000000 time: 9.24s
neon asm rev/dup		iterations: 50000000 time: 2.70s
neon intrin vzip/vuzp		iterations: 50000000 time: 9.10s
neon asm vzip/vuzp		iterations: 50000000 time: 2.31s
neon asm vld3/vst3 4x		iterations: 50000000 time: 1.53s
* math-neon not inlined		iterations: 50000000 time: 3.94s
* math-neon inlined		iterations: 50000000 time: 3.72s

torpor · Jul 23, 2010

Thanks for the details .. I'm for one very interested in your results. What do you conclude is the best approach, now?

Adventus · Jul 23, 2010

Haya, I wrote math-neon. Yeah not too surprised by the result, my routines are loading/storing the vectors from/to memory every time. When you use intrinsics r0-r3 are used to pass the first argument (the 2nd argument is loaded though) and if you inline the compiler will optimise out the return stall (normally the result is passed back via r0-r3). So if your happy to use neon intrinsics through out your code and inline (this is important, without inlining your single cross products are slower than math-neon), then you'll out perform math-neon's vec operations pretty easily.

JesseT · Jul 25, 2010

torpor said:
Thanks for the details .. I'm for one very interested in your results. What do you conclude is the best approach, now?

Well, if you're building a vector math library for 2D/3D geometry and your primary aim is performance, it appears best to make your vector types just typedefs of the NEON intrinsic vector types found in arm_neon.h (float32x4_t, int32x4_t, etc.), make sure all of your basic functions are inlined, and hopefully the compiler will keep your data in NEON registers as much as possible. This means you can't wrap your vector types in classes or structs, so you won't be able to have any fancy instance methods. You also can't pass your vector types by reference to functions, you need to pass vectors by value. If you pass by reference, or if you use float arrays or class/struct types that wrap the intrinsic types, the compiler will end up copying data to/from memory and the NEON register bank. It doesn't make sense to have to load a vector from memory, perform an operation on it, store it back to memory, only to then load it back in for the next operation.

These ideas are also applicable when implementing a vector math library using SSE, SPU, or AltiVec/VMX128 instrinsics.

I've found that if you aim for a simple API that is modeled after say GLSL/HLSL or OpenCL, you end up with something that is fairly comprehensive and won't require a huge mental context switch as you switch between writing code that runs on a either the CPU or a GPU. GCC's vector types also have overloaded operators built-in, so things like +, -, *, / etc. map to the correct SIMD instructions, and the generated code is pretty decent.

But with GCC, you need to avoid intrinsics that cause unnecessary vmov instructions, so routines with a lot of things like vget_low/high_f32, vcombine_f32, vreinterpret_s32_f32, etc. may need to be handwritten in inline assembly. I'm not sure if Clang/LLVM targeting ARM/NEON is any better at optimizing away these unnecessary vmov instructions when using the intrinsics. You'll probably also want to hand code important processing kernels in assembly as well, like a NEON accelerated vertex skinner for example.

Another thing, with NEON is that you don't need to preprocess your data into SoA order (Structure of Arrays) for batch processing, as can be seen in the 4-way cross product implementation above. This is pretty killer, because it means you don't need to reorder data meant for dynamic vertex buffers to be used with OpenGL.

Adventus said:
Haya, I wrote math-neon. Yeah not too surprised by the result, my routines are loading/storing the vectors from/to memory every time.

I didn't know you were on this forum! Yeah, it's a real shame that the ABI doesn't allow for the first 3 or 4 NEON vector arguments of a function to be passed by value directly in NEON registers.

csigg · Jun 22, 2011

Thanks for posting your results. It was very helpful to get me started implementing my own cross product.
I thought I would post my version. It looks simpler/faster, but I haven't done any profiling.
The 4th component of the result is garbage.

Code:

float32x4_t cross3( const float32x4_t& v0, const float32x4_t& v1 )
{
	float32x4_t result;

	asm volatile( 
		"vext.8 d6, %e2, %f2, 4 \n\t" // y1, z1
		"vext.8 d7, %e1, %f1, 4 \n\t" // y0, z0
		"vmul.f32 %e0, %f1, %e2 \n\t" // z0*x1, ?
		"vmul.f32 %f0, %e1, d6  \n\t" // x0*y1, y0*z1
		"vmls.f32 %e0, %f2, %e1 \n\t" // z0*x1-x0*z1, ?
		"vmls.f32 %f0, %e2, d7  \n\t" // x0*y1-x1*y0, y0*z1-z0*y1
		"vext.8 %e0, %f0, %e0, 4    " // y0*z1-z0*y1, z0*x1-x0*z1
		: "+w" (result.f4) 
		: "w" (v0), "w" (v1)
		: "d6", "d7" );

	return result;
}

Cheers, Christian

Fast Neon 3-Term Cross Product

JesseT

Still Fresh

JesseT

Still Fresh

JesseT

Still Fresh

torpor

hack hack hack, the little machines fight back

Adventus

GP Mania

JesseT

Still Fresh

csigg

Still Fresh