5 ms·

I actually had a fantastic experience with Fortran lately. I ported a compute kernel from python/numpy to Fortran 2018, partially due to the GIL and partly so I

I actually had a fantastic experience with Fortran lately. I ported a compute kernel from python/numpy to Fortran 2018, partially due to the GIL and partly so I could use Intel's compiler. The performance improvement was tremendous. Several times faster per core, then multiplying further because I could take advantage of threading. In all, the 3 day project increased actual throughput 450x. (I considered JAX, but the code in question was not amenable to a compute graph. Another option was to thread by fork, and use IPC.) I liked the language itself more than expected. You have something like "generics" with tensors. Suppose you pass a parameter, N, and you also would like to pass a tensor, and you would like to specify the tensor's shape (N, N). You can do this; the parameter type constraints can reference other parameters. Tensors and various operations are first-class types, so the compiler can optimise operations easily for the system you're building on. In my case, I got 80% improvement from ifx over gfortran. Invocation from Python was basically the same as a C library. Both Python and Fortran have facilities for C interop, and Numpy can be asked to lay out tensors in a Fortran compatible way. Part of what eased the port was that Numpy seems to be a kind of "Fortran wrapper". The ergonomics on tensor addressing, slicing and views is identical.

prennert 10mo ago
I did something similar many years ago. I was amazed that Fortran was not more discussed as an option to write performant code within a Python / numpy codebase. At the time everyone seems to default to using C instead. But Fortran is so much easier! It even has slicing notations for arrays and the code looked so much like Numpy as you say.
3uruiueijjj 10mo ago
I've never found anything to back this up, but my impression was that both the Python / Numpy and Fortran 90 slicing operations were directly inspired by MATLAB (although most of the ideas go back to at least Algol 68). It also helps that Fortran compatibility is a must for pretty much anything that expects to use BLAS.
- thatjoeoverthr 10mo ago
  I'm learning just now how old MATLAB is. Wow!
foxglacier 10mo ago
Yea, Fortran is nice to use for so-called "scientific" computing. It has high performance as well as some handy intrinsic functions like DOT_PRODUCT and TRANSPOSE but the best features to me are colon array slicing syntax like Python/Numpy and arrays being indexed from 1 which makes converting math equations into code more natural without constantly worrying about off-by-one errors. I wouldn't call multi-dimensional arrays tensors though. That's a bit of a bastardization of the term that seemed to be introduced by ML guys. It wasn't until I started using Fortran that I realized how similar it is to BASIC which must have been a poor-man's Fortran.
- PaulHoule 10mo ago
  e.g. "tensors" are like "vectors" in they transform in a specific way when the coordinate system changes; what felt so magic about vectors as an undergrad was that they embody "the shape of space" and thus simplify calculations. If you didn't have vectors, Maxwell's equations would spill all over the place. Tensors on the other hand are used in places like continuum mechanics and general relativity where something more than vectors are called for but you're living in the same space(/time) with the same symmetries.
  - foxglacier 10mo ago
    Sure, but not all arrays are tensors. In continuum mechanics, you often represent tensors as matrices or vectors to exploit symmetry and then it's really confusing to keep calling them tensors because they don't behave like their tensor equivalent.
  - Joker_vD 10mo ago
    > If you didn't have vectors, Maxwell's equations would spill all over the place. What do you mean, "would": they did! :) The original equations had 20 separate equations, although Maxwell himself tried to reformulate them in quaternions. But if you look e.g. at works of Lorentz, or Einstein's famous 1905 paper, you'll see the fully-expanded version of them. The vector form really didn't fully catch until about the middle of the XX century.
pantsforbirds 10mo ago
If your problem fits into arrays/matrices/vectors as the only required datastructures, Fortran is a VERY good language.
drnick1 10mo ago
Any reason why you couldn't do this in C/C++? I use the Eigen matrix library as a replacement for Numpy when I need more performance and the code looks surprisingly very similar.
- thatjoeoverthr 10mo ago
  Of course I could in C. Intel and Nvidia are both offering both C and Fortran compilers, so I was looking at both. I know C well but I decided to not look at it as a presumed default. When I used C like this in the past, the intrinsics were very low-level, e.g. wrapping specific Altivec or SSE instructions or whatever. I see see it has OpenMP intrinsics, though, which I’m sure I’ll try later. If I use a library, I’m breaking up the operations and don’t give the optimizing compiler an opportunity to take operations into account together. With Fortran, I can give the operations directly to the compiler, tell it the exact chip I’m working with and it deals with it. It would be fun, when I have some time, to go rewrite it in C and see how it compares.
