CPU
vs
GPU

import time
import random

N = 512  # keep small — this is O(N³)

# Build two NxN matrices filled with random floats
A = [[random.random() for _ in range(N)] for _ in range(N)]
B = [[random.random() for _ in range(N)] for _ in range(N)]
C = [[0.0]  * N       for _ in range(N)]

t0 = time.perf_counter()

for i in range(N):
    for j in range(N):
        for k in range(N):
            C[i][j] += A[i][k] * B[k][j]

elapsed = time.perf_counter() - t0
print(f"Pure Python: {elapsed:.2f}s  (N={N})")
# Typical output: ~70–120 seconds for N=512

import numpy as np
import time

N = 4096  # now we can afford a much bigger matrix

A = np.random.rand(N, N).astype(np.float32)
B = np.random.rand(N, N).astype(np.float32)

# warm-up (first call can be slower due to BLAS init)
_ = np.dot(A[:64,:64], B[:64,:64])

t0 = time.perf_counter()
C = np.matmul(A, B)           # one line — BLAS handles everything
elapsed = time.perf_counter() - t0

print(f"NumPy (CPU): {elapsed*1000:.1f} ms  (N={N})")
# Typical output: ~600–1200 ms on a modern 8-core CPU

import torch
import time

N = 4096

# ── host (CPU) tensors ──────────────────────────────────
A_cpu = torch.rand(N, N, dtype=torch.float32)
B_cpu = torch.rand(N, N, dtype=torch.float32)

# ── move to device (GPU VRAM) ───────────────────────────
device = "cuda" if torch.cuda.is_available() else "cpu"
A = A_cpu.to(device)
B = B_cpu.to(device)

# warm-up — first kernel launch has JIT overhead
_ = torch.matmul(A[:64,:64], B[:64,:64])
torch.cuda.synchronize()   # wait for GPU to finish

t0 = time.perf_counter()
C = torch.matmul(A, B)
torch.cuda.synchronize()   # ← CRITICAL: GPU is async, must sync before timing
elapsed = time.perf_counter() - t0

print(f"PyTorch GPU: {elapsed*1000:.1f} ms  (N={N})")
# Typical output: ~4–12 ms on a modern GPU (RTX 3090 / A100)
# That's ~100× faster than NumPy CPU

import threading, time, random

N = 128
A = [[random.random() for _ in range(N)] for _ in range(N)]
B = [[random.random() for _ in range(N)] for _ in range(N)]
C = [[0.0] * N for _ in range(N)]

def compute_row(i):
    for j in range(N):
        s = 0.0
        for k in range(N):
            s += A[i][k] * B[k][j]
        C[i][j] = s

t0 = time.perf_counter()
threads = [threading.Thread(target=compute_row, args=(i,)) for i in range(N)]
for t in threads: t.start()
for t in threads: t.join()
print(f"Threading: {time.perf_counter()-t0:.2f}s")
# ⚠ No speedup — GIL serialises all threads on pure-Python loops
# May even be SLOWER than serial due to lock contention

import os, time, random
from concurrent.futures import ProcessPoolExecutor

N = 256
A = [[random.random() for _ in range(N)] for _ in range(N)]
B = [[random.random() for _ in range(N)] for _ in range(N)]

# Each worker gets ONE row index + the full A, B matrices
# B is pickled and sent to every worker — the main overhead
def compute_row(args):
    i, A_w, B_w, n = args
    row = [0.0] * n
    for j in range(n):
        s = 0.0
        for k in range(n):
            s += A_w[i][k] * B_w[k][j]
        row[j] = s
    return i, row

n_cores = os.cpu_count()   # use ALL cores on this machine
print(f"Launching {n_cores} workers")

t0 = time.perf_counter()
C = [[0.0] * N for _ in range(N)]

with ProcessPoolExecutor(max_workers=n_cores) as ex:
    args = [(i, A, B, N) for i in range(N)]
    # chunksize batches rows to reduce pickling round-trips
    for i, row in ex.map(compute_row, args, chunksize=N//n_cores):
        C[i] = row

elapsed = time.perf_counter() - t0
print(f"ProcessPool ({n_cores} cores): {elapsed:.2f}s")
# Typical speedup: ~(n_cores - 1)× after pickling overhead
# 8 cores → ~5–6× faster than serial pure-Python

import numpy as np
import time

# See what BLAS backend numpy is using
np.show_config()   # look for openblas_info or blas_mkl_info

N = 4096
A = np.random.rand(N, N).astype(np.float32)
B = np.random.rand(N, N).astype(np.float32)

t0 = time.perf_counter()
C = np.matmul(A, B)
# ↑ Internally: OpenBLAS tiles A&B, spawns N pthreads,
#   each thread works a band of tiles with AVX2 SIMD.
#   No GIL. No pickling. All in C/assembly.
elapsed = time.perf_counter() - t0
print(f"NumPy matmul: {elapsed*1000:.1f} ms")

# Control how many cores OpenBLAS uses:
# OMP_NUM_THREADS=4 python matmul_numpy_parallel.py
# Or at runtime:
import os
os.environ["OMP_NUM_THREADS"] = "4"   # must be set BEFORE importing numpy