Hacker News Reader

Comment by Animats | Hacker News Reader

Fortran is still better at multidimensional arrays. Other than Matlab, few other languages do multidimensional arrays at all. C-derived languages use arrays of arrays, and use templates, macros, or objects for multidimensional arrays. Fortran optimizers have lots of multidimensional array-oriented optimizations which other compilers usually lack.

I tried to talk the Go and Rust people into putting in multidimensional arrays, but the discussion degenerated into bikeshedding, with arguments for fancy but slower representations so you could take arbitrary N-dimensional slices from a multidimensional array.

View on HN

As I member of the Rust team, I personally would like Rust to be a suitable replacement for fortran, and I suspect most Rust core contributors would feel the same. Rust already has competitive performance so it feels like a small step to be able to fit into fortran's niche. I don't know the particular multidimensional array discussion you experienced, but evolving a language takes slow time, and meanwhile people experiment with what is possible out-of-tree, like with the ndarray library, which seems to be well respected.

Right now there is a focus on SIMD in Rust, which is crucial for this domain, so there is progress being made.

> I don't know the particular multidimensional array discussion you experienced, but evolving a language takes slow time, and meanwhile people experiment with what is possible out-of-tree, like with the ndarray library, which seems to be well respected.

A Fortran-like multidimensional library could effectively happen out of tree, but ndarray is closer to Python's numpy.ndarray than to Fortran arrays. And this means that ndarray can not support negative indexing (see here: https://github.com/bluss/rust-ndarray/issues/152). I started a prototype of Fortran like array here: https://github.com/Luthaf/mudi, if anyone want to look at the code or even work on this.

My only problem at the time was that `Index` must always return a reference, and cannot return a new struct. From my understanding, this could be solved with the HKT/ATC work.

> Right now there is a focus on SIMD in Rust, which is crucial for this domain, so there is progress being made.

This is very exciting! And knowing the Rust community, we may even get a cross-platform cross-architecture SIMD library with a nice API!

The SIMD stuff is extremely complex. The focus for now seems to be to _not_ work on a nice API, and tackle the problems with just getting the absolute basics done in an acceptable way. Higher level stuff will come later.

So in my previous life as a physics student I did talk to many physicists about why Fortran was popular (and later when I discovered Rust, discussed Rust with them).

The main blocker is just that fortran has a wonderful ecosystem of scientific computing packages and you can find almost any operation you need out there. This is not true of even Matlab or Mathematica. I have often had to reimplement random algorithms in Mathematica. Fortunately mathematica gives some pretty good high-level tools, so it's not that bad, but it's still annoying.

Multidimensional arrays was brought up as an issue with C/++. AIUI using macros in C++ makes it tolerable, but not great. I did ask about `arr[(1, 2, 3)]` and folks seemed okay with it. Many have used mathematica where indexing is done with double braces (`arr[[1,2,3]]`) so it isn't too bad.

Implicit bounds checking was something that was frowned upon. Could even be a dealbreaker. While usually bounds checking can get elided by llvm, scientific computing does tons of indexing and the overhead might crop up in the profile. I'm not actually sure of how well this maps to the real world, and never got a chance to profile this, but it's a concern.

Libraries like ndarray could turn it off but it would basically be unsafe :|

-------

I did use Rust for some scientific computing. At the time, IIRC numeric types had to be more explicit, which was annoying, but that was pretty much my only issue with it (I didn't need any existing scientific libraries). I didn't know of any good visualization libraries so I would usually output to json/csv and visualize in Mathematica. But I have done the same for Fortran code; Mathematica is just awesome at visualization and analysis.

I recall having access to a scientific computing cluster where the only up to date thing was ifort. It had an ancient python, and icc (IIRC it was ancient too, but I'm not sure). Had to compile a bunch of compilers myself to get my own code to run. Asked about this and apparently everyone just used fortran. Python is common in post-processing but a lot of the tools they were using were ancient or maintained by people with similarly ancient toolchains so stuff still worked on it.

> I tried to talk the Go and Rust people into putting in multidimensional arrays, but the discussion degenerated into bikeshedding, with arguments for fancy but slower representations so you could take arbitrary N-dimensional slices from a multidimensional array.

This is the most frustrating aspect of Rust, to me, as a developer of scientific software. Rust has great potential to challenge Fortran in this domain, but numerics is mostly an afterthought in Rust.

> Rust has great potential to challenge Fortran in this domain, but numerics is mostly an afterthought in Rust.

A common problem with Rust: many people from different backgrounds sees it as a godsend (for some good reasons) and they'd like to replace their former language with it, but the language can't evolve fast enough to satisfy everybody at the same time.

It started as a language for desktop softwares (like a web browser, it's first use-case), and currently the focus is to evolve it to be suitable for back-end development (with asynchronous IO). I hope one day scientific software will become a first class citizen too.

While this is true in some sense, there's also a difference between the language changing to support a use-case, and the library ecosystem. The language doesn't _need_ to change for the async io stuff, though there is one feature that will be extremely helpful, "impl Trait". That feature is also needed for other reasons too, so it's not solely motivated by the async io stuff.

It is totally true that there are many areas that can be improved, and it's gonna take time and effort to knock them out. It's tough! My personal pet area is OS dev...

Can you go into a bit about what specifically is missing?

A lot of the people who work on the compiler/library support for this stuff aren't the primary consumers of it. I personally like working on this kind of code when I have time (I've done a bit of work on the rust bignum library), but it helps to know what's actually going to be useful to people.

Basically, OPs article gives a good overview on what scientists want out of a language used for large scale numeric simulations.

1.) Multi-dim arrays supporting array math, slicing and flexible indexing, sizeof always known (internally solved through lookup tables).

2.) Fast-by-default rather than safe-by-default. A safe language is useless if it gives 2x or more runtime overhead. Make it safe where it doesn't give you any performance penalty (i.e. compile-time checks), otherwise make the default performant or leave the choice to the programmer. Example: Pass-by-reference as default, no GC.

