[hansvm]

> In what world would anyone implement their own matrix operations?

In any world where their matrices are "special" -- big, small, scattered across the internet, high read/write ratio, informationally independent columns with 0.01% error tolerable, poorly conditioned, blockwise sparse, used as a sub-component in differentiating anything hairy (like a function of eigenvalues), ....

You can often get halfway decent results stringing together existing matrix primitives, but throw in a 50 clock cycle latency budget, a $5k per-job compute budget, or all sorts of other constraints and you're just as likely to find >>10x performance that you left lying on the table (assuming your matrices are at all special, which is common but not at all guaranteed in arbitrary domains).