Re Rust gotchas: the standard array container breaks with lengths > 32 - alright, it doesn't "break", but half of the traits, including really basic stuff - more basic than "add" - aren't implemented. So if you write something at small scale, then scale it up slightly - gotcha! And you can't easily reimplement them, either. Basically, the easiest way I've found to do it that way is to fork rustc.
That leaves using dynamic allocation, which I wanted to avoid, or external libraries (nalgebra), which aren't necessarily mature yet. And since type parameters are sometimes namespace separated, and sometimes not, and type specification is sometimes necessary, and sometimes not - examples work, but any tiny modification of them does not, because the example relied on type inference, and what looks like a type, and is declared as a type, is actually a generic.
It also, incidentally, compiled to ninety-three megabytes. To multiply two 4096x4096 f64's.
My point was that in comparing numpy to yodelshady's FORTRAN code, we're probably just comparing two different FORTRAN implementations, one much better than the other (perhaps using a different algorithm, perhaps just better optimised).
(I am, obviously, taking inspiration from "There's plenty of room at the top": https://news.ycombinator.com/item?id=23442123. Naive FORTRAN's slow even compared to C there, but it's not the same machine anymore so not apples to apples. Also, as a small point - I initialised all matrices with random numbers. Below half a second and that might start to be a significant time demand.)
I commented on this thread, because it definitely seems to be in the realm where the compiler matters. So there is opportunity.
┌────────────────┬─────────────────┬──────┐
│ Language │ Detail │ Time │
├────────────────┼─────────────────┼──────┤
│ Python3 and C │ numpy 1.18.4 │ 1.3s │
│ FORTRAN │ gfortran 9.3.0 │ 6.0s │
│ Rust │ 1.43.1, debug │ >60s │
│ Rust │ 1.43.1, release │ 4.0s │
└────────────────┴─────────────────┴──────┘I'd assumed numpy used FORTRAN, but you're right, its heavy-lifting is done in C.
> every single way I tried had "gotchas" and utterly astonishing behaviour in it
As someone who doesn't know much about Rust (or FORTRAN for that matter), what kinds of gotchas? I'd expect numeric code like this to be quite straightforward, as it presumably doesn't lean heavily on memory-management cleverness in the language, or anything like that.
[0] https://en.wikipedia.org/wiki/Strassen_algorithm
edit For clarity, I made a table of your data. Seems surprising that Rust outperformed FORTRAN.