nimcity/doc/part04/ch45.md

Chapter 45. Parallel for each and reduce

There is a form of sharing memory that is particularly simple and effective, and that requires neither locks nor atomic instructions:

1. “Parallel for each”: If a simple operation should be applied to each element of an array and if the iteration order does not matter, the loop can run in parallel.
2. “Parallel reduce”: A sum or product over an array of numbers can be computed by splitting up the array in disjoint slices, then computing the sum or product of every slice and then combining the intermediate results. The slices we work on need to be sendable to the thread pool’s internal task queue, and without causing a copy of the data. Only a (pointer, length) pair should be transmitted. Unfortunately, Nim’s openArray type is not sendable between threads because the compiler cannot guarantee safe access nor safe lifetimes. Instead, we use the type ptr system.UncheckedArray[T] which is exactly what its name suggests: A raw unsafe pointer to an array of unspecified size. There is no index checking. We compute the address of an array element with a custom operator @! in order to make the rest of the code more pleasing to the eye: template @!T: untyped = 1 castptr UncheckedArray[T] 2 1 An operator for array element address computation, also known as pointer arithmetic. 2 The address of the i-th array element is addr data[i] but it needs to be casted into the type ptr UncheckedArray as we will access the successive elements data[i], data[i+1], data[i+2], ... with it. 273 45.1. parMap A “for each” operation is also commonly known as a map. We call our parallel map parMap: import malebolgia template parMapT = 1 proc worker(a: ptr UncheckedArray[T]; until: int) = 2 for i in 0 ..< until: op a[i] 3 var m = createMaster() m.awaitAll: var i = 0 while i+bulkSize <= data.len: 4 m.spawn worker(data@!i, bulkSize) 5 i += bulkSize 6 if i < data.len: 7 m.spawn worker(data@!i, data.len-i) 8 1 parMap takes an openArray, a bulkSize and the operation op to perform. The bulkSize is crucial to make the tasks big enough to amortize the overhead of sending the task to a different CPU. 2 The worker operates on a (pointer, length) slice. 3 The worker applies op to every element of the slice. 4 As long as a slice of bulkSize exists... 5 ...run op on data[i ..< i+bulkSize]. 6 Advance the run index i by bulkSize and proceed with the next slice. 7 The final slice might have fewer elements than bulkSize and needs to be special cased. 8 The final slice has length data.len-i. For parallelization to be worthwhile, the input data array has to be of a sufficient length and the bulkSize must not be too small: 274 var testData: seq[int] = @[] for i in 0 ..< 10_000: testData.add i 1 proc mul4(x: var int) = x *= 4 2 parMap(testData, 600, mul4) 3 for i in 0 ..< 10_000: assert testData[i] == i*4 4 1 Creates test data, an array of 10_000 elements with the values 0, 1, 2, ... 2 mul4 takes a number and multiplies it by 4 and stores the result back into the var parameter. This is passed to the parMap template which mutates the array in place. 3 Calls parMap with the test data, a block size of 600 and the mul4 routine. mul4 could also have been declared as a template and parMap would accept it. 4 Tests that the seq was successfully mutated, the i-th element should have the value i * 4. The power of structured concurrency combined with raw memory accesses and Nim’s template mechanism cannot be underestimated. Of course, these dangerous mechanisms should only be used behind the curtain of a safe abstraction, but that is parMap's purpose. At the same time exposing the bulkSize is crucial for performance tweaking and should not be hidden. Providing a good default value for bulkSize is basically impossible as it depends on the cost of the op snippet that is run on every array element. 275 45.2. parReduce There are few if any languages besides Nim that can express implementations of parMap and parReduce as concisely. parReduce can be implemented like this: template parReduceT: untyped = 1 proc reduceTx: Tx = 2 result = default(Tx) 3 for i in 0 ..< until: op(result, a[i]) 4 var m = createMaster() var res = newSeqint 5 var r = 0 6 m.awaitAll: var i = 0 while i+bulkSize <= data.len: 7 m.spawn reduce(data@!i, bulkSize) -> res[r] 8 r += 1 i += bulkSize if i < data.len: 9 m.spawn reduce(data@!i, data.len-i) -> res[r] 10 r += 1 reduce(res@!0, r) 11 1 parReduce takes an openArray, a bulkSize and which operation op to perform. The bulkSize is crucial to make the tasks big enough to amortize the overhead of sending the task to a different CPU. 2 reduce is a helper proc that takes a slice and accumulates a result. Note that due to a current compiler limitation the inner generic type cannot be named T and so Tx was used. 3 Instead of 0 or 0.0 we use default(Tx) to keep it generic. 4 The actual reduction step. If op is += then op(result, a[i]) is transformed into +=(result, a[i]) which is the same as result += a[i]. 5 For a parallel reduction we need a helper container that keeps the intermediate results. It is named res here. It is very important not to grow the seq after construction! An add might cause a reallocation of the seq which would be disastrous! The binding → res[r] passes the address of res[r] to a worker thread and so it must be stable. 6 r is a helper variable that keeps the currently used length of res. 276 7 As long as a slice of bulkSize exists... 8 ...reduce data[i ..< i+bulkSize] via op. 9 The final slice might have fewer elements than bulkSize and needs to be special cased. 10 The final slice has length data.len-i. 11 All intermediate results need to be reduced after the spawned tasks have been completed. This reduction is also the final result that the template “returns”. Now we need to test our parReduce implementation: var numbers: seq[int] = @[] for i in 0 ..< 10_000: numbers.add i 1 let sum = parReduce(numbers, 600, +=) 2 assert sum == 49995000 3 1 Creates test data, an array of 10_000 elements with the values 0, 1, 2, ... 2 Calls parReduce with the test data, a block size of 600 and the += builtin operator for integers. 3 Ensures the sum has the correct value. Notice how this example works even though += is not a real proc but a builtin thanks to the substitution rules of a template. 277 45.3. parFind Quite analogous to reduce, a search parFind can be written. The task is to return the minimal index of an element that fulfills some criterion or predicate. The helper proc that runs serially over the slice needs to know the offset so that later on the minimum of the results can be taken: template parFindT: int = 1 proc linearFindTx: int = 2 for i in 0 ..< until: if predicate(a[i]): return i + offset 3 return -1 4 var m = createMaster() var res = newSeqint 5 var r = 0 6 m.awaitAll: var i = 0 while i+bulkSize <= data.len: 7 m.spawn linearFind(data@!i, bulkSize, i) -> res[r] 8 r += 1 i += bulkSize if i < data.len: 9 m.spawn linearFind(data@!i, data.len-i, i) -> res[r] 10 r += 1 var result = -1 11 for i in 0 ..< r: if res[i] >= 0: 12 result = res[i] break result 13 1 parFind takes an openArray, a bulkSize and the predicate to search for. It returns the smallest index of an element that fulfills predicate. It produces the value -1 if no such element exists. 2 linearFind is a helper proc that takes a slice and performs a linear search. The offset parameter is used to adjust the index so that it refers to the real position within the openArray and not within the slice. Note that due to a current compiler limitation the inner generic type cannot be named T and so Tx was used. 3 Returns on a successful search. 4 For an unsuccessful search we return -1. 278 5 For a parallel search we need a helper container that keeps the intermediate results. It is named res here. It is very important not to grow the seq after construction! An add might cause a reallocation of the seq which would be disastrous! The binding → res[r] passes the address of res[r] to a worker thread and so it must be stable. 6 r is a helper variable that keeps the currently used length of res. 7 As long as a slice of bulkSize exists... 8 ...search data[i ..< i+bulkSize]. 9 The final slice might have fewer elements than bulkSize and needs to be special cased. 10 The final slice has length data.len-i. 11 Keep in mind that a template does not have an implicitly declared result variable. So we need to declare one here ourselves. 12 Iterates over the intermediate results and stops as soon as a valid index was found. 13 The final value that is produced by the template is result. This time we pass a helper template to parFind for our testing purposes: var haystack: seq[int] = @[] for i in 0 ..< 10_000: haystack.add i 1 template predicate(x): untyped = x == 1000 2 let idx = parFind(haystack, 600, predicate) 3 assert idx == 1000 4 1 Creates test data, an array of 10_000 elements with the values 0, 1, 2, ... 2 predicate takes the current array element and compares it to the number
   1. This is helper template that is then passed to parFind. 3 Calls parFind with the test data, a block size of 600 and the predicate. 4 Ensures that idx has the correct value.