# VGPU Persistent Compute Pattern ## Compile Once, Dispatch Many: A Compiled-Kernel Cache for GLSL/CUDA **Version 1.0** --- ## 0. The Core Insight GPU programs are objects you **construct once** (compile, link, validate) and **dispatch many times** (set uniforms, upload data, dispatch). The data changes; the program never re-compiles. Three layers: | Layer | Lifetime | Operation | |---|---|---| | **Compilation** | Once per (source, signature) | `ctx.compute_shader(source)` | | **Binding** | Once per kernel instance | `A_buf.bind_to_storage_buffer(0)` | | **Dispatch** | Per call | Upload data → set uniforms → `program.run(...)` | The pathological anti-pattern (slow): ```python def compute(A, B): program = ctx.compute_shader(source) # recompile! program['N'] = ... # set buf_a = ctx.buffer(A.tobytes()) # allocate, bind, upload buf_a.bind_to_storage_buffer(0) buf_b = ctx.buffer(B.tobytes()) # same again buf_b.bind_to_storage_buffer(1) program.run(...) return result # tiny timing benefit ``` The correct pattern (fast): ```python program = ctx.compute_shader(source) # COMPILE ONCE, KEEP ALIVE buf_a = ctx.buffer(reserve=N*N*4); buf_a.bind_to_storage_buffer(0) # BIND ONCE buf_b = ctx.buffer(reserve=N*N*4); buf_b.bind_to_storage_buffer(1) buf_c = ctx.buffer(reserve=N*N*4); buf_c.bind_to_storage_buffer(2) def compute(A, B): # Just data updates buf_a.write(A.tobytes()) # fast buf_b.write(B.tobytes()) program.run(...) return np.frombuffer(buf_c.read(), ...).copy() ``` A function call here costs only: 1. SSBO write (DMA, GPU-bound) 2. Uniform updates (cheap) 3. Dispatch (one CPU-GPU sync) 4. Readback (DMA, GPU-bound) No shader compilation. No buffer allocation. No driver re-binding. **10-100× faster on repeat calls** vs recompile-each-time. --- ## 1. The Matmul Kernel: GLSL Compute Shader ```glsl #version 430 core // =================================================================== // C = A @ B (matmul: A[N,N] · B[N,N] = C[N,N]) // Persistent kernel: source never changes between calls. // =================================================================== layout(local_size_x = 16, local_size_y = 16) in; layout(std430, binding = 0) readonly buffer MatrixA { float A[]; }; layout(std430, binding = 1) readonly buffer MatrixB { float B[]; }; layout(std430, binding = 2) writeonly buffer MatrixC { float C[]; }; uniform int u_N; void main() { uint row = gl_GlobalInvocationID.y; uint col = gl_GlobalInvocationID.x; if (row >= uint(u_N) || col >= uint(u_N)) return; float sum = 0.0; for (uint k = 0; k < uint(u_N); k++) { sum += A[row * uint(u_N) + k] * B[k * uint(u_N) + col]; } C[row * uint(u_N) + col] = sum; } ``` > **Note:** this is the conventional CUDA-equivalent kernel. A more sophisticated version with shared-memory tiling (16×16 tile per workgroup) reaches 80-95% of cuBLAS performance. For pedagogical clarity we use the direct global-memory version. Same persistence pattern applies. --- ## 2. The VGPU PersistentKernel Class ```python """ persistent_kernel.py — VGPU pattern: compile-once, dispatch-many. """ import moderngl import numpy as np from typing import Tuple # ==================================================================== # GLSL SOURCE — Compile once, share across instances of same N # ==================================================================== MATMUL_SOURCE = r""" #version 430 core layout(local_size_x = 16, local_size_y = 16) in; layout(std430, binding = 0) readonly buffer MatrixA { float A[]; }; layout(std430, binding = 1) readonly buffer MatrixB { float B[]; }; layout(std430, binding = 2) writeonly buffer MatrixC { float C[]; }; uniform int u_N; void main() { uint row = gl_GlobalInvocationID.y; uint col = gl_GlobalInvocationID.x; if (row >= uint(u_N) || col >= uint(u_N)) return; float sum = 0.0; for (uint k = 0; k < uint(u_N); k++) { sum += A[row * uint(u_N) + k] * B[k * uint(u_N) + col]; } C[row * uint(u_N) + col] = sum; } """ # ==================================================================== # PERSISTENT KERNEL — compiled once, reused forever # ==================================================================== class PersistentKernel: """ A GLSL compute kernel that stays compiled and bound for its lifetime. Data is uploaded via SSBO writes; only uniforms and dispatch change. Idempotent: calling `feed_data` repeatedly with new matrices replaces the previous content without recompiling. """ def __init__(self, ctx: moderngl.Context, N: int): self.ctx = ctx self.N = N # COMPILE ONCE self.program = ctx.compute_shader(MATMUL_SOURCE) self.u_N = self.program['u_N'] self.u_N.value = N # ALLOCATE BUFFERS ONCE size_bytes = N * N * 4 self.A_buf = ctx.buffer(reserve=size_bytes) self.B_buf = ctx.buffer(reserve=size_bytes) self.C_buf = ctx.buffer(reserve=size_bytes) # BIND ONCE self.A_buf.bind_to_storage_buffer(0) self.B_buf.bind_to_storage_buffer(1) self.C_buf.bind_to_storage_buffer(2) # Dispatch parameters (N unchanged) self._groups = (max(1, (N + 15) // 16), max(1, (N + 15) // 16), 1) self.dispatch_count = 0 self.total_runtime_ns = 0 def feed_and_dispatch(self, A: np.ndarray, B: np.ndarray) -> np.ndarray: """ Upload new matrices and dispatch. NEITHER compiles nor rebinds. Returns C = A · B. """ assert A.shape == (self.N, self.N), f"A must be ({self.N}, {self.N})" assert B.shape == (self.N, self.N), f"B must be ({self.N}, {self.N})" assert A.dtype == np.float32 and B.dtype == np.float32, "FP32 only" # 1. Upload data (DMA-bound; ~0.01–1ms for typical N) self.A_buf.write(A.tobytes()) self.B_buf.write(B.tobytes()) # 2. Dispatch (zero CPU overhead; GPU runs in parallel) self.program.run(group_x=self._groups[0], group_y=self._groups[1], group_z=self._groups[2]) # 3. Readback (DMA-bound) result_bytes = self.C_buf.read() result = np.frombuffer(result_bytes, dtype=np.float32).reshape( self.N, self.N).copy() self.dispatch_count += 1 return result def __call__(self, A, B): return self.feed_and_dispatch(A, B) # ==================================================================== # TEST / BENCHMARK # ==================================================================== if __name__ == "__main__": import time ctx = moderngl.create_standalone_context() N = 16 A = np.random.randn(N, N).astype(np.float32) B = np.random.randn(N, N).astype(np.float32) # Reference: pure-CPU matmul expected = A @ B # --- VGPU PERSISTENT PATH --- kernel = PersistentKernel(ctx, N=N) # First call: just discovery cost (no compile, since already done in __init__) t0 = time.perf_counter_ns() C = kernel(A, B) t1 = time.perf_counter_ns() print(f"Call 1 (warm path): {(t1-t0)/1e6:.3f} ms") print(f" Max abs diff: {np.max(np.abs(C - expected)):.2e}") assert np.allclose(C, expected, atol=1e-3), "MISMATCH" # 100 more calls — show constant-time dispatch t_call1 = time.perf_counter_ns() for _ in range(100): A = np.random.randn(N, N).astype(np.float32) B = np.random.randn(N, N).astype(np.float32) C = kernel(A, B) t_callN = time.perf_counter_ns() print(f"100 calls (mean): {(t_callN-t_call1)/100/1e6:.3f} ms each") print(f"\nTotal dispatches: {kernel.dispatch_count}") print(f"Object still alive & reusable (no recompilation needed).") ``` **Performance characteristic:** - `__init__` cost: dominated by shader compile (~50-300 ms first time) - Each `feed_and_dispatch`: ~0.1-1 ms for N=16 (mostly DMA transfer) - **Crucially: zero shader recompilation overhead per call** --- ## 3. The Kernel Database — When You Have Many Kernels For 5+ different kernel variants, you need a cache: ```python class KernelDatabase: """ Cache of compiled PersistentKernels, keyed by (kernel_id, shape). Recompiles ONLY when shape changes (NOT when data changes). """ def __init__(self): self.ctx = moderngl.create_standalone_context() self.kernel_table = {} # (kernel_id, N) -> PersistentKernel self.invalidation_rules = { 'matmul': lambda N: ('N', N), } self.hits = 0 self.misses = 0 def get(self, kernel_id: str, N: int) -> PersistentKernel: """Get a compiled kernel, compiling it if absent.""" key = (kernel_id, N) if key in self.kernel_table: self.hits += 1 return self.kernel_table[key] self.misses += 1 if kernel_id == 'matmul': k = PersistentKernel(self.ctx, N=N) else: raise KeyError(f"Unknown kernel_id: {kernel_id}") self.kernel_table[key] = k return k def report(self): return f"KernelDB: {self.hits} hits, {self.misses} compiles" def evict(self, kernel_id: str = None): """Drop kernels for memory recovery.""" if kernel_id is None: self.kernel_table.clear() else: self.kernel_table = {k: v for k, v in self.kernel_table.items() if k[0] != kernel_id} # Usage: db = KernelDatabase() # First call to N=16 matmul: compile mm16 = db.get('matmul', 16) A = np.random.randn(16, 16).astype(np.float32) B = np.random.randn(16, 16).astype(np.float32) C1 = mm16(A, B) # compile + dispatch # Reuse cached: skip compile mm16_cached = db.get('matmul', 16) # cache hit C2 = mm16_cached(A, B) # just dispatch # N=32: NEW compile required mm32 = db.get('matmul', 32) # compile + dispatch print(db.report()) # KernelDB: 1 hits, 2 compiles ``` **Cache invalidation rule:** recompile **only when**: - Kernel source changes - Shape parameter changes (e.g., N value) - Buffer binding topology changes (very rare) Never recompile because data changed. --- ## 4. The CUDA Graph Analogy CUDA has a built-in equivalent: **stream capture + graph instantiation:** ```cuda // Capture once cudaStream_t stream; cudaGraph_t graph; cudaGraphExec_t instance; cudaStreamBeginCapture(stream, cudaStreamCaptureModeGlobal); my_kernel<<>>(A_dev, B_dev, C_dev, N); // any kernel sequence cudaStreamEndCapture(stream, &graph); cudaGraphInstantiate(&instance, graph, NULL, NULL, 0); // Re-instantiate any number of times with new data cudaGraphLaunch(instance, stream); // almost zero CPU overhead cudaGraphLaunch(instance, stream); // same kernel, different data ``` The data pointers (`A_dev`, `B_dev`, `C_dev`) are bound to the graph instance. To swap data: - Same buffers, different content: just write into `A_dev`, `B_dev` - Different buffers: graph node update via `cudaGraphExecKernelNodeSetParams` **VGPU / GLSL does the same thing conceptually** but the OpenGL driver handles it automatically — every `program.run()` after the first is essentially equivalent to `cudaGraphLaunch`. --- ## 5. Numba Alternative (for CPU-Fallback GLSL) If you don't have a real GPU but want "compiled once, dispatch many" semantics on CPU: ```python from numba import njit, prange @njit(cache=True, parallel=True, fastmath=True) def matmul_kernel(A, B, C, N): """Compiled once, cached on disk, JIT-evaded on subsequent imports.""" for i in prange(N): for j in range(N): s = 0.0 for k in range(N): s += A[i, k] * B[k, j] C[i, j] = s # First call: JIT compile (~3 seconds) A = np.random.randn(16, 16) B = np.random.randn(16, 16) C = np.zeros_like(A) matmul_kernel(A, B, C, 16) # Subsequent calls in same Python session: zero compile cost matmul_kernel(A, B, C, 16) # fast # Subsequent program calls: still fast — `cache=True` saves .nbi/.nbc files # And the next run will load compiled artifact instantly: # ~/__pycache__/_nbc_python_xxx.nbc ``` Cache key includes signature `(argtypes, parallel, fastmath)`. If you change parameters, you recompile. For PyTorch, similar story with `torch.compile(model, mode="reduce-overhead")` plus `torch.cuda.graphs.CUDAGraph`. --- ## 6. WebGL Constraints & Workarounds WebGL 2 (browser-side) has stricter limits than desktop OpenGL: | Feature | Desktop GL 4.3+ | WebGL 2 | |---|---|---| | Compute shaders | ✓ | ✓ (with WEBGL2) | | SSBO | ✓ | Limited (no atomics, smaller max) | | Uniform updates | Cheap | Cheap (faster than SSBO writes) | | Dispatch rate | ~10⁶/sec | ~10⁵/sec | **WebGL2 strategy: small-N uniforms, large-N SSBO** ```python class WebGLMatMul: def __init__(self, ctx, N=16): # Same as PersistentKernel; works on WebGL2 with fallback ... ``` For very large matrices, WebGL falls back to texture-based compute or JS-side matrices (effectively no GPU). For N≤256, the moderngl pattern transfers to WebGL directly. --- ## 7. Vertically Stratified Kernel Cache A practical pattern for production use: ```python class VGPUComputingBackend: """ Three-tier cache: L1: In-RAM (PersistentKernel already constructed) L2: On-disk numba cache (when CPU-only path) L3: Precompiled shader files (*.spv SPIR-V for Vulkan) Picks the fastest path on first call. """ def __init__(self): self.L1_persistent = {} # RAM cache self.L2_numba_cache = {} # on-disk numba artifacts self.L3_spirv_files = {} # precompiled for Vulkan def matmul(self, A, B, device='auto'): N = A.shape[0] if device == 'auto': device = self._detect_best_backend() if device == 'gpu': return self._gpu_matmul(A, B, N) elif device == 'cpu': return self._cpu_matmul(A, B, N) elif device == 'vulkan': return self._vulkan_matmul(A, B, N) def _gpu_matmul(self, A, B, N): key = ('matmul', N) if key not in self.L1_persistent: ctx = moderngl.create_standalone_context() self.L1_persistent[key] = PersistentKernel(ctx, N) return self.L1_persistent[key](A, B) def _cpu_matmul(self, A, B, N): from numba import njit, prange @njit(cache=True, parallel=True) def _matmul(A, B, C, N): for i in prange(N): for j in range(N): s = 0.0 for k in range(N): s += A[i, k] * B[k, j] C[i, j] = s C = np.zeros_like(A) _matmul(A.astype(np.float32), B.astype(np.float32), C, N) return C def _vulkan_matmul(self, A, B, N): # Load precompiled SPIR-V if available # Otherwise, compile GLSL to SPIR-V offline ... def _detect_best_backend(self): try: import moderngl self.ctx_test = moderngl.create_standalone_context() return 'gpu' except Exception: return 'cpu' ``` --- ## 8. Why This Pattern Matters for VGPU In the VGPU framework, this pattern closes the loop: | VGPU Concept | Map to Pattern | |---|---| | Kernel `F(h, t)` with fixed (D, T) | Same field across many `h₀` queries | | "Compile once" | Persistent shader program | | "Feed data" | Update `h₀` via SSBO write → new initial condition | | "Dispatch many" | `program.run(...)` repeatedly | | Conditional collapse | New shader compiled or runtime branch | **For neural-network inference:** ```python nn_kernel = PersistentKernel(ctx, N=model_input_dim) for sample in dataset: # 60,000 iterations output = nn_kernel(sample, model_weights) # each: ~0.5ms # Total: 30 seconds, zero shader compilations ``` vs recompile-each-time: ```python for sample in dataset: # 60,000 iterations program = ctx.compute_shader(...) # EACH: ~100ms to recompile output = ... # else bug # Total: 100 minutes, mostly wasted ``` **100× speedup** just from keeping the program alive. --- ## 9. Performance Characteristics | Operation | Cost (N=16, RTX class GPU) | Notes | |---|---|---| | First compile (init) | 50–300 ms | One-time per (source, N) | | Subsequent init (cache hit) | <1 ms | Just lookup | | SSBO write | 0.01–1 ms | DMA bound | | Uniform update | <0.001 ms | PCIe register write | | Dispatch overhead | 0.005 ms | CPU-GPU sync | | GPU kernel execution | 0.05–0.5 ms | Compute bound | | SSBO readback | 0.01–1 ms | DMA bound | **Total per dispatch:** ~0.1-3 ms. The persistent pattern preserves this entire cheapness; the recompile path adds ~100 ms to each first call. --- ## 10. Putting it Together A minimal but complete production-ready VGPU cache: ```python """ vgpu_cache.py — Production VGPU cache. Drop this in, replace your ad-hoc shader creation with VGPUCache. """ import moderngl import numpy as np from typing import Optional _MATMUL_GLSL = r""" #version 430 core layout(local_size_x=16, local_size_y=16) in; layout(std430, binding=0) readonly buffer A { float _A[]; }; layout(std430, binding=1) readonly buffer B { float _B[]; }; layout(std430, binding=2) writeonly buffer C { float _C[]; }; uniform int _N; void main() { uint r = gl_GlobalInvocationID.y; uint c = gl_GlobalInvocationID.x; if (r >= uint(_N) || c >= uint(_N)) return; float s = 0.0; for (uint k = 0; k < uint(_N); k++) s += _A[r * uint(_N) + k] * _B[k * uint(_N) + c]; _C[r * uint(_N) + c] = s; } """ class VGPUCache: """Single-instance VGPU kernel cache, lazy-compiling on first use.""" def __init__(self): try: self.ctx = moderngl.create_standalone_context() self.gpu_available = True except Exception: self.ctx = None self.gpu_available = False self._kernels = {} # (op, N) -> PersistentKernel self._stats = {'compiles': 0, 'dispatches': 0, 'cache_hits': 0} def matmul(self, A, B): """Drop-in replacement for np.matmul / A @ B.""" N = A.shape[0] if self.gpu_available: key = ('matmul', N) if key not in self._kernels: self._kernels[key] = PersistentKernel(self.ctx, N) self._stats['compiles'] += 1 else: self._stats['cache_hits'] += 1 result = self._kernels[key](A.astype(np.float32), B.astype(np.float32)) self._stats['dispatches'] += 1 return result else: return A.astype(np.float32) @ B.astype(np.float32) def report(self): return (f"VGPUCache: {self._stats['compiles']} compiles, " f"{self._stats['cache_hits']} cache hits, " f"{self._stats['dispatches']} dispatches") # ==================================================================== # PRODUCTION USAGE # ==================================================================== if __name__ == "__main__": cache = VGPUCache() N = 16 num_calls = 100 print("Running matmul 100x with cache...") for i in range(num_calls): A = np.random.randn(N, N).astype(np.float32) B = np.random.randn(N, N).astype(np.float32) C = cache.matmul(A, B) if i == 0 or i == num_calls - 1: print(f" Call {i+1}: C[0,0]={C[0,0]:.4f}") print(cache.report()) ``` Expected output (with GPU): ``` Running matmul 100x with cache... Call 1: C[0,0]=... Call 100: C[0,0]=... VGPUCache: 1 compiles, 99 cache hits, 100 dispatches ``` --- ## 11. Summary | Anti-pattern | Correct pattern | |---|---| | `compute_shader(source)` per call | One persistent `compute_shader()` per `(source, N)` | | Buffer allocation per call | Pre-allocated, reused buffers | | Buffer binding per call | Bind once after compile | | Compiling when data changes | Compiling only when source/shape changes | | CPU↔GPU sync per launch | Batched dispatches between sync points | The VGPU specification's "dispatch many trajectories" is implemented as "compile once, bind once, dispatch many" on real GPU. The math is unchanged; the speedup is 10-100×. **One line of truth:** keep your compiled shader as a Python attribute, and call `program.run(...)` directly when you want to crunch new data — don't recompile. --- *— End of VGPU Persistent Compute Pattern v1.0 —* ``` A few practical notes on choosing your library: - **moderngl** — what the artifact uses. Cross-platform (Linux/Windows/macOS), works on any GPU with OpenGL 4.3+. Cleanest API for this pattern. - **PyOpenGL** — works but verbose (manual wrapper setup, no high-level buffer objects). - **PyGLet / GLFW** + raw GL — fastest but more boilerplate. - **CUDA via PyCUDA / CuPy** — direct CUDA equivalent, but requires NVIDIA toolkit. - **Numba `@njit(cache=True)`** — if no GPU; saves compiled JIT to disk, loads on next run. For your `Emulate_CUDA` directory, the natural next step is `persistent_kernel.py` plus `vgpu_cache.py` from this artifact. Both fit beside `ref2.py` / `ref3.py` and the GLSL shader from AI_VOXEL Stage 6 follows the same pattern. The CUDA Graph wrapper from §4 maps directly to PyCUDA's `cudaGraph*` API if you want the analogous native-CUDA version without moderngl.