4 ms·

Modern fortran is surpisingly pleasant to write. The last time optimisation came up, I tried the author's problem - multiplying two 4096x4096 matrices - in num

Modern fortran is surpisingly pleasant to write. The last time optimisation came up, I tried the author's problem - multiplying two 4096x4096 matrices - in numpy, FORTRAN and Rust, on a standard laptop. Supposedly it took 9 hours in naive Python (I didn't try). Results were: Python3 (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 What surprised me is that python and fortran solutions were written in minutes. The Rust solution took hours and multiple forum posts for help. There's no obvious way, and every single way I tried had "gotchas" and utterly astonishing behaviour in it, particularly with simple, statically-allocated arrays. Fortran could do with proper namespaces, and properly dropping some of its older conventions, though. Reading it feels like archaeology.

MaxBarraclough 6y ago
So to be clear, numpy gives you a highly optimised matrix multiplication algorithm, and your FORTRAN solution was just something you threw together yourself, right? Little surprise that numpy roundly outperforms it, especially if it's using a different algorithm, such as Strassen rather than naive matrix multiplication. [0] > 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 https://en.wikipedia.org/wiki/Strassen_algorithm edit For clarity, I made a table of your data. Seems surprising that Rust outperformed FORTRAN. ┌──────────┬─────────────────┬──────┐ │ Language │ Detail │ Time │ ├──────────┼─────────────────┼──────┤ │ Python3 │ 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 │ └──────────┴─────────────────┴──────┘
- Someone 6y ago
  The main reason for numpy’s speed is that it uses highly optimized code written in another language (for matrix multiplications, quite possibly FORTRAN).
  - MaxBarraclough 6y ago
    Yes, of course numpy's heavy-lifting code isn't written in Python. 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).
- yodelshady 6y ago
  Oh no, simpler than that even. The FORTRAN solution was also using somebody else's algorithm (matmul), which is intrinsic, and, I believe, not multithreaded. I'm not employed to develop matrix multiplication algorithms, but it will benefit me to know which tools work best. 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.
