simd/approaches

Source

main.odin ¶

Wrap linesTab width 12345678289 linesSource

1package simd_approaches
2
3import "base:intrinsics"
4import "core:fmt"
5import "core:math/rand"
6import "core:simd"
7import "core:time"
8
9// The number of elements to use for SIMD vectors (in cases where the width is variable).
10WIDTH :: #config(WIDTH, 8)
11
12// The number of objects to use in benchmarking.
13NUM_OBJECTS :: #config(NUM_OBJECTS, 1_000_000)
14
15// Extra padding to add in to the Object. As this increases, *any* AoS solution will suffer due to
16// the decreasing effectiveness of memory caching.
17PADDING :: #config(PADDING, 0)
18
19// Straightforward data layout.
20Object :: struct {
21	pos, vel: [3]f32,
22	_: [PADDING]f32,
23}
24
25// A straightforward implementation using array programming.
26step_aos_scalar :: proc (data: []Object, dt: f32) {
27	for &obj in data {
28		obj.pos += obj.vel * dt
29	}
30}
31
32/*
33We can use SIMD within an object, where a SIMD vector is used similarly to a mathematical vector.
34This can provide moderate speedup without requiring the layout of your data to change
35significantly, but doesn't necessarily scale. Even if your hardware can support wider SIMD
36(e.g. 16 f32s with AVX-512), this approach will only allow you to SIMD up to a single
37(mathematical) vector's worth of values.
38
39However, most any SIMD hardware can handle 128-bit vectors, and as such will benefit from this.
40While it may not technically be the fastest, it can still be a significant speedup.
41
42Using masked loads and stores, you can potentially benefit from SIMD without needing to change
43the data layout at all. Note that on amd64, this approach is comparable to the scalar approach
44with default settings, but becomes significantly more effective with AVX enabled
45(-target-features:avx or -microarch:x86-64-v3). AVX is available on any remotely-recent amd64
46system.
47*/
48
49// Loads a [3]f32 into a #simd[4]f32 using masking. The fourth element is filled with 0. This
50// approach allows for SIMD to be used with minimal alterations to existing data layouts.
51load_vec :: proc (src: ^[3]f32) -> #simd[4]f32 {
52	mask := #simd[4]u32{ 0..<3 = max(u32) }
53	return simd.masked_load(cast(^#simd[4]f32)src, cast(#simd[4]f32)0, mask)
54}
55
56// Stores the first three elements of a #simd[4]f32 into a [3]f32 using masking. This approach
57// allows for SIMD to be used with minimal alterations to existing data layouts.
58store_vec :: proc (dst: ^[3]f32, src: #simd[4]f32) {
59	mask := #simd[4]u32{ 0..<3 = max(u32) }
60	simd.masked_store(cast(^#simd[4]f32)dst, src, mask)
61}
62
63// Using the above load_vec and store_vec, we can SIMD the update without needing to make any
64// changes to the data layout at all. The code itself stays pretty simple too.
65step_aos_within_mask :: proc (data: []Object, dt: f32) {
66	for &obj in data {
67		pos := load_vec(&obj.pos)
68		vel := load_vec(&obj.vel)
69		store_vec(&obj.pos, pos + vel * dt)
70	}
71}
72
73/*
74As another potential improvement, we can actually change the layout of the object so that the
75mathematical vectors are stored as SIMD vectors. This makes loading them easier and possibly
76faster--though this depends on the target.
77
78This also has the downside that it affects the layout of your structs, resulting in a larger
79alignment and a bit of padding after each vector. Additionally, any code which *does* need to
80deal with the individual components of a vector will have a harder time doing so since #simd
81vectors can't be indexed directly.
82*/
83Object_Simd :: struct {
84	pos, vel: #simd[4]f32,
85	_: [PADDING]int,
86}
87
88step_aos_within_simd :: proc (data: []Object_Simd, dt: f32) {
89	for &obj in data {
90		obj.pos += obj.vel * dt
91	}
92}
93
94/*
95What if we want to take advantage of the full width of the SIMD vectors? Modern hardware can have
96as many as 16-float wide vectors (AVX-512), it sure seems restrictive to only be able to use 3.
97`core:simd` has gather and scatter intrinsics that we can use to load and store each element of
98a vector from arbitrary locations, provided as pointers. In theory, this allows us to take
99advantage of the full width of the SIMD vector without needing to rearrange the struct at all!
100
101There's just one problem: this approach is often very slow--possibly even significantly slower
102than the scalar version, depending on your hardware.
103
104On amd64, the gather instruction is only available with AVX2, so it can't be used with the
105default compilation settings. Even then LLVM will tend to avoid using the actual gather
106instruction in the x86-64-v3 microarch, and with good reason--it's quite slow on many CPUs. This
107isn't just on old CPUs, either--even on a recent Ryzen 9950X, the alternative code that LLVM
108generates to load the values in scalar fashion and load them into a vector is *still* faster than
109the version it generates that uses the actual gather instruction. This may vary depending on the
110hardware, so test on your target hardware if you know what it will be (and avoid gather if you
111don't). Scatter is supported by even less hardware.
112
113For amd64, to enable the use of hardware gather instructions, use
114-target-features:avx2,fast-gather . Hardware scatter requires AVX-512
115( -target-features:avx512f,avx512vl ), but that's rare and may not fare much better even on
116hardware that has it.
117*/
118step_aos_gather :: proc (data: []Object, dt: f32) {
119	data := data
120
121	do_add :: #force_inline proc (base: ^Object, dt: f32, mask: #simd[WIDTH]u32 = max(u32)) {
122		index := iota(#simd[WIDTH]uintptr)
123		step := index * size_of(Object)
124
125		// Generate a pointer to each value of interest
126		px_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.pos.x) + step)
127		vx_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.vel.x) + step)
128		px := simd.gather(px_ptr, cast(#simd[WIDTH]f32)0, mask)
129		px += simd.gather(vx_ptr, cast(#simd[WIDTH]f32)0, mask) * dt
130		simd.scatter(px_ptr, px, mask)
131
132		py_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.pos.y) + step)
133		vy_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.vel.y) + step)
134		py := simd.gather(py_ptr, cast(#simd[WIDTH]f32)0, mask)
135		py += simd.gather(vy_ptr, cast(#simd[WIDTH]f32)0, mask) * dt
136		simd.scatter(py_ptr, py, mask)
137
138		pz_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.pos.z) + step)
139		vz_ptr := cast(#simd[WIDTH]rawptr)(uintptr(&base.vel.z) + step)
140		pz := simd.gather(pz_ptr, cast(#simd[WIDTH]f32)0, mask)
141		pz += simd.gather(vz_ptr, cast(#simd[WIDTH]f32)0, mask) * dt
142		simd.scatter(pz_ptr, pz, mask)
143	}
144
145	for len(data) >= WIDTH {
146		do_add(raw_data(data), dt)
147		data = data[WIDTH:]
148	}
149
150	if len(data) > 0 {
151		mask := simd.lanes_lt(iota(#simd[WIDTH]u32), cast(#simd[WIDTH]u32)len(data))
152		do_add(raw_data(data), dt, mask)
153	}
154}
155
156/*
157Finally, if we want to go all-out and take full advantage of the hardware, we need to rearrange
158our data layout. Rather than storing all of the data for each object together, in
159Array-of-Structs form, we separate them so that each field's data is stored in a separate array.
160We call this Struct-of-Arrays (SoA) form, as its in memory layout more closely resembles a struct
161where each field is an array, rather than an array where each value is a complete struct.
162
163Odin's #soa tag helps significantly with this, allowing you to write code that accesses SoA data
164but still looks mostly like AoS data. We do, unfortunately, have to sacrifice the fixed arrays
165for pos/vel as otherwise the compiler will still group their components together.
166
167For this particular update procedure, this data layout can be extremely fast. All of the data it
168uses is stored consecutively with no gaps, so the SIMD loads and stores can operate directly on
169it and also take advantage of the full width of the SIMD vector. Try playing with WIDTH in
170conjunction with AVX and/or AVX512 (if you have it)!
171
172Because each field's data is stored separately, it's also unaffected by other data that may be
173stored in the struct that isn't being used here. Try increasing PADDING--the other approaches
174will tend to slow down as the data they operate on becomes spaced further apart, but this one
175will not since the position and velocity will always be tightly-packed. However, for random
176access of the data in the array it can end up being slower, due to the possibility of being
177multiple cache misses instead of just one (not shown in these benchmarks).
178*/
179Object_Split :: struct {
180	px, py, pz: f32,
181	vx, vy, vz: f32,
182	_: [PADDING]int,
183}
184
185step_soa :: proc (data: #soa[]Object_Split, dt: f32) {
186	data := data
187
188	do_add :: #force_inline proc (base: #soa[]Object_Split, first: int, dt: f32, mask: #simd[WIDTH]u32 = max(u32)) {
189		px := intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.px[first:])
190		px += intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.vx[first:]) * dt
191		intrinsics.unaligned_store(cast(^#simd[WIDTH]f32)base.px[first:], px)
192
193		py := intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.py[first:])
194		py += intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.vy[first:]) * dt
195		intrinsics.unaligned_store(cast(^#simd[WIDTH]f32)base.py[first:], py)
196
197		pz := intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.pz[first:])
198		pz += intrinsics.unaligned_load(cast(^#simd[WIDTH]f32)base.vz[first:]) * dt
199		intrinsics.unaligned_store(cast(^#simd[WIDTH]f32)base.pz[first:], pz)
200	}
201
202	i: int
203	for i = 0; i+WIDTH <= len(data); i += WIDTH {
204		do_add(data, i, dt)
205	}
206
207	if i < len(data) {
208		left := len(data) - i
209		mask := simd.lanes_lt(iota(#simd[WIDTH]u32), cast(#simd[WIDTH]u32)left)
210		do_add(data, i, dt, mask)
211	}
212}
213
214main :: proc () {
215	if ODIN_OPTIMIZATION_MODE <= .Minimal {
216		fmt.println("WARNING: For best results, run benchmarks in an optimized build!")
217	}
218
219	aos_data := make([]Object, 1_000_000)
220	for &o in aos_data {
221		o = make_object()
222	}
223
224	aos_simd_data := make([]Object_Simd, NUM_OBJECTS)
225	for &o in aos_simd_data {
226		temp := make_object()
227		o = {
228			pos = {temp.pos.x, temp.pos.y, temp.pos.z, 0},
229			vel = {temp.vel.x, temp.vel.y, temp.vel.z, 0},
230		}
231	}
232
233	soa_data := make(#soa[]Object_Split, NUM_OBJECTS)
234	for &o in soa_data {
235		temp := make_object()
236		o = {
237			px = temp.pos.x, py = temp.pos.y, pz = temp.pos.z,
238			vx = temp.vel.x, vy = temp.vel.y, vz = temp.vel.z,
239		}
240	}
241
242	bench_2("AoS Scalar", step_aos_scalar, aos_data, 1.0/60.0)
243	bench_2("AoS Within (Masking)", step_aos_within_mask, aos_data, 1.0/60.0)
244	bench_2("AoS Within (Vector Storage)", step_aos_within_simd, aos_simd_data, 1.0/60.0)
245	bench_2("AoS Across (Gather)", step_aos_gather, aos_data, 1.0/60.0)
246	bench_2("SoA Across", step_soa, soa_data, 1.0/60.0)
247}
248
249make_object :: proc () -> Object {
250	return {
251		pos = {
252			rand.float32_range(-1000, 1000),
253			rand.float32_range(-1000, 1000),
254			rand.float32_range(-1000, 1000),
255		},
256
257		vel = {
258			rand.float32_range(-10, 10),
259			rand.float32_range(-10, 10),
260			rand.float32_range(-10, 10),
261		},
262	}
263}
264
265bench_2 :: proc (name: string, p: proc($A, $B), a: A, b: B) {
266	fmt.print(name, "... ", sep="")
267
268	iterations, done := 1, 0
269	start := time.tick_now()
270	for time.tick_since(start) < time.Second {
271		for _ in 0..<iterations {
272			p(a, b)
273		}
274
275		done += iterations
276		iterations += iterations
277	}
278	elapsed := time.tick_since(start)
279	ips := f64(done) / f64(time.duration_seconds(elapsed))
280	fmt.println(ips, "per sec")
281}
282
283iota :: proc ($T: typeid/#simd[$N]$E) -> (result: T) {
284	for i in 0..<N {
285		result = simd.replace(result, i, E(i))
286	}
287	return
288}
289