COMPASS  5.0.0
End-to-end AO simulation tool using GPU acceleration
test_magma.py
1 import numpy as np
2 import carmaWrap as ch
3 import numpy.testing as npt
4 import time
5 
6 m = 1024
7 n = 1024
8 min_mn = min(m, n)
9 
10 dec = 3
11 prec = 10**(-dec)
12 
13 c = ch.context.get_instance()
14 
15 print("precision: ", prec)
16 
17 
18 def test_float_gemv():
19 
20  #function gemv
21  # testing: y=A.x
22  # x and y are vector, A a matrix
23 
24  alpha = 2
25  beta = 1
26 
27  Mat = ch.obj_float(c, np.random.randn(m, n))
28  MatT = ch.obj_float(c, np.random.randn(n, m))
29  # Mat.random_host(seed, 'U')
30  # MatT.random_host(seed + 2, 'U')
31 
32  Vectx = ch.obj_float(c, np.random.randn(n))
33  Vecty = ch.obj_float(c, np.random.randn(m))
34  # Vectx.random_host(seed + 3, 'U')
35  # Vecty.random_host(seed + 4, 'U')
36 
37  A = np.array(Mat)
38  AT = np.array(MatT)
39  x = np.array(Vectx)
40  y = np.array(Vecty)
41 
42  y = alpha * A.dot(x) + beta * y
43  y2 = alpha * A.dot(x)
44  y3 = alpha * AT.T.dot(x)
45 
46  # Vecty = ch.obj_float(c, np.random.randn((m)))
47 
48  Mat.magma_gemv(Vectx, alpha, 'N', Vecty, beta)
49  Vecty_2 = Mat.magma_gemv(Vectx, alpha, 'N')
50  Vecty_3 = MatT.magma_gemv(Vectx, alpha, 'T')
51 
52  npt.assert_array_almost_equal(x, np.array(Vectx), decimal=dec - 1)
53  npt.assert_array_almost_equal(y, np.array(Vecty), decimal=dec - 1)
54  npt.assert_array_almost_equal(y2, np.array(Vecty_2), decimal=dec - 1)
55  npt.assert_array_almost_equal(y3, np.array(Vecty_3), decimal=dec - 1)
56 
57 
58 def test_double_gemv():
59 
60  #function gemv
61  # testing: y=A.x
62  # x and y are vector, A a matrix
63 
64  alpha = 2
65  beta = 1
66 
67  Mat = ch.obj_double(c, np.random.randn(m, n))
68  MatT = ch.obj_double(c, np.random.randn(n, m))
69  # Mat.random_host(seed, 'U')
70  # MatT.random_host(seed + 2, 'U')
71 
72  Vectx = ch.obj_double(c, np.random.randn(n))
73  Vecty = ch.obj_double(c, np.random.randn(m))
74  # Vectx.random_host(seed + 3, 'U')
75  # Vecty.random_host(seed + 4, 'U')
76 
77  A = np.array(Mat)
78  AT = np.array(MatT)
79  x = np.array(Vectx)
80  y = np.array(Vecty)
81 
82  y = alpha * A.dot(x) + beta * y
83  y2 = alpha * A.dot(x)
84  y3 = alpha * AT.T.dot(x)
85 
86  # Vecty = ch.obj_double(c, np.random.randn((m)))
87 
88  Mat.magma_gemv(Vectx, alpha, 'N', Vecty, beta)
89  Vecty_2 = Mat.magma_gemv(Vectx, alpha, 'N')
90  Vecty_3 = MatT.magma_gemv(Vectx, alpha, 'T')
91 
92  npt.assert_array_almost_equal(x, np.array(Vectx), decimal=dec - 1)
93  npt.assert_array_almost_equal(y, np.array(Vecty), decimal=dec - 1)
94  npt.assert_array_almost_equal(y2, np.array(Vecty_2), decimal=dec - 1)
95  npt.assert_array_almost_equal(y3, np.array(Vecty_3), decimal=dec - 1)
96 
97 
98 def test_magma_float_svd():
99 
100  mat = np.random.rand(m, n).astype(np.float32)
101 
102  h_mat = ch.host_obj_float(mat, ch.MA_PAGELOCK)
103  h_eig = ch.host_obj_float(np.random.randn(min_mn), ch.MA_PAGELOCK)
104  h_U = ch.host_obj_float(np.random.randn(m, m), ch.MA_PAGELOCK)
105  h_VT = ch.host_obj_float(np.random.randn(n, n), ch.MA_PAGELOCK)
106 
107  npt.assert_array_equal(mat, np.array(h_mat))
108 
109  ch.magma_svd_cpu_float(h_mat, h_eig, h_U, h_VT)
110 
111  # expected: U.S.V=mat
112  # U = np.array(h_U)
113  # S = np.random.randn((m, n), dtype=np.float32)
114  # S = np.diag(h_eig)
115  # VT = np.array(h_VT)
116 
117  # res = np.dot(U, S)
118  # res = np.dot(res, VT.T)
119 
120  # iErr = np.argmax(np.abs(res - mat))
121  # err = np.abs(np.abs(res.item(iErr) - mat.item(iErr)))
122  # if (mat.item(iErr) != 0):
123  # err = err / mat.item(iErr)
124 
125  # npt.assert_almost_equal(err, 0., decimal=dec)
126  # print("")
127  # print(err)
128 
129  p_U, p_S, p_V = np.linalg.svd(mat)
130  npt.assert_array_almost_equal(np.array(h_eig), p_S, decimal=dec)
131 
132 
133 def test_magma_float_getri_cpu():
134 
135  a = np.random.rand(m, m).astype(np.float32)
136  a = np.dot(a, a.T)
137 
138  a += np.identity(m, dtype=np.float32)
139 
140  h_mat = ch.host_obj_float(a, ch.MA_PAGELOCK)
141  h_mat2 = ch.host_obj_float(a, ch.MA_PAGELOCK)
142 
143  ch.magma_getri_cpu_float(h_mat2)
144 
145  a_1 = np.array(h_mat2)
146  #expected: a.a_1=Id
147  res = np.dot(a, a_1)
148 
149  err = np.amax(np.abs(res - np.identity(m)))
150 
151  npt.assert_almost_equal(err, 0., decimal=dec)
152 
153 
154 def test_magma_float_potri_cpu():
155 
156  a = np.random.rand(m, m).astype(np.float32)
157  a = np.dot(a, a.T)
158  a += np.identity(m, dtype=np.float32)
159 
160  h_mat = ch.host_obj_float(a, ch.MA_PAGELOCK)
161  h_mat2 = ch.host_obj_float(a, ch.MA_PAGELOCK)
162 
163  ch.magma_potri_cpu_float(h_mat2)
164 
165  a_1 = np.array(h_mat2)
166  #expected: a.a_1=Id
167  res = np.dot(a, a_1)
168 
169  err = np.amax(np.abs(res - np.identity(m)))
170  print("")
171  print(err)
172 
173  npt.assert_almost_equal(err, 0., decimal=dec)
174 
175 
176 def test_magma_float_getri_gpu():
177 
178  d_mat = ch.obj_float(c, np.random.randn(m, m))
179  d_mat.random(np.int32(time.perf_counter() * 1e3))
180  a = np.array(d_mat)
181  a = a.T
182  a.reshape(a.T.shape)
183 
184  identity = ch.obj_float(c, np.identity(m))
185 
186  d_res = d_mat.gemm(d_mat, 'n', 't', 1, identity, 1)
187 
188  b = np.dot(a, a.T)
189  mat = np.array(d_res)
190 
191  d_mat = ch.obj_float(c, d_res)
192  ch.magma_getri_float(d_res)
193 
194  d_id = d_mat.gemm(d_res)
195  res_id = np.array(d_id)
196 
197  err = np.amax(np.abs(res_id - np.identity(m)))
198 
199  npt.assert_almost_equal(err, 0., decimal=dec)
200 
201 
202 def test_magma_float_potri_gpu():
203 
204  d_mat = ch.obj_float(c, np.random.randn(m, m))
205  d_mat.random(np.int32(time.perf_counter() * 1e3))
206  a = np.array(d_mat)
207  a = a.T
208  a.reshape(a.T.shape)
209 
210  identity = ch.obj_float(c, np.identity(m))
211 
212  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t', beta=1, matC=identity)
213 
214  b = np.dot(a, a.T)
215  mat = np.array(d_res)
216 
217  d_mat = ch.obj_float(c, d_res)
218  ch.magma_potri_float(d_res)
219 
220  d_id = d_mat.gemm(d_res)
221  res_id = np.array(d_id)
222 
223  err = np.amax(np.abs(res_id - np.identity(m)))
224  print("")
225  print(err)
226 
227  npt.assert_almost_equal(err, 0, decimal=dec)
228 
229 
230 def test_magma_float_syevd():
231 
232  d_mat = ch.obj_float(c, np.random.randn(m, m))
233  d_U = ch.obj_float(c, np.random.randn(m, m))
234  h_EV = ch.host_obj_float(np.random.randn(m))
235  h_EV2 = ch.host_obj_float(np.random.randn(m))
236 
237  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t')
238 
239  ch.magma_syevd_float(d_res, h_EV, d_U)
240 
241  U = np.array(d_U)
242  Mat = np.array(d_res)
243  EV = np.diag(h_EV)
244 
245  npt.assert_almost_equal(
246  np.dot(np.dot(U, EV), U.T).astype(np.float32), Mat, decimal=dec - 1)
247 
248  err = np.amax(np.abs(Mat - np.dot(np.dot(U, EV), U.T)))
249 
250  print("")
251 
252  print("out of place, compute U")
253  print(err)
254 
255  npt.assert_almost_equal(err, 0., decimal=dec - 1)
256 
257  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t')
258  ch.magma_syevd_float(d_res, h_EV2, computeU=False)
259 
260  err = np.amax(np.abs(np.array(h_EV) - np.array(h_EV2)))
261 
262  print("in place, U not computed")
263  print(err)
264  npt.assert_almost_equal(np.array(h_EV), np.array(h_EV2), decimal=dec - 2)
265 
266 
267 def test_magma_double_svd():
268 
269  mat = np.random.rand(m, n).astype(np.float32)
270 
271  h_mat = ch.host_obj_double(mat, ch.MA_PAGELOCK)
272  h_eig = ch.host_obj_double(np.random.randn(min_mn), ch.MA_PAGELOCK)
273  h_U = ch.host_obj_double(np.random.randn(m, m), ch.MA_PAGELOCK)
274  h_VT = ch.host_obj_double(np.random.randn(n, n), ch.MA_PAGELOCK)
275 
276  npt.assert_array_equal(mat, np.array(h_mat))
277 
278  ch.magma_svd_cpu_double(h_mat, h_eig, h_U, h_VT)
279 
280  # expected: U.S.V=mat
281  # U = np.array(h_U)
282  # S = np.random.randn((m, n), dtype=np.float32)
283  # S = np.diag(h_eig)
284  # VT = np.array(h_VT)
285 
286  # res = np.dot(U, S)
287  # res = np.dot(res, VT.T)
288 
289  # iErr = np.argmax(np.abs(res - mat))
290  # err = np.abs(np.abs(res.item(iErr) - mat.item(iErr)))
291  # if (mat.item(iErr) != 0):
292  # err = err / mat.item(iErr)
293 
294  # npt.assert_almost_equal(err, 0., decimal=dec)
295  # print("")
296  # print(err)
297 
298  p_U, p_S, p_V = np.linalg.svd(mat)
299  npt.assert_array_almost_equal(np.array(h_eig), p_S, decimal=2 * dec - 1)
300 
301 
302 def test_magma_double_getri_cpu():
303 
304  a = np.random.rand(m, m).astype(np.float32)
305  a = np.dot(a, a.T)
306 
307  a += np.identity(m, dtype=np.float32)
308 
309  h_mat = ch.host_obj_double(a, ch.MA_PAGELOCK)
310  h_mat2 = ch.host_obj_double(a, ch.MA_PAGELOCK)
311 
312  ch.magma_getri_cpu_double(h_mat2)
313 
314  a_1 = np.array(h_mat2)
315  #expected: a.a_1=Id
316  res = np.dot(a, a_1)
317 
318  err = np.amax(np.abs(res - np.identity(m)))
319 
320  npt.assert_almost_equal(err, 0., decimal=2 * dec)
321 
322 
323 def test_magma_double_potri_cpu():
324 
325  a = np.random.rand(m, m).astype(np.float32)
326  a = np.dot(a, a.T)
327  a += np.identity(m, dtype=np.float32)
328 
329  h_mat = ch.host_obj_double(a, ch.MA_PAGELOCK)
330  h_mat2 = ch.host_obj_double(a, ch.MA_PAGELOCK)
331 
332  ch.magma_potri_cpu_double(h_mat2)
333 
334  a_1 = np.array(h_mat2)
335  #expected: a.a_1=Id
336  res = np.dot(a, a_1)
337 
338  err = np.amax(np.abs(res - np.identity(m)))
339  print("")
340  print(err)
341 
342  npt.assert_almost_equal(err, 0., decimal=2 * dec)
343 
344 
345 def test_magma_double_getri_gpu():
346 
347  d_mat = ch.obj_double(c, np.random.randn(m, m))
348  d_mat.random(np.int32(time.perf_counter() * 1e3))
349  a = np.array(d_mat)
350  a = a.T
351  a.reshape(a.T.shape)
352 
353  identity = ch.obj_double(c, np.identity(m))
354 
355  d_res = d_mat.gemm(d_mat, 'n', 't', 1, identity, 1)
356 
357  b = np.dot(a, a.T)
358  mat = np.array(d_res)
359 
360  d_mat = ch.obj_double(c, d_res)
361  ch.magma_getri_double(d_res)
362 
363  d_id = d_mat.gemm(d_res)
364  res_id = np.array(d_id)
365 
366  err = np.amax(np.abs(res_id - np.identity(m)))
367 
368  npt.assert_almost_equal(err, 0., decimal=2 * dec)
369 
370 
371 def test_magma_double_potri_gpu():
372 
373  d_mat = ch.obj_double(c, np.random.randn(m, m))
374  d_mat.random(np.int32(time.perf_counter() * 1e3))
375  a = np.array(d_mat)
376  a = a.T
377  a.reshape(a.T.shape)
378 
379  identity = ch.obj_double(c, np.identity(m))
380 
381  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t', beta=1, matC=identity)
382 
383  b = np.dot(a, a.T)
384  mat = np.array(d_res)
385 
386  d_mat = ch.obj_double(c, d_res)
387  ch.magma_potri_double(d_res)
388 
389  d_id = d_mat.gemm(d_res)
390  res_id = np.array(d_id)
391 
392  err = np.amax(np.abs(res_id - np.identity(m)))
393  print("")
394  print(err)
395 
396  npt.assert_almost_equal(err, 0, decimal=2 * dec)
397 
398 
399 def test_magma_double_syevd():
400 
401  d_mat = ch.obj_double(c, np.random.randn(m, m))
402  d_U = ch.obj_double(c, np.random.randn(m, m))
403  h_EV = ch.host_obj_double(np.random.randn(m))
404  h_EV2 = ch.host_obj_double(np.random.randn(m))
405 
406  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t')
407 
408  ch.magma_syevd_double(d_res, h_EV, d_U)
409 
410  U = np.array(d_U)
411  Mat = np.array(d_res)
412  EV = np.diag(h_EV)
413 
414  npt.assert_almost_equal(
415  np.dot(np.dot(U, EV), U.T).astype(np.float32), Mat, decimal=2 * dec - 2)
416 
417  err = np.amax(np.abs(Mat - np.dot(np.dot(U, EV), U.T)))
418 
419  print("")
420 
421  print("out of place, compute U")
422  print(err)
423 
424  npt.assert_almost_equal(err, 0., decimal=2 * dec - 1)
425 
426  d_res = d_mat.gemm(d_mat, op_a='n', op_b='t')
427  ch.magma_syevd_double(d_res, h_EV2, computeU=False)
428 
429  err = np.amax(np.abs(np.array(h_EV) - np.array(h_EV2)))
430 
431  print("in place, U not computed")
432  print(err)
433  npt.assert_almost_equal(np.array(h_EV), np.array(h_EV2), decimal=2 * dec - 2)
test_magma.test_magma_double_getri_gpu
def test_magma_double_getri_gpu()
Definition: test_magma.py:345
test_magma.test_magma_float_potri_cpu
def test_magma_float_potri_cpu()
Definition: test_magma.py:154
test_magma.test_magma_float_getri_gpu
def test_magma_float_getri_gpu()
Definition: test_magma.py:176
test_magma.test_magma_float_potri_gpu
def test_magma_float_potri_gpu()
Definition: test_magma.py:202
test_magma.test_magma_double_svd
def test_magma_double_svd()
Definition: test_magma.py:267
test_magma.test_magma_float_getri_cpu
def test_magma_float_getri_cpu()
Definition: test_magma.py:133
test_magma.test_magma_double_getri_cpu
def test_magma_double_getri_cpu()
Definition: test_magma.py:302
test_magma.test_magma_double_potri_cpu
def test_magma_double_potri_cpu()
Definition: test_magma.py:323
test_magma.test_float_gemv
def test_float_gemv()
Definition: test_magma.py:18
test_magma.test_magma_double_potri_gpu
def test_magma_double_potri_gpu()
Definition: test_magma.py:371
test_magma.test_magma_float_svd
def test_magma_float_svd()
Definition: test_magma.py:98
test_magma.test_magma_double_syevd
def test_magma_double_syevd()
Definition: test_magma.py:399
test_magma.test_double_gemv
def test_double_gemv()
Definition: test_magma.py:58
test_magma.test_magma_float_syevd
def test_magma_float_syevd()
Definition: test_magma.py:230