Rewriting Fortran Software in Rust | Hacker News Reader

Rewriting Fortran Software in Rust

mckeogh.tech

151 pointsfanf276 comments

I used to maintain a large FORTRAN codebase - it seemed like every time new grad students showed up, everyone got a wise idea to re-build everything in a modern language.

This question is revisited a lot, but the general consensus has been that FORTRAN is fast, simple, easy to understand and code, and the compilers used are super optimized (Intel's FORTRAN compiler was really a gem - it managed to do a ton of automatic parallelization on modern hardware)

I'll see if I can find it, but I remember attending a conference talk at AGU in ~'08-10 called "Do-over or make-do" [1] that analyzed how many people-hours it would take to update all the climate modeling software (earth and space sciences are really huge on legacy FORTRAN code), and the conclusion of that talk was: use modern code/tools for glue, keep your climate models in FORTRAN.

There's a ton relating to precision - we had a huge ordeal converting our 32-bit space weather model to 64-bit because the change in precision changed our results, so we couldn't publish (in good conscious) without making sure the results weren't within a good margin of error. Anyhow.

This is always a fun thing to revisit now that I'm really far away from having to maintain any FORTRAN code any more :)

edit found it! [1]: https://www.easterbrook.ca/steve/2010/12/agu-session-on-soft... [2]: slides http://www.cs.toronto.edu/~sme/presentations/Easterbrook-AGU... Also randomly found this: http://www.moreisdifferent.com/2015/07/16/why-physicsts-stil...

I came here to ask whether your typical Fortran program would carry the coding security risks that Rust is designed to mitigate; I couldn't think of any. Your comment seems to confirm my doubts.

It probably wouldn't have affected my work at all. We were running our code as a monolith - it was the only thing that ran on any given machine, and its output was stored on a large filesystem on NFS.

That's not to say there isn't someone out there with this concern :)

https://www.fortran.io/

For example.

More seriously, I could imagine a security flaw in numpy wreaking all sorts of havok. But it's also somewhat difficult for me to imagine that there are any security vulnerabilities left in BLAS. But then I remember Heartbleed.

I think „coding security risks“ referred to language properties that increase the likelihood of introducing bugs in a general sense, not in terms of cyber security/hacking.

Memory errors also affect correctness.

Forgetting to initialize that last entry or initializing one too many entries by accident can cause random, very hard to track down, issues.

Security? Perhaps if your code is pulling chunks of data from a cloud source?

You are right, I forgot to mention correctness.

> we had a huge ordeal converting our 32-bit space weather model to 64-bit because the change in precision changed our results, so we couldn't publish (in good conscious) without making sure the results weren't within a good margin of error

This is a really interesting dimension to the transition that I've never considered before. An actual change in correctness, not just compatibility. I would have thought for floats it would just add more decimal points

FWIW, I didn't bother getting too deep into an analysis (I was an undergrad and not really savvy enough to deep dive, but maybe I could now if I went back :) )

Our model was a space weather model of the solar wind speed and density - it was very iterative, and I think that iterative nature is what contributed to quirks in the transition - think of the issues when summing say, 0.1 ten times in Python (or other languages) - ultimately, we didn't need to do any real work to show that our converged model was within a margin of error and still presented the same results in terms of predicting solar flares - but it was still a lingering curiosity that the numbers weren't exactly the same (as you say!) - most of the work was done in finding the right compiler parameters to not break existing code, while taking advantage of the larger memory space.

Ah that makes sense

One of our numerical analysis(NA) professors had worked at NASA, at IBM and had enormous real-world experience in NA. Several weeks into the class he challenged us to write a correct square root routine in FORTRAN. Over the next two weeks submissions were presented in class, where our professsor showed how each had multiple critically fatal errors. He let it stew awhile: two weeks before semester's end he presented and explicated to us a correct implementation of the square root, which turned out to be far more lengthy and involved than we had imagined. He said there was always the possibility of more bugs being present. That was really an eye-opening class.

Post the code here, and educate us.

musl sqrt is here, about hundred lines of code with nice explanation on top. Considering that it has 1993 copyright stamp, I'd imagine it is fairly canonical implementation

https://git.musl-libc.org/cgit/musl/tree/src/math/sqrt.c

