This page was generated from unit-5.6-gpu/matrixfree.ipynb.
5.6.2 Matrix-free operator application on the GPU¶
For operator application we do not need the assembled sparse matrix: the bilinear form
factors into the element-wise differential operator \(B\) (universal on the reference element), a point-wise operation \(D\) built from the coefficient and the geometry, and the transposed test-function operator \(B^T\). The mf option of the BilinearForm sets up this factorization instead of assembling a matrix, and CreateDeviceMatrix compiles one fused GPU kernel for the whole operator - gather, \(B\), \(D\), \(B^T\), scatter.
The kernel is written in the common GPU language of unit 5.6.1, so the same code runs on CUDA, Metal and the host reference backend. Block sizes, warps and the launch size are set by MFOpts, with defaults that suit both NVIDIA and Apple gpus (the group count adapts to the device’s compute units).
[1]:
from ngsolve import *
from ngsolve.comp import MFOpts
from ngsolve.gpu import backend # registers cuda, metal or the host reference device
from time import time
print ("gpu backend:", backend)
gpu backend: host
[2]:
with TaskManager():
mesh = Mesh(unit_cube.GenerateMesh(maxh=0.3))
for l in range(1): # increase for serious timings
mesh.Refine()
mesh.ngmesh.OrderElements()
mesh = Mesh(mesh.ngmesh)
print ("elements:", mesh.ne)
elements: 1648
Correctness and timing¶
For each space and order we assemble the operator twice: matrix-free (mf=MFOpts()) and as a sparse matrix. The matrix-free device operator is applied to a random vector and compared against the sparse product. Then we time
the sparse matrix on the host (with
TaskManager),the sparse matrix on the device (cuda only),
the matrix-free device operator.
Throughput is reported in GDofs/sec, counting the dofs of the discontinuous version of the space (the actual data volume the operator moves).
[3]:
def Sync (vec):
try: vec.WaitUntilCompleted() # metal
except AttributeError: s = vec.Norm() # any backend: forces completion
def TimeApply (apply_once, syncvec, runs):
for j in range(min(runs,20)): apply_once()
Sync(syncvec)
ts = time()
for j in range(runs): apply_once()
Sync(syncvec)
return (time()-ts)/runs
[4]:
def Bench (space, order, form, runs=100):
fes = space(mesh, order=order)
u,v = fes.TnT()
equ = form(u,v)
ndof = Discontinuous(fes).ndof
bfmf = BilinearForm(equ, mf=MFOpts()).Assemble()
gpumat = bfmf.mat.CreateDeviceMatrix()
bfsp = BilinearForm(equ).Assemble()
# correctness
gfu = GridFunction(fes)
gfu.vec.SetRandom()
vy = (bfsp.mat * gfu.vec).Evaluate()
xdev = gpumat.CreateRowVector(); xdev.data = gfu.vec
ydev = gpumat.CreateColVector(); ydev.data = gpumat * xdev
yh = bfsp.mat.CreateColVector(); yh.data = ydev
relerr = Norm(vy-yh)/Norm(vy)
# timings
x = gfu.vec.CreateVector(); x.SetRandom()
y = gfu.vec.CreateVector()
def apply_host(): y.data = bfsp.mat * x
with TaskManager():
t_host = TimeApply(apply_host, y, max(runs//5,5))
t_spdev = None
if backend == "cuda":
spdev = bfsp.mat.CreateDeviceMatrix()
xs = x.CreateDeviceVector(copy=True)
ys = (spdev*xs).Evaluate()
def apply_spdev(): ys.data = spdev * xs
t_spdev = TimeApply(apply_spdev, ys, runs)
def apply_mf(): ydev.data = gpumat * xdev
t_mf = TimeApply(apply_mf, ydev, runs)
G = lambda t: ndof/t*1e-9 if t else None
return ndof, relerr, G(t_host), G(t_spdev), G(t_mf)
[5]:
def RunTable (space, form, orders=range(1,6)):
print (f"{space.__name__}: order ndof relerr sparse-host sparse-dev MF-dev [GDofs/s]")
for order in orders:
ndof, relerr, gh, gs, gm = Bench(space, order, form)
gs = f"{gs:9.2f}" if gs else " -"
print (f" {order} {ndof:9d} {relerr:.1e} {gh:9.2f} {gs} {gm:8.2f}")
RunTable (H1, lambda u,v: grad(u)*grad(v)*dx)
H1: order ndof relerr sparse-host sparse-dev MF-dev [GDofs/s]
1 6592 4.5e-16 0.49 - 0.02
2 16480 5.2e-16 0.36 - 0.01
3 32960 6.5e-16 0.19 - 0.01
4 57680 1.1e-15 0.09 - 0.01
5 92288 3.3e-15 0.05 - 0.01
[6]:
RunTable (L2, lambda u,v: u*v*dx)
L2: order ndof relerr sparse-host sparse-dev MF-dev [GDofs/s]
1 6592 1.7e-16 0.25 - 0.01
2 16480 2.4e-16 0.22 - 0.01
3 32960 3.2e-16 0.13 - 0.01
4 57680 4.1e-16 0.06 - 0.01
5 92288 5.3e-16 0.04 - 0.00
Measured results¶
Throughput in GDofs/sec on a 3-times refined unit cube (approx. 10^5 elements), 2026-09-01.
Apple M4 Pro (Metal, matrix-free in fp32):
space |
order |
sparse host |
MF device |
|---|---|---|---|
H1 |
1 |
15.2 |
2.7 |
H1 |
2 |
3.2 |
7.4 |
H1 |
3 |
1.2 |
7.8 |
H1 |
4 |
0.58 |
7.8 |
H1 |
5 |
0.29 |
3.9 |
L2 |
1 |
3.4 |
7.0 |
L2 |
2 |
1.4 |
5.9 |
L2 |
3 |
0.75 |
8.3 |
L2 |
4 |
0.44 |
7.0 |
L2 |
5 |
0.30 |
2.2 |
NVIDIA RTX 5090 (CUDA):
space |
order |
sparse host |
sparse device |
MF device fp64 |
MF device fp32 |
|---|---|---|---|---|---|
H1 |
1 |
4.3 |
16.7 |
11.0 |
26.1 |
H1 |
2 |
0.58 |
25.0 |
11.2 |
37.0 |
H1 |
3 |
0.18 |
11.0 |
5.8 |
29.6 |
H1 |
4 |
0.13 |
5.7 |
4.4 |
22.8 |
H1 |
5 |
0.07 |
3.4 |
2.7 |
15.2 |
L2 |
1 |
1.4 |
17.0 |
15.1 |
30.9 |
L2 |
2 |
0.48 |
10.3 |
13.1 |
35.8 |
L2 |
3 |
0.26 |
6.4 |
12.5 |
42.4 |
L2 |
4 |
0.15 |
3.9 |
7.4 |
28.5 |
L2 |
5 |
0.10 |
2.3 |
2.7 |
9.4 |
Observations:
Beyond the lowest order the matrix-free operator clearly beats the assembled matrix - it moves element data instead of the sparse matrix, whose size grows with the square of the local dofs.
On consumer NVIDIA hardware fp64 arithmetic is slow (1/64 of fp32), so there the matrix-free operator wants
MFOpts(fp32=True); data-center gpus do not have this restriction. On Metal fp32 is the only choice.The drop at the highest orders is a shared-memory/occupancy limit of the fused kernel - the block sizes in
MFOptsare the tuning knobs.Dof numbering matters for the gather/scatter:
Reorder(fes)(Morton element order with first-touch dof numbering) gives another 5-15% on the gpu.