<jupyter_start><jupyter_text>--- **LAB 4 - CUDA in Python**---  ▶️ CUDA setup<jupyter_code>!nvcc --version
!nvidia-smi<jupyter_output>Tue Mar 29 13:50:17 2022       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 460.32.03    Driver Version: 460.32.03    CUDA Version: 11.2     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  Tesla K80           Off  | 00000000:00:04.0 Off |                    0 |
| N/A   38C    P8    28W / 149W |      0MiB / 11441MiB |      0%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+
                                                                               
+-------[...]<jupyter_text>[GPU Compute Capability](https://developer.nvidia.com/cuda-gpus)<jupyter_code>from numba import cuda

# getting device information
cuda.detect();
from numba import cuda

attribs = ['COMPUTE_CAPABILITY', 'MULTIPROCESSOR_COUNT', 'MAX_REGISTERS_PER_BLOCK', 
         'MAX_SHARED_MEMORY_PER_BLOCK', ]
device = cuda.get_current_device()
#attribs = [s for s in dir(device) if s.isupper()]        # all attribs
for attr in attribs:
    print(attr, '=', getattr(device, attr))
#@title Bash setup
%%writefile /root/.bashrc

# If not running interactively, don't do anything
[ -z "$PS1" ] && return

# don't put duplicate lines in the history. See bash(1) for more options
# ... or force ignoredups and ignorespace
HISTCONTROL=ignoredups:ignorespace

# append to the history file, don't overwrite it
shopt -s histappend

# for setting history length see HISTSIZE and HISTFILESIZE in bash(1)
HISTSIZE=10000
HISTFILESIZE=20000

# check the window size after each command and, if necessary,
# update the values of LINES and COLUMNS.
shopt -s checkwinsize

# make less more friendly for non-text input files, see lesspipe(1)
[ -x /usr/bin/lesspipe ] && eval "$(SHELL=/bin/sh lesspipe)"

PS1='\[\033[01;34m\]\w\[\033[00m\]\$ '

# enable color support of ls and also add handy aliases
if [ -x /usr/bin/dircolors ]; then
    test -r ~/.dircolors && eval "$(dircolors -b ~/.dircolors)" || eval "$(dircolors -b)"
    alias ls='ls --color=auto'
    #alias dir='dir --color=auto'
    #alias vdir='vdir --color=auto'

    alias grep='grep --color=auto'
    alias fgrep='fgrep --color=auto'
    alias egrep='egrep --color=auto'
fi

# some more ls aliases
alias ll='ls -lF'
alias la='ls -A'
alias l='ls -CF'

# path setup
export PATH="./:/usr/local/cuda/bin:$PATH"

!source /root/.bashrc<jupyter_output>Overwriting /root/.bashrc<jupyter_text>Clone GPUcomputing site on github...<jupyter_code>!git clone https://github.com/giulianogrossi/GPUcomputing.git<jupyter_output>Cloning into 'GPUcomputing'...
remote: Enumerating objects: 158, done.[K
remote: Counting objects: 100% (158/158), done.[K
remote: Compressing objects: 100% (106/106), done.[K
remote: Total 158 (delta 63), reused 133 (delta 38), pack-reused 0[K
Receiving objects: 100% (158/158), 1.19 MiB | 21.37 MiB/s, done.
Resolving deltas: 100% (63/63), done.<jupyter_text>🐍 Add kernel<jupyter_code>from numba import cuda
import numpy as np
# CUDA kernels written with `@cuda.jit` do not return values
@cuda.jit
def add_kernel(x, y, z):
  """Perform vector sum z = x * y
  """
    
  # The actual values of the following CUDA-provided variables for thread and block indices,
  # like function parameters, are not known until the kernel is launched.
  
  # This calculation gives a unique thread index within the entire grid (see the slides above for more)
  idx = cuda.grid(1)      # 1 = one dimensional thread grid, returns a single value.
                          # This Numba-provided convenience function is equivalent to
                          # `cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x`

  # This thread will do the work on the data element with the same index as its own
  # unique index within the grid.
  z[idx] = x[idx] + y[idx]<jupyter_output><empty_output><jupyter_text>Define the kernel parameters...<jupyter_code># Because of how we wrote the kernel above, we need to have a 1 thread to one data element mapping,
# therefore we define the number of threads in the grid (128*32) to equal n (4096).
threads_per_block = 128
blocks_per_grid = 32<jupyter_output><empty_output><jupyter_text>Generate data and device setup/transfert...<jupyter_code># populate data
n = 4096
x = np.arange(n).astype(np.int32) # [0...4095] on the host
y = np.ones_like(x)               # [1...1] on the host

d_x = cuda.to_device(x) # Copy of x on the device
d_y = cuda.to_device(y) # Copy of y on the device
d_out = cuda.device_array_like(d_x) # Like np.array_like, but for device arrays<jupyter_output><empty_output><jupyter_text>launch the kernel...<jupyter_code>%%time
add_kernel[blocks_per_grid, threads_per_block](d_x, d_y, d_out)
cuda.synchronize()
print(d_out.copy_to_host()) # Should be [1...4096]<jupyter_output>[   1    2    3 ... 4094 4095 4096]
CPU times: user 140 ms, sys: 4.01 ms, total: 144 ms
Wall time: 147 ms<jupyter_text>🐍 Mat multiplication kernel<jupyter_code>from numba import cuda
import numpy as np
@cuda.jit
def matmul_gpu(A, B, C):
  """Perform square matrix multiplication of C = A * B
  """

  row, col = cuda.grid(2)
  if row < C.shape[0] and col < C.shape[1]:
    tmp = 0.
    for k in range(A.shape[1]):
      tmp += A[row, k] * B[k, col]
    C[row, col] = tmp

# matrix multiplication using the cpu and C-like programming style (** strongly discouraged **)
def matmul_cpu(A, B):
  """Perform square matrix multiplication of C = A * B
  """
  C = np.zeros((A.shape[1],B.shape[1]))
  for i in range(A.shape[1]):
    for j in range(B.shape[0]):
      tmp = 0.                            
      for k in range(A.shape[1]):
        tmp += A[i, k] * B[k, j]   # multiply elements in row i of A and column j of B and add to temp
      C[i, j] = tmp    
# generate random vals
np.random.seed(42)
SIZE = 1000
A = np.ones((SIZE,SIZE)).astype('float32')  # mat 1
B = np.ones((SIZE,SIZE)).astype('float32')  # mat 2
C = np.zeros((SIZE,SIZE)).astype('float32')                       # mat where we store answer 

d_A = cuda.to_device(A) # Copy of A on the device
d_B = cuda.to_device(B) # Copy of B on the device
d_C = cuda.device_array_like(A) # malloc on the device

# data type
d_A.dtype
!nvidia-smi
print(C)<jupyter_output>[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]<jupyter_text>Kernel parameters...<jupyter_code>block = (16, 16)  # each block will contain 16x16 threads, typically 128 - 512 threads/block
grid_x = int(np.ceil(C.shape[0] / block[0]))
grid_y = int(np.ceil(C.shape[1] / block[1]))
grid = (grid_x, grid_y)  # we calculate the gridsize (number of blocks) from array
print(grid)
print(f"The kernel will be executed up to element {block[0]*grid_x}")
%%time
# execution of the kernel
matmul_gpu[grid, block](d_A, d_B, d_C)
cuda.synchronize()
C = d_C.copy_to_host()
print(C)
%%time
# execution of the python function in C-like programming style (strongly discouraged)
D = matmul_cpu(A, B)
%%time
# execution of the kernel
C = np.dot(A, B)
print(C)<jupyter_output>[[1000. 1000. 1000. ... 1000. 1000. 1000.]
 [1000. 1000. 1000. ... 1000. 1000. 1000.]
 [1000. 1000. 1000. ... 1000. 1000. 1000.]
 ...
 [1000. 1000. 1000. ... 1000. 1000. 1000.]
 [1000. 1000. 1000. ... 1000. 1000. 1000.]
 [1000. 1000. 1000. ... 1000. 1000. 1000.]]<jupyter_text>🐍 Histogram<jupyter_code>from numba import cuda
import numpy as np<jupyter_output><empty_output><jupyter_text>🔴 TODO<jupyter_code>@cuda.jit
def BMP_hist(I, H):
  row, col = cuda.grid(2)
  if row < I.shape[0] and col < I.shape[1]:
      rgb = I[row, col]
      cuda.atomic.add(H,(0, I[row,col,0]),1)  #H[0, I[row, col, 0]] += 1
      cuda.atomic.add(H,(1, I[row,col,1]),1)  #H[1, I[row, col, 1]] += 1
      cuda.atomic.add(H,(2, I[row,col,2]),1)  #H[2, I[row, col, 2]] += 1<jupyter_output><empty_output><jupyter_text>Gestione immagine e kernel...<jupyter_code>from PIL import Image

#image.save('beach1.bmp')
img = Image.open('GPUcomputing/images/dog.bmp')
img_rgb = img.convert('RGB')
img
# converto to numpy (host)
I = np.asarray(img)
print(f"Image size W x H x ch = {I.shape}")
# device data setup
d_I = cuda.to_device(I)
H = np.zeros((3,256)).astype(np.float32)
d_H = cuda.to_device(H)

print(H[0,I[0,0,0]])
H[0, I[0,0,0]] += 1
print(H[0,I[0,0,0]])
H[0, I[0,0,0]] += 1
print(H[0,I[0,0,0]])<jupyter_output>Image size W x H x ch = (524, 1024, 3)
0.0
1.0
2.0<jupyter_text>Kernel parameters...<jupyter_code>block = (16, 16)  # each block will contain 16x16 threads, typically 128 - 512 threads/block
grid_x = int(np.ceil(I.shape[0] / block[0]))
grid_y = int(np.ceil(I.shape[1] / block[1]))
grid = (grid_x, grid_y)  # we calculate the gridsize (number of blocks) from array
print(grid)
print(f"The kernel will be executed up to element {block[0]*grid_x}")

# kernel launch

BMP_hist[grid, block](d_I, d_H)
hist = d_H.copy_to_host()
print(hist.shape)

import matplotlib.pyplot as plt

plt.bar(np.arange(256),hist[0])
plt.title('Histogram (R)')
plt.show()
plt.bar(np.arange(256),hist[1])
plt.title('Histogram (G)')
plt.show()
plt.bar(np.arange(256),hist[2])
plt.title('Histogram (B)')
plt.show()<jupyter_output><empty_output><jupyter_text>🔴 TODO Calcolare il prodotto di matrici MQDB con kernel CUDA in Python<jupyter_code>import numpy as np
import random
from numba import cuda



#n = dimension of the matrix
#k = number of submatrices
#const_ro_rand = how to fill the submatrices
#const_value = the constant value (if needed)
def generate_mqdb(n, k, const_or_rand, const_value = 0):

  if k > n: return

  mat = np.zeros((n,n))

  sub_dims = []

  tmp = 0
  remaining_values = k
  for i in range(k-1):
    curr_subdim = random.randint(1, n-remaining_values+1 - tmp)
    sub_dims.append(curr_subdim)
    for j in range(tmp, tmp + curr_subdim):
      for k in range(tmp, tmp + curr_subdim):  
        if const_or_rand == "C":
          mat[j,k] = const_value
        else: 
          mat[j,k] = random.uniform(0,1)
        
    tmp += curr_subdim
    remaining_values -= 1

  final_dim = n - tmp
  sub_dims.append(final_dim)
  for j in range(tmp, tmp + final_dim):
    for k in range(tmp, tmp + final_dim):  
      if const_or_rand == "C":
        mat[j,k] = const_value
      else: 
         mat[j,k] = random.uniform(0,1)

  return (mat, sub_dims)


class my_mqdb:
  def __init__(self, n, k, const_or_rand, const_value = 0):
    (self.mat, self.sub_matrices_dims) = generate_mqdb(n,k, const_or_rand, const_value)


#A,B = input matrices
#C = output matrices
#cur_k = dimension of the current submatrix
#offset = from where to start the calculus
@cuda.jit
def BMP_hist(I, H):
  row, col = cuda.grid(2)
  if row < I.shape[0] and col < I.shape[1]:
      rgb = I[row, col]
      H[0, rgb[0]] += 1
      H[1, rgb[1]] += 1
      H[2, rgb[2]] += 1

@cuda.jit
def mqdb_smart(A, B, C, cur_k, offset):
  x, y = cuda.grid(2)
  if x < cur_k and y < cur_k:
    sum = 0
    for r in range(cur_k):
      sum = sum + A[x+offset, r+offset] * B[r+offset, y+offset]
    C[x+offset,y+offset] = sum


#A = my_mqdb(100,4,"C", 5)
#B = my_mqdb(100,4,"C", 5)
#B.mat = np.copy(A.mat)
#B.sub_matrices_dims = A.sub_matrices_dims

n = 10
k = 3

A, k_arr = generate_mqdb(n,k,"C", 5)
B = np.copy(A)
C = np.zeros((n,n))
d_A = cuda.to_device(A)
d_B = cuda.to_device(B)
d_C = cuda.to_device(C)



block = (16, 16)  # each block will contain 16x16 threads, typically 128 - 512 threads/block
off = 0

for elem in k_arr:
  grid_x = int(np.ceil(elem / block[0]))
  grid_y = int(np.ceil(elem / block[1]))
  grid = (grid_x, grid_y)
  print("" + str(grid) + ", " + str(elem))
  mqdb_smart[grid, block](d_A, d_B, d_C, elem, off)
  off += elem

C = d_C.copy_to_host()
print(A)
print(B)
print(C)<jupyter_output>(1, 1), 7
(1, 1), 1
(1, 1), 2
[[5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 5. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 5. 5.]
 [0. 0. 0. 0. 0. 0. 0. 0. 5. 5.]]
[[5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [5. 5. 5. 5. 5. 5. 5. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 5. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 5. 5.]
 [0. 0. 0. 0. 0. 0. 0. 0. 5. 5.]]
[[175. 175. 175. 175. 175. 175. 175.   0.   0.   0.]
 [175. 175. 175. 175. 175. 175. 175.   0.   0.   0.]
 [175. 175. 175. 175. 175. 175. 175.   0.   0.   0.]
 [175. 175. 175. 175. 175. 175. 175.   0.   0.   0.]
 [175. 175. 175. 175. 175. 175. 175.   0.   0.   0.]
 [175. 175. 175. 175. 175. 175. 175.   0.  [...]
{'fraction_non_alphanumeric': 0.1822366435255056, 'fraction_numerical': 0.07032900694234832, 'mean_word_length': 2.5984251968503935, 'pattern_counts': {'":': 2, '<': 57, '<?xml version=': 0, '>': 50, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '647985', 'license': 'MIT', 'n_tokens_mistral': 6006, 'n_tokens_neox': 5216, 'n_words': 1550, 'path': 'GPU_computing/CUDA_lab4_TODO.ipynb', 'repo_name': 'katema-official/Universita_magistrale_anno_1.2', 'size': '13252'}