I saw that done years ago in assembler for UNIVAC mainframes. Tear apart the exponent and mantissa. Shift mantissa if exponent is odd. Shift exponent right one bit to get square root of exponent. Look up starting guess from table using a few high bits of mantissa. Run two iterations of Newton's method. Reassemble mantissa and exponent. Normalize. Done.

Is it possible to keep ulp precision in all cases after doing all that separation, normalization, and recombination?

Although Fortran compilers may produce fast code (they have to - there's no chance to play with SIMD intrinsics or templates unless you want to use C++), I don't find them very good. I've personally hit several bugs in gfortran. The code I maintain now has some extremely ugly hacks to work around unfixed compiler bugs. I recently tried compiling it with Intel Fortran and got a compiler crash on one input file. Maybe they work well with Fortran 77 level code, but the Fortran 2003 level support seemed rather buggy. I suppose this is the problem with a minority language.

I won't recommend Fortran unless the code has to be written by scientists who know nothing else. It has horrible string support, no data structure or algorithm library (like STL), no templates, few non-numeric libraries (e.g. http, json...) and (as above) buggy compilers. Furthermore, C++ gives you much more control over numerical speed when you need it.

"It has horrible string support", this is untrue. Learn Fortran 2003, 2008, 2018.

What? Where is this good string support in Fortran 2003, 2008 and 2018? I must have missed it in the manual. Now you can allocate a string of a non fixed size - simply amazing technology! Actually, one of my compiler crashes was trying to use this revolutionary feature.

> Now you can allocate a string of a non fixed size - simply amazing technology!

That’s more than what C gives you out of the box ;)

Many languages forbid mutable strings. They aren't really a very good idea, imo. Even less of a good idea are extensible strings. More complex data structures are required for efficient extensible strings, and this complexity shouldn't be hidden in a simple interface. Better to have a Buffer type or something that is distinct from the String type so programmers are made aware of the performance/allocation characteristics.

Though, this is kind of ridiculous to even talk about, since the CPU and compiler itself are responsible for the greatest sins of hiding something unimaginably complex behind a simple interface.

Ah, the joys of the 6502 and the Commodore 64!

Delphi had mutable xor aliased strings. That is really great. Each function can change the string as fast as a mutable string, but once it shared with other functions, they can be certain that it will not be changed randomly

Or just bite the bullet and learn how format statements really work. ;)

We could never use gfortran, and other open source compilers were not really great for us either. We used Intel's fortran compiler, and it stood out as a pure gem. At the time it was free for educational/research purposes, which is why we used it, but if we had to pay (I think the license was like $1k/yr) we would have done it in a heartbeat.

Our code was entirely FORTRAN 66 and 77, so that may explain why it worked so well for us.

Have you tried Flang? It’s LLVM's fortan compiler and was written from the ground up to support FORTRAN-2018

Nope, I left the field in 2012 so haven’t kept up to date at all, sounds cool though!

I wonder if when global warming starts eating coastal cities will the government have to seek out centenarian fortran programmers to analyze the climate models?

(like the new jersey unemployment system in cobol)

I'll do the same thing they did and ask for a shit ton more money :)

Though the beauty of fortran to me is really how easy it is to learn and make working code. That said - it's not elegant at all, probably a feature.

Ditto this, I've run into some strange compiler bugs I've had to work around. Intel Fortran 2019 actually fails to compile a code at work that's worked for ~40 years, Intel Fortran 2018 compiles it but produces numerical errors when the code runs. Intel Fortran 2016 actually works perfectly, as do several earlier versions.

With gfortran I run into issues primarily when using the "newer" features of Fortran 2003/2018.

Julia is taking some inroads in climate modelling: "when the coder-climatologists of the Climate Modeling Alliance (CliMA) — a coalition of US-based scientists, engineers and mathematicians — set out to build a model from the ground up, they opted for a language that could handle their needs. They opted for Julia."

https://www.nature.com/articles/d41586-019-02310-3

https://clima.caltech.edu/

That's awesome, and not at all a surprise! It makes a lot of sense when building from the ground up to choose a great tool, and Julia has had some serious adoption in academia.

I think long term it has enormous potential to take over as _the_ language for modelling.

The reason is that for every hyper optimized climate model there are a thousand small or medium models. For these Julia offers sufficient performance and excellent developer experience.

