Het is een bijverschijnsel van ouder worden dat je vaker terugdenkt aan vroeger. Dus toen ik recent een nieuwe rekencomputer aan het samenstellen was, moest ik terugdenken aan mijn eerste PC. Deze heb ik in 1993 gekocht, een 386DX-40 met 4 MB RAM en een harde schijf van 170 MB – een enorme vooruitgang ten opzichte van de Commodore 64 waarmee ik het daarvoor moest toen. En toch was deze computer, ondanks een tussentijdse update naar een 486DX-50, na zes jaar verschrikkelijk achterhaald. Ik heb hem daarom in 1999 vervangen door een Pentium III 450 met 128 MB RAM en een harde schijf van 10 GB, de computer die mij in mijn studie begeleid heeft. Qua besturingssysteem ging het van MS-DOS 6.22 met Windows 3.11 naar Windows 98 SE en SuSE Linux.
Wat een verschil met nu – want mijn oude rekencomputer, nu ook al meer dan zes jaar oud, doet ook na velen Terabytes aan verwerkte puntenwolken nog steeds trouw dienst. Hij is geschikt voor Windows 11 en het is voornamelijk de grafische kaart die aan een upgrade toe is – is er dus de laatste jaren werkelijk zo weinig vooruitgang op het gebied van rekenkracht geweest?
Om dat te kunnen beoordelen moet je het meten en in de historische context plaatsen. Voor het meten heb je benchmarks – populair zijn bijv. Cinebench en Geekbench. Maar die draaien niet op oude hardware, waardoor een vergelijking over vele generaties computers heen niet mogelijk is.
In de wereld van supercomputers en wetenschappelijke berekeningen is vooral belangrijk hoe snel een computer numerieke berekeningen met zwevendekommagetallen kan uitvoeren. Deze rekensnelheid wordt uitgedrukt in FLOPS, wat staat voor floating point operations per second. Voor supercomputers wordt dit gemeten met de LINPACK benchmark, die lineaire vergelijkingsstelsels opbouwt en oplost – zoals wij geodeten het uit de vereffeningsrekening kennen. LINPACK meet daarbij de snelheid bij het verwerken van 64 bit getallen, dit staat bekend als double precision – 32 bit is single precision.
Het leuke is dat voor de FLOPS en LINPACK resultaten veel historische data beschikbaar is die een vergelijking door de jaren heen mogelijk maakt. Voor PCs kun je dat eigenlijk pas vanaf de 486 processor doen, want dat was de eerste met een ingebouwde rekeneenheid voor zwevendekommagetallen (floating point unit of FPU). De 386 en eerdere x86 CPUs moesten deze in software nabootsen, wat echt heel traag was – tenzij je een dure coprocessor kocht, voor de 386 was dat de 387.
Laten we daarom eerst naar de supercomputers kijken om zo verder terug te kunnen gaan in de tijd. De eerste computers zoals de Zuse Z3 en ENIAC waren ontwikkeld voor wetenschappelijke berekeningen. De Z3 deed krap een seconde over het optellen van twee getallen, ENIAC was al veel sneller met een paar duizend berekeningen per seconde. Maar het onderscheid tussen supercomputers en “gewone” computers werd eigenlijk pas vanaf de jaren 60 gemaakt, toen het gebruik van computers (en dan hebben we het over mainframes, dus computers die hele kamers vulden) enorm toenam, met de IBM System/360 als de meest bekende computerfamilie uit die tijd.
In 1964 bracht Control Data Corporation zijn CDC 6600 uit die vaak als de eerste supercomputer wordt beschouwd. Deze behaalde 3 MFLOPS, dus drie miljoen berekeningen per seconde.
De CDC 6600 was ontworpen door ene Seymour Cray. Nu gaat er misschien een belletje rinkelen. Want in 1972 begon hij zijn eigen bedrijf, Cray Research. Sindsdien is de naam Cray bijna een synoniem voor supercomputers.
De Cray X-MP was een succesvol model en van 1983 tot 1985 de snelste computer ter wereld met 800 MFLOPS snelheid – dus bijna 300 keer zoveel als de CDC 6600. Daarbij was de X-MP een relatief compact systeem dat bijna als een meubel leek. De opvolger Y-MP kraakte de GFLOPS barriere met meer dan 2 miljard berekeningen per seconde. Typische kenmerken voor de Crays uit deze tijd waren het gebruik van vectorprocessoren (die met niet met enkele getallen rekenen, maar met meerdere tegelijk) en vloeistofkoeling.
Ook NEC uit Japan maakte met zijn SX-reeks gebruik van vectorprocessoren. Maar in de jaren 90 kwam een nieuwe trend: niet meer het gebruik van maatwerk vectorprocessoren, maar in plaats daarvan grote aantallen van veel goedkopere “gewone” CPUs. Cray speelde op deze trend in met de T3E die gebruik maakte van de Alpha processoren van DEC. De snelste incarnatie van de T3E haalde 891 GFLOPS en was daarmee 1000 keer sneller dan de X-MP van 15 jaar eerder. Sinds 1993 wordt twee keer per jaar een lijst van de snelste supercomputers ter wereld volgens de LINPACK gepubliceerd, de TOP500 lijst. In juni 2000 stond deze snelste versie van de Cray T3E op de 7e plek.
Mijn eerste echte contact met supercomputers was tijdens een door mij georganiseerd bezoek voor geodesiestudenten aan het Höchstleistungsrechenzentrum Stuttgart (HLRS) in 2003. Het HLRS had toen zowel een NEC SX-5 als een Cray T3E. De laatste stond bij Daimler, want daar moest per vierkante meter vloeroppervlak betaald worden, wat dus handig was met de compacte vloeistofgekoelde T3E. En waarom twee systemen? Massief parallelle system moeten anders geprogrammeerd worden dan vectorsystemen, waardoor sommige applicaties zich beter voor de ene dan de andere architectuur lenen.
Inmiddels was ook een derde architectuur populair geworden: clusters, opgebouwd uit gewone servers, maar wel vaak met een snellere netwerkverbinding. Zonder het complexe interne netwerk van een “echte” supercomputer waren deze een stuk goedkoper, maar wederom niet geschikt voor elke applicatie.
Voor mijn afstudeeropdracht wilde ik heel graag iets met supercomputers doen – in de geodesie beland je dan bijna automatisch bij de zwaartekrachtmodellering. En zo mocht ik in 2004 vijf maanden lang algoritmes voor het berekenen van lineaire vergelijkingsstelsels parallelliseren en op de systemen van het HLRS testen. Die hadden net bij NEC een nieuw systeem besteld en als tijdelijke oplossing een SX-6 gekregen waarop wij konden rekenen. Maar ja, er was uiteraard een queue, een wachtrij. Dus gingen we op het TX-7 frontend systeem rekenen dat met 16 Itanium II CPUs ook aardig vlot was…
NEC was in deze tijd sowieso een grote naam, want zij hadden in 2002 de Earth Simulator geleverd. Met een snelheid van 36 TFLOPS was deze vijf keer zo snel als de nummer 2 op de TOP500 lijst en meer dan twee jaar lang de snelste computer ter wereld. De TFLOPS barrière was trouwens al in 1997 doorbroken, door de massief parallelle ASCI Red met Intel Pentium Pro CPUs.
Na mijn afstuderen in Stuttgart ging ik bij de TU Delft aan de slag. Wij konden toen op de nationale supercomputers TERAS en ASTER rekenen. Veel leuker was echter dat wij onze eigen cluster konden aanschaffen. Dit door mij CLEOPATRA gedoopte systeem bestond uit 33 servers en haalde een snelheid van ruim 600 GFLOPS. Vijf jaar eerder was dat nog voldoende geweest voor een plek in de bovenste 10% van de TOP500 lijst.
Tegen het einde van mijn promotieonderzoek was CLEOPATRA te klein geworden voor mij, ik had vooral meer geheugen nodig. Ik paste een filtertechniek toe waarover een luisteraar bij een presentatie opmerkte “mij was verteld dat dit niet mogelijk was”. Wel dus – als je er genoeg rekenkracht tegenaan gooide. Die kwam in 2009 van de nieuwe nationale supercomputer Huygens 2 met 60 TFLOPS performance – 100 keer zo snel als onze CLEOPATRA.
Nu, ruim 15 jaar later, is 60 TFLOPS niet zoveel meer. Het “traagste” systeem op de TOP500 lijst van juni 2025 haalt 2,4 PFLOPS. Maar we zijn zelfs al in het “Exascale” tijdperk aangekomen – de snelste computer van dit moment, El Capitan, haalt 1,7 EFLOPS – en is gebouwd door Cray, nu een dochter van HP. Supercomputers zijn inmiddels ook modulair opgebouwd en bestand altijd uit verschillende gekoppelde systemen, niet meer uit een enkel eenheid zoals tot in de jaren ’90.
Maar terug naar de computers voor gewone mensen. Mijn 486DX-50 zal ongeveer 2MFLOPS behaald hebben – iets minder dus dan de CDC 6600 van 30 jaar eerder. De Pentium III 450 haalde meer dan 60 MFLOPS. De toename in snelheid was niet alleen het resultaat van de toename in kloksnelheid, maar ook door het sneller uitvoeren van de berekeningen door verbetering van de FPU. De 800 MFLOPS van de Cray X-MP werden ingehaald door de in 2002 verschenen 3 GHz versie van de Pentium 4 – in slechts drieënhalf jaar tijd was Intel van 450 Mhz naar 3 GHz en van 60 naar 800 MFLOPS gegaan.
De Pentium 4 was geen efficiënt ontwerp, maar bracht aan het einde van zijn loopbaan nog twee belangrijke innovaties die nu gemeengoed zijn: de door AMD ontwikkelde 64-bit extensie en met de Pentium D de eerste dual-core processor, waarbij twee bijna onafhankelijke processoren op één chip zitten. En zo’n Pentium D had ik als eerste bij de faculteit Luchtvaart & Ruimtevaarttechniek van de TU Delft onder mijn bureau staan, waarmee ik een tijdje genoeg rekenkracht voor mijn onderzoek had.
Tot nu toe heb ik het telkens over LINPACK resultaten gehad als benchmark voor het bepalen van de rekensnelheid in FLOPS. Ik wil echter iedereen de kans geven om zelf een benchmark te draaien, maar de volledige LINPACK benchmark is misschien een beetje veel voor de gemiddelde computer en moeilijk te installeren. Een belangrijk onderdeel hiervan is de matrix-matrix-multiplicatie. Deze is goed te optimaliseren waardoor veel computers redelijk in de buurt (en dan bedoel ik 50% of beter) van de theoretische pieksnelheid komen. De bovengenoemde Pentium D kon per core twee berekeningen per cyclus uit voeren, dus in theorie 12 GFLOPS – en ik heb oude benchmark resultaten van mij gevonden met bijna 11 GFLOPS snelheid, dus rond de 90% van de pieksnelheid. Mij is niet helemaal duidelijk waarom de LINPACK resultaten voor de gewone Pentium 4 zoveel trager zijn, misschien omdat het veel minder efficiënte oplossen van het vergelijkingsstelsel hierbij een rol speelt. Die 11 GFLOPS waren in ieder geval 10 jaar daarvoor nog goed geweest voor de 111e plek op de TOP500 lijst.
De afgelopen 20 jaar is de ontwikkeling uiteraard gewoon doorgegaan. Oké, qua kloksnelheid zijn er grenzen en zitten we als turbo voor korte tijd net iets onder de 6 GHz, dus krap het dubbele van de Pentium 4 van ruim 20 jaar geleden. Maar twee dingen zijn wel enorm gegroeid: het aantal cores (een AMD Ryzen 9 9950X heeft 16 identieke cores, een Intel Core Ultra 9 285 zelfs 24 stuks in twee verschillende soorten, performance en efficiency) en het aantal operaties dat per cyclus uitgevoerd kan worden – dat zijn er inmiddels 32.
Zo’n Ryzen 9 9950X komt dus theoretisch boven de 2 TFLOPS uit – toevallig heb ik er eentje onder mijn bureau staan, en met een in Python geschreven benchmarkprogramma meet ik 1,4 TFLOPS – meer dan 100x zoveel als de Pentium D van 20 jaar geleden. Mijn oude rekencomputer met een Intel Core i5-9600K haalt 250 GFLOPS, dus iets meer dan een zesde hiervan.
Maar ho maar, rekent niet iedereen tegenwoordig op grafische kaarten of daarvan afgeleide rekenversnellers, omdat die zo snel zijn? Laten we eens kijken, wat haalt Nvidia’s paradepaardje voor gamers, de RTX 5090? In mijn benchmark precies 10% meer dan de Ryzen 9 9950X, ruim 1,5 TFLOPS – wat een teleurstelling!
De reden? Nvidia’s consumer-GPUs en de daarvan afgeleide workstation GPUs zijn verschrikkelijk traag als het om berekeningen met 64 bit double-precision getallen gaat. Daarom kun je de benchmark ook met 32 bit single-precision getallen berekenen. Dan haalt de Ryzen 9 9950X 3,2 TFLOPS, de RTX 5090 25,6 TFLOPS – 16 keer sneller dan met double precision. Als je ook in double precision snel wilt rekenen, dan vraagt Nvidia veel meer geld voor een van de high-end CPUs. De RTX 2070 in de oude PC haalt in single precision trouwens ook al 5 TFLOPS, dus hier is de vooruitgang nog kleiner.
Dus nee, de snelheid neemt niet meer met dezelfde stappen toe als 25 jaar geleden, en de supercomputers hebben de afstand weer vergroot. Waar de Pentium D dus 10 jaar voor zijn tijd redelijk bovenaan de TOP500 lijst had gestaan geldt dat voor de Ryzen 9 9950X niet – in 2015 had je al 200 TFLOPS nodig om überhaupt in de TOP500 lijst te staan.
Maar is dat erg? We zijn ook op een punt aangekomen waar een computer van een paar jaar oud nog prima geschikt is voor de dagelijkse taken die zich toch vaak in een webbrowser afspelen. Het zijn vooral de gamers die nog steeds meer rekenkracht willen. We zien dus dat naast de pure rekensnelheid ook de efficiëntie steeds meer aandacht krijgt. Een laptop die een hele dag op één acculading meegaat, een paar jaar geleden nog een ijdele droom, is inmiddels werkelijkheid. Zelfs bij de supercomputers is de strijd om meer efficientie begonnen – er is inmiddels ook een GREEN500 lijst, waarbij de GFLOPS per watt tellen.
Bronnen
https://en.wikipedia.org/wiki/Floating_point_operations_per_second
http://www.roylongbottom.org.uk/linpack%20results.htm
Wittwer, T.: An Introduction to Parallel Programming. PDF.
Benchmark code
Deze is door Github Copilot geschreven en door mij uitgebreid met het gebruik van CuPy. Hiervoor heb je een geschikte Nvidia GPU, de CUDA Toolkit en uiteraard CuPy nodig.
import argparse
import time
from typing import Tuple, Callable, Dict
import numpy as np
import cupy as cp
#!/usr/bin/env python3
"""
Matrix-matrix multiplication benchmark.
Usage examples:
python bench_matmul.py --sizes 256 512 --repeats 3 --methods numpy python einsum
python bench_matmul.py --size 1024 --repeats 5
"""
def parse_args():
p = argparse.ArgumentParser(description="Matrix-matrix multiplication benchmark")
p.add_argument("--sizes", "-s", nargs="+", type=int, default=[256, 512, 1024, 8192],
help="Square matrix sizes to benchmark (space separated).")
p.add_argument("--repeats", "-r", type=int, default=3, help="Number of repetitions per test.")
p.add_argument("--dtype", "-d", choices=["float64", "float32"], default="float64")
p.add_argument("--seed", type=int, default=42, help="Random seed")
p.add_argument("--methods", "-m", nargs="+",
choices=["numpy", "einsum", "python", "cupy"],
default=["numpy", "einsum","cupy"],
help="Which implementations to run.")
return p.parse_args()
def make_matrices(n: int, dtype=np.float64, seed: int = 0) -> Tuple[np.ndarray, np.ndarray]:
rng = np.random.default_rng(seed + n)
A = rng.standard_normal((n, n), dtype=dtype)
B = rng.standard_normal((n, n), dtype=dtype)
return A, B
def ops_count(n: int) -> int:
# For square n x n: 2*n^3 floating-point operations (mul + add)
return 2 * (n ** 3)
def bench(func: Callable[[np.ndarray, np.ndarray], np.ndarray],
A: np.ndarray, B: np.ndarray, repeats: int) -> Tuple[float, np.ndarray]:
# warm up
C = func(A, B)
t0 = time.perf_counter()
for _ in range(repeats):
C = func(A, B)
t1 = time.perf_counter()
elapsed = (t1 - t0) / repeats
return elapsed, C
def impl_cupy(A: np.ndarray, B: np.ndarray) -> np.ndarray:
return cp.asnumpy(cp.matmul(cp.asarray(A),cp.asarray(B)))
def impl_numpy(A: np.ndarray, B: np.ndarray) -> np.ndarray:
return np.matmul(A,B)
#return A.dot(B)
def impl_einsum(A: np.ndarray, B: np.ndarray) -> np.ndarray:
return np.einsum("ik,kj->ij", A, B, optimize=True)
def impl_python(A: np.ndarray, B: np.ndarray) -> np.ndarray:
# Pure-Python triple loop on lists. Very slow for large n.
# Convert to nested lists for faster indexing in Python.
a = A.tolist()
b = B.tolist()
n = len(a)
# initialize zero matrix
c = [[0.0] * n for _ in range(n)]
for i in range(n):
ai = a[i]
ci = c[i]
for k in range(n):
aik = ai[k]
bk = b[k]
# unroll inner add across j
for j in range(n):
ci[j] += aik * bk[j]
return np.array(c, dtype=A.dtype)
def main():
args = parse_args()
dtype = np.float64 if args.dtype == "float64" else np.float32
impls: Dict[str, Callable[[np.ndarray, np.ndarray], np.ndarray]] = {}
if "numpy" in args.methods:
impls["numpy"] = impl_numpy
if "einsum" in args.methods:
impls["einsum"] = impl_einsum
if "python" in args.methods:
impls["python"] = impl_python
if "cupy" in args.methods:
impls["cupy"] = impl_cupy
print(f"Matrix-matrix multiplication benchmark")
print(f"Sizes: {args.sizes}, repeats: {args.repeats}, dtype: {args.dtype}, methods: {list(impls.keys())}")
print(f"{'method':10s} {'n':6s} {'time(s)':>10s} {'GFLOPS':>10s} {'max_err':>12s}")
for n in args.sizes:
A, B = make_matrices(n, dtype=dtype, seed=args.seed)
# reference result (numpy) for correctness
C_ref = impl_numpy(A, B)
ops = ops_count(n)
for name, fn in impls.items():
# warn if python impl with large n
if name == "python" and n > 512:
# still run, but it may be extremely slow
pass
try:
elapsed, C = bench(fn, A, B, repeats=args.repeats)
except Exception as e:
print(f"{name:10s} {n:6d} ERROR during run: {e}")
continue
gflops = ops / elapsed / 1e9
max_err = float(np.max(np.abs(C - C_ref)))
print(f"{name:10s} {n:6d} {elapsed:10.6f} {gflops:10.3f} {max_err:12.3e}")
if __name__ == "__main__":
main()
