simd/basic-sum
Code
main.odin ¶213 linesSource
1package main
2
3import "base:intrinsics"
4import "core:fmt"
5import "core:math/rand"
6import "core:simd"
7import "core:time"
8
9MISALIGN :: #config(MISALIGN, false) // Whether to misalign the data, for testing purposes
10NUM_DATA :: #config(NUM, 10_000_000) // The amount of data points to generate and process
11NUM_REPETITIONS :: #config(REP, 100) // The number of times to run each proc, for performance measurement
12
13// Calculates the sum of an array of f32, the basic way.
14// Due to LLVM math settings that Odin doesn't provide a way to set, this won't (currently) be
15// autovectorized.
16sum_scalar :: proc (s: []f32) -> (sum: f32) {
17 for x in s {
18 sum += x
19 }
20 return
21}
22
23// Calculates the sum of an array of f32, the basic way, using f64s to hold the sum.
24// This gives a more precise result, as floating-point precision results in a cumulative error
25// when adding many numbers together. The performance cost is negligible compared to summing in f32
26// (on amd64).
27sum_scalar_wide :: proc (s: []f32) -> f32 {
28 sum : f64
29 for x in s {
30 sum += f64(x)
31 }
32 return f32(sum)
33}
34
35// The number of elements to use in SIMD vectors.
36// The best value for this will depend on the available SIMD instructions, and possibly the hardware
37// itself.
38//
39// On amd64, the default target only uses SSE4, which has 128-bit SIMD registers. This means that
40// #simd[4]f32 would be the native vector size on that target--but that doesn't always give the
41// fastest results, as larger vectors allow for better instruction-level parallelism. For larger
42// vectors, LLVM will automatically spread the data over multiple SIMD registers, but if the SIMD
43// vector is too large this starts to become a detriment to performance. Try comparing the results
44// between 4, 8, and 16! This particular example isn't complex enough to suffer from too large of
45// vectors, even with a WIDTH of 64.
46//
47// You will likely see a performance boost, particularly in the f64 case, by enabling AVX, and
48// around 97% of PCs support that (according to the October 2024 Steam hardware survey). You can
49// enable that by building with -target-features:avx (or -microarch:x86-64-v3, for a number of CPU
50// features available on many modern systems). However, note that doing so will cause the program to
51// crash on systems that don't support these features!
52WIDTH :: #config(WIDTH, 16)
53
54// Calculates the sum of an array of f32 using SIMD.
55// Like with the scalar version, this is susceptible to limitations of floating-point precision--
56// however, you'll notice that the result is often different! This is why LLVM won't auto-vectorize
57// that version. By default, LLVM handles floating-point math in a strict fashion, where it will
58// only perform optimizations that don't change the order of the math operations (as that can change
59// the result).
60sum_simd :: proc (s: []f32) -> f32 {
61 s := s
62
63 vec_sum : #simd[WIDTH]f32
64 for len(s) >= WIDTH {
65 chunk_ptr := cast(^#simd[WIDTH]f32)raw_data(s)
66 s = s[WIDTH:]
67
68 // SIMD vectors are sensitive to alignment, so an unaligned load is used here.
69 // A f32 has an alignment of 4 bytes, whereas SIMD registers often have an alignment of
70 // 16/32/64 bytes. The input slice may not meet this alignment, so an unaligned load is
71 // used. A regular dereference will load the vector with strict alignment, and if the
72 // pointer isn't correctly aligned that will crash the program!
73 //
74 // There are other ways of dealing with this too (e.g. process in scalar
75 // until you reach the appropriate alignment), but using an unaligned load is
76 // simple and has a negligible cost on most modern systems.
77 chunk := intrinsics.unaligned_load(chunk_ptr)
78
79 // While it may be intuitive to reduce the intermediate vector to a single value here, that
80 // will actually reduce the parallelism and hurt performance! Generally, for best SIMD
81 // performance, stay wide as long as possible.
82 vec_sum += chunk
83 }
84
85 // Reduces all of the elements in the vector to a single value by adding them (also known as a
86 // horizontal add). In this case, this gives the sum of all the values processed so far.
87 sum := simd.reduce_add_ordered(vec_sum)
88
89 // Since the vectorized part above worked in chunks of size WIDTH, any leftover needs to be
90 // handled separately. There are multiple ways to do this; in this case, we just process the
91 // remaining few values in scalar fashion.
92 //
93 // Note that for more complex logic, this could result in a duplication of logic! See
94 // sum_simd_masked for an alternative.
95 for x in s {
96 sum += x
97 }
98 return sum
99}
100
101// Calculates the sum of an array of f32 using SIMD, using f64s to hold the sum.
102// This gives a more precise result, as floating-point precision results in a cumulative error
103// when adding many numbers together. Unlike with the scalar procs, though, the performance cost vs.
104// the pure f32 version is much more significant (but still much faster than the non-SIMD variants).
105// This is because half as many f64s can fit in a single SIMD register, reducing parallelism. It's
106// still significantly faster than the scalar versions, though!
107sum_simd_wide :: proc (s: []f32) -> f32 {
108 s := s
109
110 vec_sum : #simd[WIDTH]f64
111 for len(s) >= WIDTH {
112 chunk_ptr := cast(^#simd[WIDTH]f32)raw_data(s)
113 s = s[WIDTH:]
114
115 chunk := cast(#simd[WIDTH]f64)intrinsics.unaligned_load(chunk_ptr)
116 vec_sum += chunk
117 }
118
119 sum := simd.reduce_add_ordered(vec_sum)
120 for x in s {
121 sum += f64(x)
122 }
123 return f32(sum)
124}
125
126// Calculates the sum of an array of f32 using SIMD with masking to handle non-multiples-of-WIDTH.
127sum_simd_masked :: proc (s: []f32) -> f32 {
128 s := s
129
130 process_chunk :: proc (chunk_ptr: ^#simd[WIDTH]f32, mask: #simd[WIDTH]u32, sum: #simd[WIDTH]f32) -> #simd[WIDTH]f32 {
131 // Masked loads only load the elements that are selected in the mask. Memory locations that
132 // are not selected in the mask are effectively not touched, so it's safe for some elements
133 // to be "past the end" of the available data, as long as they aren't selected by the mask.
134 // Values that aren't selected in the mask receive their values from the corresponding
135 // element in the second parameter instead (in this case, 0).
136 //
137 // Masked loads and stores don't have the strict alignment requirements that raw dereferences do.
138 chunk := simd.masked_load(chunk_ptr, cast(#simd[WIDTH]f32)0, mask)
139
140 // In this case the operation is trivial--but for more complex operations, sharing the code
141 // can be beneficial. See the motion example for a more complex example.
142 return sum + chunk
143 }
144
145 vec_sum : #simd[WIDTH]f32
146
147 for len(s) >= WIDTH {
148 chunk_ptr := cast(^#simd[WIDTH]f32)raw_data(s)
149 s = s[WIDTH:]
150
151 mask := cast(#simd[WIDTH]u32)max(u32) // Selects every element
152 vec_sum = process_chunk(chunk_ptr, mask, vec_sum)
153 }
154
155 // Handle any leftovers via masking. This can be combined with the above loop, but that will result in worse performance
156 if len(s) > 0 {
157 chunk_ptr := cast(^#simd[WIDTH]f32)raw_data(s)
158
159 // This mask will select vector elements where the element's value in index is less than the
160 // remaining length of the slice (in this case, the elements that are within the bounds of
161 // the slice)
162 // Comparisons generate a vector with integer elements, where each element is either 0
163 // (false) or non-zero (true), depending on the result of the comparison
164 index := iota(#simd[WIDTH]i32)
165 mask := simd.lanes_lt(index, cast(#simd[WIDTH]i32)len(s))
166
167 vec_sum = process_chunk(chunk_ptr, mask, vec_sum)
168 }
169
170 return simd.reduce_add_ordered(vec_sum)
171}
172
173main :: proc() {
174 if ODIN_OPTIMIZATION_MODE <= .Minimal {
175 fmt.println("WARNING: For best results, run benchmarks in an optimized build!")
176 }
177
178 when MISALIGN {
179 data_original := make([]f32, NUM_DATA + 1, context.temp_allocator)
180 data := ([^]f32)(uintptr(raw_data(data_original)) + 1)[:NUM_DATA]
181 } else {
182 data := make([]f32, NUM_DATA, context.temp_allocator)
183 }
184
185 for &x in data {
186 x = rand.float32()
187 }
188
189 fmt.printfln("Sum (Scalar f32): %0.5f (%v)", benchmark(sum_scalar, data))
190 fmt.printfln("Sum (SIMD f32): %0.5f (%v)", benchmark(sum_simd, data))
191 fmt.printfln("Sum (Masked SIMD f32): %0.5f (%v)", benchmark(sum_simd_masked, data))
192 fmt.printfln("Sum (Scalar f64): %0.5f (%v)", benchmark(sum_scalar_wide, data))
193 fmt.printfln("Sum (SIMD f64): %0.5f (%v)", benchmark(sum_simd_wide, data))
194}
195
196benchmark :: proc (p: proc ([]f32) -> f32, s: []f32) -> (f32, time.Duration) {
197 best_elapsed := max(time.Duration)
198 result : f32
199 for _ in 0..<NUM_REPETITIONS {
200 start := time.tick_now()
201 result = p(s)
202 best_elapsed = min(time.tick_since(start), best_elapsed)
203 }
204 return result, best_elapsed
205}
206
207iota :: proc ($V: typeid/#simd[$N]$E) -> (result: V) {
208 for i in 0..<N {
209 result = simd.replace(result, i, E(i))
210 }
211 return
212}
213