It probably helps that Julia and Fortran are both column major.

When interfacing with C languages, that can be a big penalty (depending on the complexity of the algorithm).

It will probably take two generations of developers to make this happen, but well worth it. When I was studying atmospheric science I tried to understand OpenWRF (weather research forecasting) and the barrier was too high... And I had programming experience, imagine those with none!

An open source project is only as useful as it is accessible to it's average developer base.

"in good conscience"

Nice catch, thanks!

> There's a ton relating to precision - we had a huge ordeal converting our 32-bit space weather model to 64-bit because the change in precision changed our results, so we couldn't publish (in good conscious) without making sure the results weren't within a good margin of error.

That sounds like a question of numerical stability of the algorithm, rather than anything to do with the programming language, no?

Yep. Fair point. The folks I worked with were very pragmatic so didn’t get into the internals and just treated the code as the instructions. Many of them came from an earlier era where the code was a lot more specific and the concept of an R8 is not quite as specific on x86_84 as it was on a VAX because the modern CPU does a lot of tricks under your nose.

Had you done some kind of external validation of the original results? If not, surely the new results are just as valid to publish as the original ones - any bugs are as much or more likely to be in the 32-bit than the 64-bit version.

Yep. Our model was run by other organizations and incorporated into amalgamated models that had further downstream error bars that they needed.

Not that our model contributed much but we were one tiny part of the early warning system for solar flares for satellites and astronauts.

We admittedly weren’t super formal about the transition but wanted to make sure the differences wouldn’t be meaningful or propagate into other models we contributed to.

Also fwiw, we had no control or say over how others ran our model so they could have easily done their own 64 bit build of our code but I think it’s unlikely.

It’s nice to see an article that makes assumptions and then realizing that they didn’t really work out and walking back on them. It’s often good to try rewriting things, but often the thing you think is old and crusty (FORTRAN) is actually quite good, perhaps in ways you weren’t aware of. So it’s important to be able to say, “yes I tried some nice new thing, it didn’t really work out how I had expected, also some of my assumptions I started off with were not really correct either and I might not have made a fair comparison”.

