/
threading-the-sieve-in-python.html
246 lines (198 loc) · 25.1 KB
/
threading-the-sieve-in-python.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
<!---
;;
;; This softeware is Copyright (c) 2009 A.F. Haffmans
;;
;; This file is part of cl-bliky.
;;
;; cl-bliky is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; cl-bliky is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with cl-bliky. If not, see <http://www.gnu.org/licenses/>.
;;
;;
--->
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
<html>
<head>
<title> Mohegan SkunkWorks </title>
<link rel="stylesheet" href="/style.css" type="text/css"/>
</head>
<body>
<div id="blog-title">
<a href="/"> <h1> Mohegan SkunkWorks </h1> </a>
</div>
<div id="mainClm">
<div id="post">
<div id="post-type" style="visibility:hidden">
<p>"post"</p>
</div>
<div id="post-timestamp" >
<p>Sat, 12 Sep 2009 17:03:23 EST </p>
</div>
<div id="post-header">
<h2> Threading the Sieve in Python </h2>
</div>
<div id="post-intro">
<P>This is the first of two posts on threading and multiprocessing in Python. In this post I'll explore the thread module and in the second post I'll look at Python's multiprocessing module. My starting point is the multi-threaded implementation of the <A HREF="http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes">Sieve of Erasthones</A> found in this <A HREF="http://heather.cs.ucdavis.edu/~matloff/Python/PyThreads.pdf">tutorial on multi-threading in Python (pdf).</A><BR> <BR> Threading a compute-bound algorithm, like the <EM>Sieve</EM> consists of subdividing of the main task into autonomous sub-tasks which share as little state as possible. Having no shared state eliminates the overhead that inevitably comes with locking. It turns out that Python is not very good at multi-threading compute-bound processes. <A HREF="http://www.dabeaz.com/blog/dablog.html">This </A><A HREF="http://ttimo.vox.com/library/post/python-gil-threading-and-multicore-hardware.html">is </A><A HREF="http://www.grouplens.org/node/244">not a </A><A HREF="http://blog.ianbicking.org/gil-of-doom.html">surprise.</A> CPython has a global interpretor lock <A HREF="http://www.dabeaz.com/python/GIL.pdf">(GIL)</A> which prevents threads from running concurrently. <BR> <BR> Regardless, there are other lessons I learned when multi-threading the <EM>Sieve</EM> algorithm. One is that sharing state between threads may be unavoidable to achieve reasonable performance. In fact, if you <EM>don't</EM> share state, performance can become predictable <EM>worse</EM> as the number of threads of execution increases. <BR> <BR> The other is that locking can have a surprising impact on performance. It's not just the cost of locking per se, but the effect locking has on the distribution of work between the various threads. </P>
</div>
<div id="post-body">
<H2>Sieve of Erasthones</H2><P>The Sieve of Erasthones is a way to find all the prime numbers smaller than N. You start with an array of size N, with all slots initialized to 1 (i.e. to 'true'). <BR>
You start moving down the array, and if the value of the slot at your current position is 'true' (i.e. 1), you set the slots at multiples of your current position index to false (i.e. 0). At each position there only one transition, from true (1) to false (0), and not vice-versa. The largest multiplier for a position with index i is N/i. You're done when you hit the slot with index equal to the square root of N. Note, that as you make your way through the array, you zero out less and less positions. </P><P>For example, if you want to find all the primes up to 10, you start at position index 2, and zero out 4,6,8 and so on. You move to 3, and zero out 6 and9. since N = 10, you stop here. The next non-zero position is 5, as 4 was already set to false previously. All slots still flagged as true (1) are primes : 2, 3, 5 and 7. </P><H2>Testing Platform</H2><P>All tests are performed on a Toshiba A215-S4747 with an AMD Turion 64 X2 (dual-core) processor. The operating system is Ubuntu 8.04.3. I replaced the Python version shipped with Ubuntu with version 2.6.2. </P><P>For the Jython test I installed Jython 2.5, from the Jython web site in stand-alone mode. IronPython 1.1.1. was installed as well. </P><P>Each implementation I'm going to discuss below is part of <EM>prime_share.py</EM> in the <A HREF="git://github.com/fons/blog-code.git">blog-code package</A> on github. </P><P><EM>prime_share.py</EM> is designed to perform side-by-side comparisons of the various implementations. It allows you to specify a range of threads used in the distribution of the calculation. Timing is done using Python's timeit module, with garbage collection turned off. </P><H2>Various Implementations</H2><P>This section discusses four different multi-threaded implementations of the <EM>Sieve</EM>. </P><UL><LI><P>main_orig: From the <A HREF="http://heather.cs.ucdavis.edu/~matloff/Python/PyThreads.pdf">tutorial (pdf).</A></P></LI><LI><P>main_nolocks: Drastically reduces locking.</P></LI><LI><P>main_nolocks_alt: Splits the sieve across threads.</P></LI><LI><P>main_nolocks_alt2: Distributes the work more equitably amongst the threads.</P></LI></UL><P>The main difference between the implementations is the amount of locking, and the way work is distributed across multiple threads. </P><P>I'll start out <A HREF="#allresults"> by showing </A> the result of a performance test of each of these algorithms. The x-axis shows the number of threads used in the run. Each run calculates the number of primes up to 10000. The time it took to run this calculation 40 times in sucession is shown on the y-axis. Each data set is labelled 'main_xyz' where 'xyz' identifies the implementation. </P><P><A NAME="allresults"><IMG SRC="http://github.com/fons/blog-images/raw/master/thread.png"></A></P><P>You would expect the run time to decrease as the number of threads increases. That's clearly not the case. The performance the <EM>main_orig</EM> implementation is extremenly erratic, and the performance of the <EM>main_nolock_alt</EM> implementation <EM>decreases</EM> steadily as the number of threads increases.<BR>
The performance of the remaining two implementations is unaffected by the number of threads. </P><P>What's the reason for this ? Well, below are four sections discussing these four implementations. </P><H3>main_orig : the original implementation</H3><P>As I mentioned before, the implementation of the <EM>Sieve</EM> as shown in the <A HREF="http://heather.cs.ucdavis.edu/~matloff/Python/PyThreads.pdf">tutorial</A> is the starting point for subsequent implementations. The full code I used in my test is shown below : </P><PRE><CODE>def count_primes(prime) :
p = reduce(lambda x, y: x + y, prime) - 2
return p
def dowork_orig(tn): # thread number tn
global n,prime_global,nexti_global,nextilock,nstarted,nstartedlock,donelock
donelock[tn].acquire()
nstartedlock.acquire()
nstarted += 1
nstartedlock.release()
lim = math.sqrt(n)
while 1:
nextilock.acquire()
k = nexti_global
nexti_global += 1
nextilock.release()
if k > lim: break
if prime_global[k]:
r = n / k
for i in range(2,r+1):
prime_global[i*k] = 0
donelock[tn].release()
def main_orig(top, nthreads):
global n,prime_global,nexti_global,nextilock,nstarted,nstartedlock,donelock
n = int(top)
nthreads = int(nthreads)
prime_global = (n+1) * [1]
nstarted = 0
nexti_global = 2
nextilock = thread.allocate_lock()
nstartedlock = thread.allocate_lock()
donelock = []
for i in range(nthreads):
d = thread.allocate_lock()
donelock.append(d)
thread.start_new_thread(dowork_orig,(i,))
while nstarted < nthreads: pass
for i in range(nthreads):
donelock[i].acquire()
return count_primes(prime_global) </CODE></PRE><P>The variable <EM>n</EM> is the upper limit of the search and <EM>nthreads</EM> is the number of threads. The global variable <EM>prime_global</EM> is the sieve, and <EM>nexti_global</EM> is its index. </P><P>The function <EM>dowork_orig</EM> implements the sieve algorithm along the lines mentioned earlier. </P><P>All the global variables are shared amongst the threads. There are three locking variables.<BR>
The first one is <EM>nstartedlock</EM>, which is used to set a counter. The counter is used in the main to wait for all threads to start. </P><P>The second one is the <EM>donelock</EM> array. This is used to implement a mechanism to wait for all the threads to finish, similar to 'join' in other threading packages. The number of entries is equal to the number of threads. </P><P>The last lock is <EM>nexti_lock</EM>. It's purpose is to protect access to the global index variable so that each thread processes a unique index.<BR>
</P><P>Access to the 'sieve' variable <EM>prime_global</EM> is not protected by lock. There is no need to, since the value of each slot can only go from 1 (true) to 0 (false). If two threads access the same slot simultaneously they can't reach conflicting conclusions of it's final state. Obviously, the GIL effectively prevents any of this from happening, but it's good the know that if it didn't the algorithm wouldn't break. Even out-of-order execution where a higher index is processed first, is not going to lead to different final state (provided all indexes are processed). </P><P>As you can see in the <A HREF="#allresults"> side-by-side comparison </A> the performance of this implementation is very erratic, and shows no obvious dependence on the number threads. </P><P><IMG SRC="http://github.com/fons/blog-images/raw/master/thread_repeat.png"></P><P> Here I show five data sets generated by running the <EM>main_orig</EM> implementation with the same input parameters as in the <A HREF="#allresults"> side-by-side comparison </A> above. Note the how the speed varies dramatically in each data set, and between data sets. </P><P>For a single thread the performance in of all four data sets is roughly the same. The small differences are probably due to other activity on the box taking some time away from the test run. In general the calculation time increases as the number of threads increases, except for the first ('red') data set. </P><P>The reason for the wide variation in performance I believe are the <EM>nstartedlock</EM> and <EM>nextilock</EM>. </P><P>Each thread tries to acquire <EM>nstartedlock</EM> at startup. My suspicion is that this lock acquisition varies<BR>
the start of the threads sufficiently to change to way work is balanced between the threads for each subsequent run, leading to the dramatic variation in speed. That's because the amount of work a thread performs depends on the index it's processing. For example, let's say that we have run of three threads. The first thread is active, and the other two threads are dormant. The first thread will process indices two and three, and therefore do most of the work. When the other threads wake up, they 'll do less work. On a second run, a context switch occurs sooner and the second thread works on index three, leading to a more equitable distribution of the work. </P><P>In addition, there is contention on <EM>nextilock</EM>. As I pointed out, for higher index values, less work is done. So when the prime check return false (as it would most of the time for these values), the thread will try to acquire the lock again which puts it in contention with other threads. The way the lock is hit by a particular thread probably varies bit from run to run, and this results in the erratic performance profile. </P><P>The main difference between this implementation and the three others discussed below is the removal of this lock, and the lock around the index variable. As you <A HREF="#allresults"> can see </A> this has dramatic impact on performance. </P><H3>main_nolocks: remove locking; load balance naively</H3><P>In this implementation I've removed the <EM>nstartedlock</EM> as well as the <EM>nextilock</EM> lock. The <EM>donelock</EM> array is kept in place, since the thread package has no 'join' mechanism. </P><P>Here's the complete code for the variation. The other two main_nolocks* which I' m discussing below are very similar. </P><PRE><CODE>def dowork2(n, nexti_ns, prime_nl) :
k = nexti_ns[0]
lim = nexti_ns[1]
if nexti_ns[0] > nexti_ns[1] :
raise "boundaries out-of-order"
while 1 :
if not (k < lim) : break
if prime_nl[k] == 1 :
r = n / k
for i in range(2, r+1) :
prime_nl[i*k] = 0
k = k + 1
return prime_nl
def dowork_th(tn, donelock, n, nexti_ns) :
global prime_nl
prime_nl = dowork2(n, nexti_ns, prime_nl)
donelock.release()
def load_balance(s, th) :
len_s = len(s)
if len_s == 0 :
return [ s ]
base = len_s / th
rem = len_s - base * th
K = map(lambda i : i * base, range(1, th+1))
t = range(1, rem + 1 ) + (th - rem )*[rem]
K = map(lambda p : p[0] + p[1], zip(K, t))
K = zip([0] + K, K)
last = s[len_s - 1]
s.append(last+1)
K = map(lambda p : (s[p[0]], s[p[1]]), K)
return K
def start_th(fn, args) :
return thread.start_new_thread(fn, args)
def main_nolocks(top, nthreads) :
global prime_nl
n = int(top)
nthreads = int(nthreads)
prime_nl = (n + 1) * [1]
donelock = map(lambda l : l.acquire() and l,
map(lambda i : thread.allocate_lock(), range(nthreads)))
lim = int(math.sqrt(n)) + 1
nexti_ns = range(2, lim, 1)
B = load_balance(nexti_ns, nthreads)
map(lambda i : start_th(dowork_th, (i, donelock[i], n, B[i])),
range(nthreads))
map(lambda i : donelock[i].acquire(), range(nthreads) )
return count_primes(prime_nl) </CODE></PRE><P>In the original implementation the <EM>nexti_global</EM> variable was shared between the threads and used to balance the load between them. In this version that's replaced by the <EM>load_balance</EM> routine. </P><P><EM>load_balance</EM> takes as its input an array and the number of threads and returns a list of pairs which represent the start and end point of the sub arrays each thread will work on. </P><P>For example, if we want to know the number of primes smaller than 101, the indexes we consider in the sieve algorithm would run from 2 to 10. If the number of threads is 3, <EM>load_balance</EM> would return the following array : <CODE>[ (2, 5), (5, 8), (8, 11) ]</CODE> </P><P>This is naive, because the thread working on the first pair will do a lot more work than the one assigned the last pair. </P><P>The difference in performance with the original implementation is <A HREF="#allresults"> quite dramatic. </A></P><P>The big difference in performance even for a single thread between this and the previous version must be due to the difference in locking. </P><P>The original implementation acquires a lock at the startup of each thread and uses locks to control access to the global index counter. In this implementation these locks have been removed. The benefit this removal is quite obvious from the graph. </P><P>The effect of the GIL is also quite obvious: The performance does not depend on the number of threads at all. </P><H3>main_nolocks_alt: remove locking, split the sieve</H3><P>This version is a minor variation of the previous one. Previously the <EM>range</EM> of indices from 2 to sqrt(N) was split up between the threads. </P><P>In this version each thread processes the full index range, but the <EM>sieve array</EM> is split between the threads: One thread works on the first part of the sieve, the second on a second part and so on. </P><P>For example, if the size of the sieve is 100, and we have three threads, the first thread works on the sieve up to index 33, the second thread takes indexes 34 to 66 and thread 3 takes the rest. </P><P>The <EM>dowork</EM> routine changes a little bit : </P><PRE><CODE>def dowork3(n, irange, prime_nla) :
k = 2
lim = int(math.sqrt(n)) + 1
istart, iend = irange
while 1 :
if not (k < lim) : break
if not (k < iend) : break
if k < istart :
s = (istart / k ) + 1
r = (iend / k) + 1
for i in range(s, r) :
prime_nla[i*k] = 0
elif prime_nla[k] == 1 :
assert k >= istart and k <= iend
s = 2
r = (iend / k) + 1
for i in range(s, r) :
prime_nla[i*k] = 0
k = k + 1
return prime_nla </CODE></PRE><P>If the index is smaller than the starting index of the part of the sieve array under consideration, we can just zero out multiples of the index. If that's not case we proceed as before, but stop at the upper index of the array if it's smaller than sqrt(N). </P><P>If you look <A HREF="#allresults"> at the results </A> , notice that the single-threaded performance of this and the previous version are the same. That's to be expected as both versions perform the same amount of work. Less expected is that the performance of this version decreases steadily as the number of threads increases. </P><P>The reason is that the amount of work performed across all threads (i.e. on the process level) increases as the number of threads increases. That's because threads working on the 'higher' part of the sieve can't take advantage of the primes 'discovered' by the thread working on the first, 'lower' part of the sieve. </P><P>Suppose we start with N = 100, and two threads. The first thread works on the sieve range from 0 to 50. It performs the basic sieve operations for the case N = 50. The second thread basically zeros out multiples of values in the range 2 to 11. This thread can't probe the lower part of the sieve to see if it's working with a prime, so it needs to calculate multiples of every element in the range. The more disconnected pieces, the more extra work is done. </P><P>In the absence of the GIL, the additional work per thread could be balanced by the performance gain from distributing the load across multiple threads. I"ll revisit this when tmultiprocessing module is discussed, as it doesn't use a GIL. </P><H3>main_nolocks_alt2: remove locking, load balance more accurately</H3><P>In this last version I try to distribute the work more evenly amongst the threads. This is in contrast with to just splitting the range in pieces, and assigning each thread a piece. This leads to the thread assigned the lowest index range doing most of the work. </P><P>To accomplish a more equitable distribution, I estimate the number of operations per index i. I then assign each thread a set of indices so that the total number of operations per thread varies very little. </P><P>An upper limit for the number of operations for index i for a sieve of length N is : </P><PRE><CODE>1 + N/i - i </CODE></PRE><P>This assumes that the sieve array is not shared amongst the threads. If it is, than we're over-counting the number of operations, since it also counts operations on indices that are multiples of each other. </P><P>For example, when N = 200, for i = 2 the number of operations is 99. For i = 4 the that number is in fact 0 as all multiples of four are also multiples of two. </P><P>Here's the code fragment where the load balancing takes place : </P><PRE><CODE>def smp_load_balance(th , n) :
def operations(t) :
return int((n / t) + 1 - t)
def find_min(thr_alloc) :
min, lst = thr_alloc[0]
if min == 0 :
return 0
midx = 0
for index in range(1, len(thr_alloc)) :
count, lst = thr_alloc[index]
if count < min :
min = count
midx = index
return midx
lim = int(math.sqrt(n)) + 1
nexti_lb = range(2, lim, 1)
if th < 2 :
return [nexti_lb]
thr_allocs = map(lambda i : (0, [] ), range(th))
Z = map(operations, nexti_lb)
L = zip(map(operations, nexti_lb), nexti_lb)
for i in L :
ops, index = i
mindex = find_min(thr_allocs)
cnt, lst = thr_allocs[mindex]
cnt += ops
lst.append(index)
thr_allocs[mindex] = (cnt, lst)
return map(lambda p: p[1], thr_allocs) </CODE></PRE><P>The distribution between the threads is done as follows: For each index, estimate the number of operations using the equation above. Then, assign the index to the thread with the least amount of work already assigned to it. </P><P>For example, for N = 200 the estimated number of operations per index assuming no sharing is shown in <A HREF="#smpdist"> this graph </A> . </P><P><A NAME="smpdist"><IMG SRC="http://github.com/fons/blog-images/raw/master/smp_dist.png "></A></P><P>For three threads of execution the distribution of the sieve indices given </P><P>by <EM>smp_load_balance</EM> is : </P><UL><LI><P>thread 1 : [2, 9, 12, 14], with a total number of operations of 120.</P></LI><LI><P>thread 2 : [3, 6, 8, 11], with a total number of operations of 119.</P></LI><LI><P>thread 3 : [4, 5, 7, 10, 13], with a total number of operations of 123.</P></LI></UL><P>When the sieve is shared, the amount of work remains unbalanced. Consider the example above. Thread 3 is assigned index 4, but - when the sieve is shared - there's no work to be done. </P><P>The <A HREF="#allresults"> timing test </A> shows that there is no difference<BR>
when using naive load balancing as is done in <EM>main_nolocks</EM>. In fact, the performance of both is comparable. </P><P>That's because the total number of operations across all threads of execution remains the same, and the GIL effectively prevents real parallel execution of the threads. So it really doesn't matter how unbalanced the work is. </P><H3>What about Jython and IronPython ?</H3><P>I repeated the runs <A HREF="#allresults"> shown earlier </A> on Jython 2.5. I run Jython in stand-alone mode. </P><P><IMG SRC="http://github.com/fons/blog-images/raw/master/threads_jython.png"></P><P>Note that all run times are about a factor three higher than under CPython. The Jython run times are pretty much the same. However, you still see that the locking used in the original implementation creates a performance hit. In addition, the performance doesn't depend on the number of threads either. </P><P>I attempted to repeat this experiment with IronPython. However, IronPython will not run the timing tests with garbage collection disabled. </P><H2>Conclusions</H2><P>The thread module in Python is not well suited when multi-treading compute-bound algorithms. That's nothing knew. </P><P>There are other lessons I learned from the various ways the <EM>Sieve</EM> can be distributed across threads: avoid locking, consider sharing state, and be aware how much work each thread does. </P><P>Let's start with locking. It's obvious from the <A HREF="#allresults"> results </A> that locking comes at a huge cost. Removing the locks proved to be the single most effective performance booster. </P><P>The <EM>Sieve</EM> algorithm is interesting in that shared state is quite beneficial to performance. Look at the <A HREF="#allresults"><EM>main_nolocks_alt</EM> implementation. </A> That was an attempt to eliminate shared state amongst the threads. Each thread proceeds based on it's local data only. But because the prime numbers discovered by one of the threads couldn't be shared with the others, they ended up doing more work, and performance decreased. It could be that when real parallel processing is possible, this may be balanced out by the increase in performance due to better parallelism. That however remains to be seen. </P><P>An alternative approach is to remove the global sieve variable in <EM>main_nolocks_alt</EM>. Each thread would work independently on a local seive array, resulting in a partially processed seive. The 'partial' sieves need to be combined to get the final sieve array. </P><P>One way to think of this is as a 'map' operation over the inputs, where each 'map' operation works independently. The results of the mapping phase are the combined (i.e. 'reduced') to get the final result, hence <A HREF="http://www.cs.vu.nl/~ralf/MapReduce/">map-reduce</A>. </P><P>The amount of work for each index in the <EM>Sieve</EM> algorithm varies quite a bit. So care should be taken to make sure that the work load is properly balanced amongst the threads. The effect of the (lack of) balance wasn't so obvious here, since the GIL limits parallelism.<BR>
</P><P>These last two points are going to play a part in my next post where I "Al be using the <EM>Sieve</EM> to explore the multiprocessing module. </P>
</div>
<div id="post-timestamp-id" style="visibility:hidden">
<p>3461778203 </p>
</div>
</div>
<div id="footer">
<p>
Powered by <a href="http://www.github.com/fons/cl-bliky"> cl-bliky </a>
</p>
</div>
</div><!-- end of main area -->
</body>
</html>