- yodelshady 6y ago
  It does seem surprising (and I was all ready to slate rust until I remembered to build for release). If you care that much - naive FORTRAN iteration on the same machine, same compiler, took 1200 s. (I am, obviously, taking inspiration from "There's plenty of room at the top": https://news.ycombinator.com/item?id=23442123 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.
- pjmlp 6y ago
  Actually the table is more like ┌────────────────┬─────────────────┬──────┐ │ 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 │ └────────────────┴─────────────────┴──────┘
  - MaxBarraclough 6y ago
    I could have put 'host language', but I figured it was clear from context. I'd assumed numpy used FORTRAN, but you're right, its heavy-lifting is done in C.
C1sc0cat 6y ago
These days you don't use the older Fortran'isams Arithmetic IF's for example and I suspect no one uses BLOCK DATA any more. But Backwards compatibly is important here as there is a lot of legacy code still in use - yes I know arithmetic if has finally been removed.
deleted 6y ago
[deleted]
_n_b_ 6y ago
I used to essentially be a professional Fortran programmer. It's been 7-8 years since I've done that, so this was a nice trip down memory lane. Anyway. > Python3 (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 suspect that you could get the Fortran and numpy results to converge by turning on the compiler option to use the local BLAS library for MATMUL, which will normally use an OpenMP-threaded solution. Upfront note: these days "Fortran" is preferred to "FORTRAN." > Fortran could do with proper namespaces, and properly dropping some of its older conventions, though. Reading it feels like archaeology. The newer Fortran standards support pretty decent namespaces (USE module x, ONLY : y_ and a limited OO with classes and class methods. The problem, at least when I was doing Fortran work everyday, was that the compilers did not fully support the newer Fortran standards. We found the Intel compiler to be the best of the bunch, but we also found a non-trivial number of compiler bugs and standard options that weren't supported. The Intel Math Kernel Library offers pretty good linear algebra performance, too. Re: older conventions, we had a lot of modern Fortran code that was easy to read and avoided all of the old Fortran 77 nonsense that really confusing to read if you weren't from that era (e.g., the infamous arithmetic GOTO). However, it was nice being able to still incorporate "just works" code from the 80s that did complicated calculational things that nobody really wanted to touch too much. You know, the kind of subroutine that starts with the comment, "This code _seems to_ ..." I suspect we're not the only organization that wanted to keep that code around but be able to add modern Fortran around and on top of it.
deleted 6y ago
[deleted]

MaxBarraclough 6y ago

So to be clear, numpy gives you a highly optimised matrix multiplication algorithm, and your FORTRAN solution was just something you threw together yourself, right? Little surprise that numpy roundly outperforms it, especially if it's using a different algorithm, such as Strassen rather than naive matrix multiplication. [0] > 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 https://en.wikipedia.org/wiki/Strassen_algorithm edit For clarity, I made a table of your data. Seems surprising that Rust outperformed FORTRAN. ┌──────────┬─────────────────┬──────┐ │ Language │ Detail │ Time │ ├──────────┼─────────────────┼──────┤ │ Python3 │ 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 │ └──────────┴─────────────────┴──────┘

Someone 6y ago
The main reason for numpy’s speed is that it uses highly optimized code written in another language (for matrix multiplications, quite possibly FORTRAN).
- MaxBarraclough 6y ago
  Yes, of course numpy's heavy-lifting code isn't written in Python. 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).
yodelshady 6y ago
Oh no, simpler than that even. The FORTRAN solution was also using somebody else's algorithm (matmul), which is intrinsic, and, I believe, not multithreaded. I'm not employed to develop matrix multiplication algorithms, but it will benefit me to know which tools work best. 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.
yodelshady 6y ago
It does seem surprising (and I was all ready to slate rust until I remembered to build for release). If you care that much - naive FORTRAN iteration on the same machine, same compiler, took 1200 s. (I am, obviously, taking inspiration from "There's plenty of room at the top": https://news.ycombinator.com/item?id=23442123 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.
pjmlp 6y ago
Actually the table is more like ┌────────────────┬─────────────────┬──────┐ │ 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 │ └────────────────┴─────────────────┴──────┘
- MaxBarraclough 6y ago
  I could have put 'host language', but I figured it was clear from context. I'd assumed numpy used FORTRAN, but you're right, its heavy-lifting is done in C.

C1sc0cat 6y ago

These days you don't use the older Fortran'isams Arithmetic IF's for example and I suspect no one uses BLOCK DATA any more. But Backwards compatibly is important here as there is a lot of legacy code still in use - yes I know arithmetic if has finally been removed.

deleted 6y ago

[deleted]

_n_b_ 6y ago

I used to essentially be a professional Fortran programmer. It's been 7-8 years since I've done that, so this was a nice trip down memory lane. Anyway. > Python3 (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 suspect that you could get the Fortran and numpy results to converge by turning on the compiler option to use the local BLAS library for MATMUL, which will normally use an OpenMP-threaded solution. Upfront note: these days "Fortran" is preferred to "FORTRAN." > Fortran could do with proper namespaces, and properly dropping some of its older conventions, though. Reading it feels like archaeology. The newer Fortran standards support pretty decent namespaces (USE module x, ONLY : y_ and a limited OO with classes and class methods. The problem, at least when I was doing Fortran work everyday, was that the compilers did not fully support the newer Fortran standards. We found the Intel compiler to be the best of the bunch, but we also found a non-trivial number of compiler bugs and standard options that weren't supported. The Intel Math Kernel Library offers pretty good linear algebra performance, too. Re: older conventions, we had a lot of modern Fortran code that was easy to read and avoided all of the old Fortran 77 nonsense that really confusing to read if you weren't from that era (e.g., the infamous arithmetic GOTO). However, it was nice being able to still incorporate "just works" code from the 80s that did complicated calculational things that nobody really wanted to touch too much. You know, the kind of subroutine that starts with the comment, "This code _seems to_ ..." I suspect we're not the only organization that wanted to keep that code around but be able to add modern Fortran around and on top of it.

deleted 6y ago

[deleted]