I find it remarkable that one person with just over a year of Rust experience can successfully rewrite ~6000 lines of well-tested, highly tuned Fortran in Rust and end up with similar performance on benchmarks without having in-depth expertise in the domain (simulation of fluid dynamics). The benchmarks (https://mckeogh.tech/post/shallow-water/84896735-54a97580-b0...) show similar performance across all input sizes for both the old Fortran and the new Rust code, with Fortran outperforming in some cases and Rust outperforming in others.

FTA: “The result of this is that these simple operations are being performed on some very large amounts of data which I believe is resulting in a memory bandwidth bottleneck. This hypothesis explains why SIMD did not improve performance and performance is not significantly improved by parallelisation”

If that’s true, comparing system load or power usage could be a way to differentiate between the two implementations.

Rust is not a very good example here considering the amount of time the op might have spent waiting for the compiler and steep learning curve of the language...

Actually closer to 4 years of experience:

> As for my Rust experience, I’ve been using it for personal projects since ~2016 and I worked as a Rust software engineer at a startup in Berlin for a year after leaving high school.

Publish date on this blog post: July 14th, 2020

This is an interesting article. I have often wondered if there would be benefits to porting scientific simulation codes to more modern languages. In the end I've always reluctantly stayed in C or Fortran for scientific simulations. The article touched on the some of the issues (and the ever present issue of "If it isn't broke, don't fix it"). But I do think that the article ignored what I think is the biggest issue.

To me the biggest issue is that some simulations just need to run on 1000s or hundreds of 1000s of cores, and the primary technology that lets this happen is MPI (Message Passing Interface). I've been writing new Fortran code only because I want the simulation to be able to scale to this kind of size. The issue is that MPI only has official bindings for C and Fortran, so the only real choice I have is to keep programing in C or Fortran if I want to use MPI.

Any thoughts on this? Personally I would love to move towards more modern languages but without the proper libraries it is still a tough decision.

Rust can bind to C libraries, as can many other languages.

I have looked into this a bit, but the issue appears to be that MPI is not a library per se but a standard interface to more underlying function calls. Different implementations of MPI are in effect different libraries using the same API, and because of that it is difficult to write general bindings rather than just bindings to a particular version of an MPI library that is implemented on a particular supercomputer. In short, it's too site specific and fragile, unfortunately. If you change the system, the program might no longer compile. I appreciate the comment and feedback, though.

This Stack Overflow page comments more on this issue:

https://stackoverflow.com/questions/22949462/rust-on-grid-co...

MPI officially supports only C and Fortran. Anything else for any other language can be only considered as adventures of some open-source developers, that may not last forever. Take a look at the fate of the C++ bindings for MPI.

Including Fortran.

I don't know that I have any thoughts, but I do have

Consider Julia. It has parallel processing built in, and uses MPI under the hood.

A bit more specifically: IIRC it defaults to trying to use MPI when you launch distributed jobs. It can also fall back to using SSH, for "extremely limited cases" and quite frankly it's not inappropriate to use SSH if you are certain that communication is not going to be a bottleneck.

For those sorts of scientific simulations, you probably want to use the MPI.jl package rather than the built-in parallelism – at least that’s the conclusion I came to when I started the fluids code I’m currently working on. MPI.jl is a thin wrapper around the C interface so you end up writing your programs pretty much the same way as in C/Fortran. This type of MPI support is still a bit rough around the edges but it works, and our code runs as fast as the Fortran code it is based on (we usually scale to ~1000 cores). You get the advantage of working with a very nice, modern language, but it is pretty easy to accidentally do things that are very expensive and killing your performance, such as allocating memory inside a hot loop.

We use Akka in java for related distributed computation use cases.

This is a really cool post, but I wish the author had profiled his code. He makes several statements like "I think the performance issue is XYZ" but never bothers to check if his guess is correct. As Walter Bright liked to say, Use the profiler, Luke!

Nice article. I'm curious which Fortran compiler did you use? One key difference that I see in the context of this article is that Fortran is natively parallel (both shared and distributed memory) via coarrays, teams, and events, whereas Rust (and Rayon), as far as I can tell, is shared-memory only. For applications that need to run on 100s or 1000s of cores, Rust needs to bind with C MPI. Is there any other approach to distributed-memory parallelism in Rust?

Good points on the lack of certain features and ergonomic libraries. We're working hard on making this much better for Fortran [1] by developing a community-driven stdlib, package manager, and similar tools.

In fact, we'll have our monthly community call tomorrow [2] where we'll discuss the best way forward for a nice strings API in Fortran stdlib. I welcome everybody interested to join the call.

[1]: https://fortran-lang.org [2]: https://fortran-lang.discourse.group/t/fortran-monthly-call-...

This is interesting, I did write BLAS in Rust https://github.com/schultzer/libblas it's only single threaded at the moment.

I feel people forget to ask a more important question before they start a full rewrite and that is: Where is the hardware trends going and, is dedicated GPU, CUDA etc. even the future.

I feel it would be wasteful to spend a lot of time writing GPU specific code, Its not really maintainable on these scientific project (where only a select few people contribute, usually people that do the research).

The author says he's happy with the performance improvements he made, but--even at the size where the Rust version performs best--he gets much less than a 2x speedup. After a couple months of effort, that's not a whole lot to show. This reads like a case study in why it's probably a bad idea to rewrite even when you think you have compelling reason for it.

Also the author has not used state of the art compilers like XL, Intel and PGI.

Indeed, and this depends between applications. For the WRF (Weather Research and Forecasting) model, I get 3-4x speed up with Intel Fortran compiler over gfortran. I don't see where the Author mentions the Fortran compiler used, though.

> I get 3-4x speed up with Intel Fortran compiler over gfortran

Wow, even at similar optimization levels? -O3 for each?

Not quite, but at optimization levels set by default WRF configuration for each compiler (definitely not a fair comparison):

gfortran: -O2 -ftree-vectorize -funroll-loops ifort: -O3

I don't have the timing results anymore. This was in 2018 on Xeon Platinum 8168.

I recently tried replacing -O2 with -Ofast -ffast-math to gfortran settings, which gives about 18% speed up. So still far from Intel. I recently proposed it here [1].

[1] https://github.com/wrf-model/WRF/issues/1254

Have you posted about this on the gfortran mailing lists? They are generally interested in examples that show where the compiler could be improved.

I doubt they are short on inspiration for further improvements. A big part of the speed advantage of ifort is the aggressive loop unrolling, pipelining, splitting and multiversioning.

Coming up with these transformations is not the difficult part. That would be actually implementing these transformations correctly and keep the whole compiler optimization framework maintainable.

And of course there is certainly a cost for the improved runtime. Primarily in compile time and code size. For the Intel compiler it generally makes sense to sacrifice these in favor of runtime performance, because it is used on these kind of scientific codes a lot, but for gfortran the balance might different.

Already mentioned on the Reddit thread, the author apparently isn't aware of modern Fortran capabilities, IDE tooling support, NVidia and Intel offerings for GPGPU programming or HPC.

PGI's Fortran compiler, bought by NVidia a few years ago, has excellent GPU extensions. You set an attribute on an array, move stuff on there, do your computation, move stuff back. It was the easiest GPU coding experience I've had.

Nice post!

I have written several versions of a scientific simulation program in python, first using MPI (single machine) and then native multiprocessing. Using numpy and Pandas for all the heavy data operations. I really valued the ergonomics of python and native multiprocessing and found the MPI code to be not very human friendly. However, I ran into some very strange bugs relating to MKL and multiprocessing that just hung the program on large input sizes with MKL_NUM_THREADS above 1. I also had to bounce the input data to the worker processes to files due to mp not being able to pickle data when it got too large.

As soon as Apache arrow support in rust gets more feature complete I'll probably try to redo it there.

TL;DR: wanted to leverage GPUs for computation (but couldn't), the ergonomics of arrays in Rust for scientific computing without the use of ndarray are not good, had memory bandwidth issues to keep cache lines fed, performance of the two implementations is comparable

One of the things I find interesting is how the Fortran version is faster both with small and huge inputs, while the Rust version is faster for input sizes in the middle.

In my experience... Reinventing (or rewriting) the wheel is a fun exercise. It’s even more fun to write it in rust (biased). Occasionally, the new wheel is better. Mostly, I learn how some of my assumptions were wrong and/or how I didn’t even account for “this, that, and the other.”

That said, I enjoyed reading through and finding “Edit: I was wrong about X.” Learning experiences all around!

Also, might have been mentioned but is this a typo?

> I’m guessing that FORTRAN is faster on small inputs due to the threading overheads, but faster on large inputs

Anyways, thanks! Enjoyed the article!

> Having written parallel software in both C and Rust, the memory safety guarantees and easy parallelisation with Rayon offered by Rust contrasts quite sharply with my poor experience using OpenMP

While I appreciate the strive for safety, and the effort being put into the rewrite, I believe in that particular case it would be more beneficial to rewrite it in ATS with safe formally-proved C rewrites/known to be safe inlined C calls.

It seemed like this one area was a definite plus for the rewrite in Rust. I recently added some OpenMP goodness to some C++ code. Making the loops run parallel was almost trivial. The challenge was going over a couple thousand lines of code checking for safety. There were data structures being created and shared in various was but the code was written in a very functional style so it worked in the end. With Rust, all that analysis is done by the compiler. So for that code migrating to Rust is probably a good idea.

ATS is not readily usable for large programs. Among the "popular" programming languages Rust is perhaps most likely to get some support for interactive formal proof down the road, given that the safe subset of Rust basically enforces the rules of functional programming (and the unsafe subset is self-contained in a way that's not necessarily true elsewhere). Work along these lines is already happening as part of RustBelt.

“ I’m guessing that FORTRAN is faster on small inputs due to the threading overheads, but faster on large inputs because of contention between threads for cache lines and memory bandwidth"

I think the author meant FORTRAN was slower on small inputs?

Look at the graph (or the chart on the left, I had to to understand it): the green and blue lines cross twice, first the FORTRAN version is faster, then the Rust one, then the FORTRAN version again.

New engineer thinking, want to do everything themself

why? Why would you do so? Did you gain any significant improvement in your final results compared to the original software? Why would rewrite a software that already works fine, other than not having the will to learn a bit of a new language? Why???

All your questions are addressed in the article. It does turn out that some of the author's beliefs were incorrect (e.g., the state of Fortran GPU compute libraries), though, so it's possible that they would have come to a different conclusion had they known that.

Reply on news.ycombinator.com