Fixes to GPU routines, improved tests and results

Author: Gabriele-Codega

README.md

@@ -1,6 +1,6 @@
 device: cpu
-size: 512
+size: 4096
 function: 
-  routine: matmul_numba_block_serial
-  block_size: 24
+  routine: matmul_numba_serial
+  block_size: 32
 print: False

experiments/config.yaml

@@ -135,8 +135,10 @@ def main_gpu(params: dict):
     a_d = cuda.to_device(A)
     c_d = cuda.to_device(C)
 
+    # each process at each step computes a block of C of size n_loc x ncols
+    # we set parameters for the kernel accordingly
     nthreads = bs
-    blocks_per_grid = ((n_loc + nthreads-1)//nthreads,(SIZE + nthreads-1)//nthreads)
+    blocks_per_grid = ((n_loc + nthreads-1)//nthreads,(ncols + nthreads-1)//nthreads)
     threads_per_block = (nthreads, nthreads)
 
     t_tot = 0
@@ -150,6 +152,7 @@ def main_gpu(params: dict):
 
             B_block = np.empty((n_loc,ncols), dtype=np.float64)
             B_col = np.empty((SIZE,ncols), dtype=np.float64)
+            blocks_per_grid = ((n_loc + nthreads-1)//nthreads,(ncols + nthreads-1)//nthreads)
 
         # create a contiguous block from B to communicate
         create_block(B, B_block, start, ncols)

figures/pie_1node_CPU.png

@@ -4,5 +4,5 @@ rank=$OMPI_COMM_WORLD_RANK
 
 export NUMBA_NUM_THREADS=1
 
-# kernprof -lz -o "logs/time/gpu/256_rank_$rank.lprof" scripts/run.py --config experiments/config
-valgrind --tool=cachegrind --cache-sim=yes --cachegrind-out-file="logs/memory/512_naive_rank_$rank.log" python scripts/run.py --config experiments/config
+kernprof -lz -o "out_big_$rank.lprof" scripts/run.py --config experiments/config
+# valgrind --tool=cachegrind --cache-sim=yes --cachegrind-out-file="logs/memory/512_block_rank_$rank.log" python scripts/run.py --config experiments/config

figures/pie_1node_GPU.png

@@ -3,7 +3,7 @@
 import numpy as np
 from numba import cuda
 
-from matmul import matmul, matmul_numba_gpu
+from matmul import matmul, matmul_numba_block_gpu
 from matmul.utils import create_block
 
 import mpi4py
@@ -40,11 +40,12 @@ def test_parallel_cpu():
     row_offset = np.cumsum(workloads)[rank-1] if rank > 0 else 0
 
     # initialise matrices somehow
-    A = np.arange(1, SIZE*n_loc + 1, dtype=np.float64).reshape((n_loc,SIZE)) + (row_offset * SIZE)
-    B = np.zeros((n_loc,SIZE), dtype=np.float64)
+    np.random.seed(0)
+    atot = np.random.rand(SIZE,SIZE)
+    btot = np.linalg.inv(atot)
+    A = atot[row_offset:row_offset+n_loc,:]
+    B = btot[row_offset:row_offset+n_loc,:]
     C = np.zeros((n_loc,SIZE), dtype=np.float64)
-    for i in range(n_loc):
-        B[i, i+row_offset] = 1
 
     # Compute quantities for Allgatherv and allocate required memory
     ncols = workloads[0]
@@ -80,12 +81,12 @@ def test_parallel_cpu():
     rcvcounts = workloads*SIZE
     displacements = np.cumsum(rcvcounts) - rcvcounts
     if rank == 0:
-        A_tot = np.arange(1, SIZE*SIZE + 1, dtype=np.float64).reshape((SIZE,SIZE))
+        target = np.eye(SIZE)
     C_tot = np.zeros((SIZE,SIZE))
     comm.Gatherv([C, MPI.DOUBLE], [C_tot, rcvcounts, displacements, MPI.DOUBLE])
 
     if rank == 0:
-        assert np.allclose(A_tot,C_tot)
+        assert np.allclose(target,C_tot)
 
     comm.Barrier()
 
@@ -117,11 +118,12 @@ def test_parallel_gpu():
     row_offset = np.cumsum(workloads)[rank-1] if rank > 0 else 0
 
     # initialise matrices somehow
-    A = np.arange(1, SIZE*n_loc + 1, dtype=np.float64).reshape((n_loc,SIZE)) + (row_offset * SIZE)
-    B = np.zeros((n_loc,SIZE), dtype=np.float64)
+    np.random.seed(0)
+    atot = np.random.rand(SIZE,SIZE)
+    btot = np.linalg.inv(atot)
+    A = atot[row_offset:row_offset+n_loc,:]
+    B = btot[row_offset:row_offset+n_loc,:]
     C = np.zeros((n_loc,SIZE), dtype=np.float64)
-    for i in range(n_loc):
-        B[i, i+row_offset] = 1
 
     # Compute quantities for Allgatherv and allocate required memory
     ncols = workloads[0]
@@ -151,9 +153,9 @@ def test_parallel_gpu():
 
             B_block = np.empty((n_loc,ncols), dtype=np.float64)
             B_col = np.empty((SIZE,ncols), dtype=np.float64)
-
             blocks_per_grid = ((n_loc + nthreads-1)//nthreads,(ncols + nthreads-1)//nthreads)
 
+
         # create a contiguous block from B to communicate
         create_block(B, B_block, start, ncols)
         # gather all pieces of B from other processes
@@ -163,7 +165,7 @@ def test_parallel_gpu():
         b_d = cuda.to_device(B_col)
 
         # multiply
-        matmul_numba_gpu[blocks_per_grid, threads_per_block](a_d,b_d,c_d[:,start:start+ncols])
+        matmul_numba_block_gpu[blocks_per_grid, threads_per_block](a_d,b_d,c_d[:,start:start+ncols])
 
         start += ncols
 
@@ -172,12 +174,12 @@ def test_parallel_gpu():
     rcvcounts = workloads*SIZE
     displacements = np.cumsum(rcvcounts) - rcvcounts
     if rank == 0:
-        A_tot = np.arange(1, SIZE*SIZE + 1, dtype=np.float64).reshape((SIZE,SIZE))
+        target = np.eye(SIZE)
     C_tot = np.zeros((SIZE,SIZE))
     comm.Gatherv([C, MPI.DOUBLE], [C_tot, rcvcounts, displacements, MPI.DOUBLE])
 
     if rank == 0:
-        assert np.allclose(A_tot,C_tot)
+        assert np.allclose(target,C_tot)
 
     comm.Barrier()
 

figures/pie_4nodes_CPU.png

@@ -6,59 +6,65 @@
 
 def test_matmul():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     matmul(A,B,C,None)
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 def test_matmul_numba_cpu():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     matmul_numba_cpu(A,B,C,None)
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 def test_matmul_numba_serial():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     matmul_numba_serial(A,B,C,None)
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 def test_matmul_numba_block_cpu():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     matmul_numba_block_cpu(A,B,C,6)
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 def test_matmul_numba_block_serial():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     matmul_numba_block_serial(A,B,C,6)
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 @pytest.mark.skipif((not numba.cuda.is_available()), reason='Could not find any CUDA GPU')
 def test_matmul_numba_gpu():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     a_d = numba.cuda.to_device(A)
@@ -73,13 +79,14 @@ def test_matmul_numba_gpu():
 
     C = c_d.copy_to_host()
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)
 
 @pytest.mark.skipif((not numba.cuda.is_available()), reason='Could not find any CUDA GPU')
 def test_matmul_numba_block_gpu():
     size = 20
-    A = np.arange(1,size*size+1,1,dtype=np.float64).reshape((size,size))
-    B = np.eye(size,dtype=np.float64)
+    np.random.seed(0)
+    A = np.random.rand(size,size)
+    B = np.linalg.inv(A)
     C = np.zeros((size,size),dtype=np.float64)
 
     a_d = numba.cuda.to_device(A)
@@ -94,6 +101,6 @@ def test_matmul_numba_block_gpu():
 
     C = c_d.copy_to_host()
 
-    assert np.allclose(A,C)
+    assert np.allclose(np.eye(size),C)