sampo 10mo ago
> Suppose you pass a parameter, N, and you also would like to pass a tensor, and you would like to specify the tensor's shape (N, N). You can do that, and it might be cleaner and less lines of code that way. But you don't necessarily need to pass the array dimensions as a parameter, as you can call `size` or `shape` to query it inside your function. program main implicit none real :: a(2, 2) = reshape([1., 2., 3., 4.], [2, 2]) call print_array(a) contains subroutine print_array(a) real, intent(in) :: a(:, :) integer :: n, m, i, j n = size(a, 1) ; m = size(a, 2) write(*, '("array dimensions:", 2i3)') [n, m] do i = 1, n do j = 1, m write(*, '(f6.1, 1x)', advance='no') a(i, j) end do print * end do end subroutine end program
- thatjoeoverthr 10mo ago
  True! I just wanted to highlight it. (I had to do this at an interface because the ABI doesn’t pass the shape.)

prennert 10mo ago

I did something similar many years ago. I was amazed that Fortran was not more discussed as an option to write performant code within a Python / numpy codebase. At the time everyone seems to default to using C instead. But Fortran is so much easier! It even has slicing notations for arrays and the code looked so much like Numpy as you say.

3uruiueijjj 10mo ago

I've never found anything to back this up, but my impression was that both the Python / Numpy and Fortran 90 slicing operations were directly inspired by MATLAB (although most of the ideas go back to at least Algol 68). It also helps that Fortran compatibility is a must for pretty much anything that expects to use BLAS.

thatjoeoverthr 10mo ago
I'm learning just now how old MATLAB is. Wow!

foxglacier 10mo ago

Yea, Fortran is nice to use for so-called "scientific" computing. It has high performance as well as some handy intrinsic functions like DOT_PRODUCT and TRANSPOSE but the best features to me are colon array slicing syntax like Python/Numpy and arrays being indexed from 1 which makes converting math equations into code more natural without constantly worrying about off-by-one errors. I wouldn't call multi-dimensional arrays tensors though. That's a bit of a bastardization of the term that seemed to be introduced by ML guys. It wasn't until I started using Fortran that I realized how similar it is to BASIC which must have been a poor-man's Fortran.

PaulHoule 10mo ago
e.g. "tensors" are like "vectors" in they transform in a specific way when the coordinate system changes; what felt so magic about vectors as an undergrad was that they embody "the shape of space" and thus simplify calculations. If you didn't have vectors, Maxwell's equations would spill all over the place. Tensors on the other hand are used in places like continuum mechanics and general relativity where something more than vectors are called for but you're living in the same space(/time) with the same symmetries.
- foxglacier 10mo ago
  Sure, but not all arrays are tensors. In continuum mechanics, you often represent tensors as matrices or vectors to exploit symmetry and then it's really confusing to keep calling them tensors because they don't behave like their tensor equivalent.
- Joker_vD 10mo ago
  > If you didn't have vectors, Maxwell's equations would spill all over the place. What do you mean, "would": they did! :) The original equations had 20 separate equations, although Maxwell himself tried to reformulate them in quaternions. But if you look e.g. at works of Lorentz, or Einstein's famous 1905 paper, you'll see the fully-expanded version of them. The vector form really didn't fully catch until about the middle of the XX century.

pantsforbirds 10mo ago

If your problem fits into arrays/matrices/vectors as the only required datastructures, Fortran is a VERY good language.

drnick1 10mo ago

Any reason why you couldn't do this in C/C++? I use the Eigen matrix library as a replacement for Numpy when I need more performance and the code looks surprisingly very similar.

thatjoeoverthr 10mo ago
Of course I could in C. Intel and Nvidia are both offering both C and Fortran compilers, so I was looking at both. I know C well but I decided to not look at it as a presumed default. When I used C like this in the past, the intrinsics were very low-level, e.g. wrapping specific Altivec or SSE instructions or whatever. I see see it has OpenMP intrinsics, though, which I’m sure I’ll try later. If I use a library, I’m breaking up the operations and don’t give the optimizing compiler an opportunity to take operations into account together. With Fortran, I can give the operations directly to the compiler, tell it the exact chip I’m working with and it deals with it. It would be fun, when I have some time, to go rewrite it in C and see how it compares.

sampo 10mo ago

> Suppose you pass a parameter, N, and you also would like to pass a tensor, and you would like to specify the tensor's shape (N, N). You can do that, and it might be cleaner and less lines of code that way. But you don't necessarily need to pass the array dimensions as a parameter, as you can call `size` or `shape` to query it inside your function. program main implicit none real :: a(2, 2) = reshape([1., 2., 3., 4.], [2, 2]) call print_array(a) contains subroutine print_array(a) real, intent(in) :: a(:, :) integer :: n, m, i, j n = size(a, 1) ; m = size(a, 2) write(*, '("array dimensions:", 2i3)') [n, m] do i = 1, n do j = 1, m write(*, '(f6.1, 1x)', advance='no') a(i, j) end do print * end do end subroutine end program

thatjoeoverthr 10mo ago
True! I just wanted to highlight it. (I had to do this at an interface because the ABI doesn’t pass the shape.)