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
|
import numpy as np
import scipy.linalg as splinalg
import dask.array as da
import timeit
import os
import ctypes
from ctypes.util import find_library
openblas_lib = ctypes.cdll.LoadLibrary(find_library('openblas'))
def get_num_threads():
return openblas_lib.openblas_get_num_threads()
def set_num_threads(n):
openblas_lib.openblas_set_num_threads(int(n))
seed = 1234
np.random.seed(seed)
N = 1000000
p = 100
X = np.random.random(N * p).reshape((N, p), order='F')
XT = X.T.copy()
true_value=33334547.40257686
#X = da.from_array(X, chunks=(N/4, p))
old_num_threads = get_num_threads()
def test():
if not np.isclose(np.trace(X.T.dot(X)), true_value):
raise ValueError()
def test2():
if not np.isclose(np.trace(splinalg.blas.dsyrk(1., X, trans=1)), true_value):
raise ValueError()
def test3():
if not np.isclose(np.trace(XT.dot(X)), true_value):
raise ValueError()
t = timeit.timeit(test, number=5)
print("Multi threaded computation with {} threads: {}".format(old_num_threads, t))
t = timeit.timeit(test2, number=5)
print("Multi threaded computation dsyrk with {} threads: {}".format(old_num_threads, t))
t = timeit.timeit(test3, number=5)
print("Multi threaded computation dgemv with {} threads: {}".format(old_num_threads, t))
set_num_threads(1)
t = timeit.timeit(test, number=5)
print("Non multi-threaded computation:{}".format(t))
#set_num_threads(old_num_threads)
print("using dask array")
for Nchunk in [1, 2, 4]:
X = da.from_array(X, chunks=(N / Nchunk, p))
def test_dask():
if not np.isclose(np.trace(X.T.dot(X).compute()), true_value, Nchunk):
raise ValueError()
t = timeit.timeit(test_dask, number=5)
print("Dask computation {} chunk: {}".format(Nchunk, t))
|