Essentially, I don't think Fortran will be replaced any time soon - maybe strong AI that can program all of this stuff from pure math specification will be here sooner than a Fortran replacement. To everyone not understanding this, read the article again. There's been a huge amount of effort going into this language the last 20 or so years - almost more than justified by the size of the HPC market, which is visible by all the compilers lagging behind the specs. If you want to help, go and improve gfortran (e.g. with better GPU support so people don't need the lagging PGI or locked down Cray compilers anymore).

Agreed. I want to write numeric software, but don't want to use Fortran or C++. Rust looked great, but I couldn't find anything in this domain.

There's a lot of activity for D in this area. For instance, there's a BLAS replacement[1] with good performance. Supporting numerical programming is a priority for Walter and Andrei. There's more work to do, but it is adding enough functionality that it will be a realistic option for those that like the language.

[1] http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...

Why don't you want to use Fortran?

If you want to get started, here's my recommendation: Program Kernels only in Modern Fortran, learn about f2py bindings and do the other parts in NumPy. If you use Fortran purely as a number crunching language - no I/O, no setup, no config code etc. - it's pretty sweet actually - IMO there isn't much need for replacing it with something new that has 50 years of compiler work to catch up (which is near impossible anyways).

Not the OP, but a computational physicist writing a fairly decent sized codebase for a long term project, just went through the language choice and initial prototyping phase. I have extensively used Julia, Python, MATLAB, C99, C++98/11/14, Fortran and D for work and I believe in picking the right tool for the job.

I chose C++ (writing in a fairly modern style, but not eschewing aliasing-restricted pointers, etc. where beneficial) because it handily beats everything except Fortran for performance and is vastly superior to Fortran in terms of zero-cost abstraction. I'm writing code that's largely portable between CPUs and accelerators (mostly GPUs) and has a lot of interchangeable parts. With C++ I can use templates to achieve compile time polymorphism and genericity, while retaining Fortran level performance on CPUs at the cost of higher code complexity.

The Fortran story on GPUs has also been very disappointing initially, especially for those of us not interested in proprietary compilers, though I understand the situation is improving there. Ultimately, in my experience, while individual strengths of C++ can be reproduced in other languages, there's no other language that does it all. I hope that changes some day (for now, D looks like the most promising alternative to me) as writing C++ is often an exercise in frustration, but that's where we're at.

I'm in a similar situation as you, albeit with a large Fortran codebase that needs to stay - therefore I chose to write a transpiler in python that enables a Fortran codebase to be parallelised on CPU and GPU - see my profile link.

Hybrid-Fortran looks like an elegant workaround for an ugly problem, kudos for dealing with the storage order as well.

Ultimately, it seems to me the only way to justify using Fortran for scientific computing nowadays (apart from legacy issues, of course) is if you have no intention of ever using accelerators and absolutely must have optimal CPU performance. NVidia and, to a lesser extent, Intel/AMD could have maintained the language advantages Fortran offers in the accelerator space as well, but chose not to, and now you need solutions like this coupled with proprietary compilers just to approach CUDA C / OpenCL performance. Legacy projects will be around for decades, of course, but for new things I think C++ is/will increasingly be the default choice.

Julia has strong numerical and multi-dimensional array support, but it's JIT compiled. However, once a function is JITted, it has comparable performance to Go & Rust.

Seconded for Julia. I have been using it for math prototyping and it's recently replaced my old favorites like matlab / python / lua.

Yep, Julia is a very nice language with very good performances. But I abandoned it because it is not (yet?) oriented toward library development, and more usable in exploratory style.

I ran into issues a year ago when trying to organize a library with more than 8000 loc, when the "dump everything into a global namespace" started to be more a hassle and less a nice way to get stuff done. I am waiting for it to stabilize a bit before giving a second try.

Could you elaborate on why you felt the need to "dump everything into the global namespace"? For quite a long time now, you can use modules in Julia to keep things in separate namespaces. With "import MyModule" (instead of "using MyModule") things will not enter the global namespace, but will need to be qualified, like "MyModule.myfunction(...)". I feel this is very much like Python. You could even introduce abbreviations for the modules, using "const mm = MyModule" (admittedly, Python's "import...as" is cleaner).

Yes, but `using MyModule` is (was?) the default and preferred way to use a module. Which is equivalent to recommend everyone to do `from module import *` in Python.

I believe the multidimensional arrays are backended by Fortran

Not really. They're just written in a way such that the memory layout is compatible with BLAS libraries — which are frequently written in Fortran. So they can be passed directly to functions that happen to have been written in Fortran.

Additionally, invoking BLAS/LAPACK only happens for dense arrays of specific types (single/double precision real/complex numbers). The great thing in Julia is that almost all of the array and linear algebra functionality is implemented completely generically, so that algorithms can be written using these higher-level concepts, and the low-level functionality will (very efficiently) dispatch to BLAS/LAPACK if possible; if not, the generic implementations still allow your algorithm to work for e.g. your user-defined number type.

haha yeah I guess i should have been more precise with my wording. I'm currently prototyping a custom number system in julia, and obviously, the multidimensional array ops will not be backended in BLAS.

Has anyone tried numerical programming in Ada?

I haven't played with Ada arrays much, but I know it supports various different representations (so it can be binary compatible with Fortran column-first arrays and C row-first arrays), it's got nice array slicing semantics, it's got good support for generics which can operate on arbitrary sized arrays.

I don't know whether it has Fortran-like aliasing guarantees, though. You can't take an address of an array element unless the array has been declared to support it, which helps, but I don't know whether you can e.g. pass the same array by mutable reference into a function twice. (It makes sense that this isn't supported, as it's kinda meaningless, but I can't find an actual reference.)

It might have changed in 2005 or 2012 revisions but here is Ada 95's take on that:

http://www.adaic.org/resources/add_content/docs/craft/html/c...

Also check out ParaSail language designed for easy, safe parallel and concurrent programming. It's an Ada variant from Tucker Taft at AdaCore. Might have something interesting on this topic. Might not.

Rust is like C++ but with better syntax for safer sharing of resources. If you don't like C++ you're not going to like Rust.

The Julia programming language which is gaining a lot of popularity in scientific computation has in fact multidimensional arrays which are very efficient. The whole idea of Julia is to be comparable to C and Fortran in speed while having the ease of development of Python.

I'd dabbled in Go and Rust as well but don't see them as serious competition to Julia in this space. As much as I love Go, it is really poorly suited for numerical computing. No operator overloading of any sort among other things. You really want that for vector based math.

Also Rust, seems like more of a language for software developer professionals. I think Rust looks neat, but it is too complicated for scientists I think. They don't want to focus that much on their language.

In this respect I'd actually think Swift has more appeal as it is quite fast, but easier to use. By easy, I mean there are fewer concepts to learn.

But Swift will have the same problems as Go and Rust, in that it is not being designed with the needs of people doing scientific computing in mind.

That is what is nice with Julia. Out of the box it does so many of the things you want. Right out of the box without explicitly importing anything you can start doing advance matrix operations, load large CSV files, groks imaginary numbers. It actually has a number system which distinguishes between rational and irrational numbers e.g.

> I tried to talk the Go and Rust people into putting in multidimensional arrays

What facility is Rust lacking that prevents you from coding this as an extension/package rather than needing core coding support?

Native long integer (128-bit) support could be coded as a compiler extension for a long time before it got pulled in.

What facility is Rust lacking that prevents you from coding this as an extension/package rather than needing core coding support?

That sort of thinking is what drives number-crunching people away from a language. You end up with ten different representations for matrices, and can't pass matrices from one package to another.

If the numerics people can't agree on something, throwing the decision over to the core Rust folks is unlikely to result in something better.

So, the question still stands: What in Rust prevents an ambitious number cruncher with good design taste from adding the appropriate pieces to the language?

Fortran has 1st class many-dimensional arrays; can you really do that efficiently (as in "works great with a loop nest optimizer/LNO" and "has pretty syntax" and "there's one way to do it") as a Rust extension? C++ says "NOOOOOOO!!!" I don't know Rust well enough to answer the question.

You can do `arr[(1,2,3)]` in Rust (not pretty, but ok). Assuming you've marked the indexers as inline, LLVM should be able to optimize it pretty well. If it doesn't, optimizations can be added.

Rust can't enforce "one way to do it" since this would be done as a library. Currently there's only one library offering this (ndarray), so in a sense there's only one way to do it :)

What does arr[(1,2,3)] have to do with what I asked for? The code that I write for my personal scientific stuff (numerical hydrodynamics) needs loop jamming/fusion and splitting, SIMD without having to use any special notation, and other stuff which is typically provided by compiler LNO. I happen to have interchanged loops in my code, but if I goofed that up, it'd be nice for the compiler to do it for me.

No array slices needed. Plain Fortran 77 notation provides what I need.

> What does arr[(1,2,3)] have to do with what I asked for

"has pretty syntax", was just showing what the syntax would look like in Rust.

LLVM autovectorizes, and does many loop optimizations. I suspect LNO exists, but if it doesn't there's nothing preventing it from being added afaict.

Was mostly explaining what the status of this in Rust is. It's not awesome, but not horrible either.

Thanks for trying to be helpful. The place where C and C++ goes wrong is multi-dimensional arrays. Some googling tells me that narray gives Rust syntax like foo[[1,2,3]] for a 3 dimensional array. I have no guess if this optimizes well or not.

It will compile down to `foo[1][2][3]` and get optimized the same way. Rust additionally does hint a lot more about aliasability to LLVM (everything is restrict by default, basically), so there probably is more scope for optimization here.

I'm not an expert in that domain, but to my mind, the two big things are

* lack of type-level integers

* lack of explicit SIMD on stable Rust

This is already true in the Fortran world. When I was a student I was working with some condensed matter physics software and there was all kinds of madness with flipping conventions in there.

And Rust doesn't have ten different ways to do n-dimensional arrays. There's one major library out there (ndarray). I just learned about a second (https://github.com/Luthaf/mudi) which seems to use the same conventions.

There are multiple ways to do matrices, but these are in graphics-oriented libraries.

Well, it can work if one package gets dominant - e.g. the Python ecosystem works quite well with having all this outside of core language in the numpy package.

But that happening is just as likely as the fragmentation that exists in C++. Getting OpenCV to play nice with Eigen can be a hassle.

Nothing prevents a person from doing that, but that isn't what folks doing scientific computing are interested in, generally speaking. Some are, sure, but in aggregate as a matter of what the "day to day" is, it's just not of much interest.

I mean, maybe I was an outlier, but when I was in academia I wanted to have my models compute as accurately and quickly as possible, and I generally didn't have much time for or interest in the computer coding beyond the necessary implementation details required to conduct my research.

Telling me "implement an extension", then, would have made me roll my eyes and go back to fortran.

I don't think the idea is necessarily that you personally implement an extension, but that "someone" needs to implement it and the bar for "someone" to do so is lower if it can be done via extensions or libraries rather than having to be done in the core. That someone should share their work, of course, and that someone may find they have a feature or two they have to get into core, but even that's easier than sticking a whole whackload of changes into core.

You don't have to implement an extension. It already exists.

https://bluss.github.io/rust-ndarray/master/ndarray/index.ht...

It would be nice to have access to an overloadable index operator that works on multiple dimensions at once (like C#'s [1, 2, 5] syntax). It's sometimes not easy to write something to be efficient if it has to go one-by-one like [1][2][5].

The ndarray library indexes via tuples, i.e. [(1,2,5)], which should be just as efficient as a built in syntax (if perhaps a bit uglier).

Doesn't Rust have macros that let you modify the syntax? Maybe you could do this.

One of the nice things about D is that it offers multidimensional indexing. For instance, you can write x[1, 2..6], and you can do it for user defined structs. In my case it makes it really nice for interfacing with R.

Fortran is also better at complex numbers. In Fortran these are a native type with every bit as much support as for real floating point types. Not so in the top competing languages (C/C++).

C has native support for complex arithmetic since C99. What kind of support do you think is missing?

http://en.cppreference.com/w/c/numeric/complex

Last time I looked at C, the libraries had yet to be ported to C99. You had to choose between z = u + 2 * i * w, or zs = fft(ws); you couldn't have both. Therefore Fortran kept winning.

It is kind of a "too little, too late" problem. C added complex numbers a decade after C++ spun off, so the two languages have different approaches. And C99 complex number support was slow to be added to compilers. On the UNIX side, it was never added to Solaris 9 (EOL 2014); Microsoft finally added Visual Studio support in 2013. And even if your compiler supported it, code generation quality has been historically bad. Some things were never added to the standard, like printf support, so complex double is still a lesser native type than double.

In C, I've seen people do "multidimensional arrays" by A[n * i + j], with A a double[m * n], rather than by arrays of arrays. Especially if the array is being passed as a parameter, and the array's dimensions are given by other arguments, and the array is to be used with BLAS.

Is the compiler smart enough to do strength reduction and get rid of the multiply when indexing? FORTRAN compilers have been doing that for multidimensional arrays since 1954.

Reply on news.ycombinator.com