Glad this was posted. First thing I did on loading the comment section was ^F "alias". Fortran is still used for numerics because it makes writing fast numerical code easy. C simply wasn't designed for it, just as Fortran wasn't designed for low level system/OS programming.

In C99, you can use the "restrict" keyword to tell the compiler an array is not aliased, and even without restrict it's often possible to get the same optimisations in C or C++. When I was still doing HPC in academia, all our C++ code was maxing out the memory bandwidth. I think in practice the difference rarely matters these days.

Insomuch as I've seen HPC code written by science academics, both when I was one and later when I became an HPC guy, as much as I cringe at some of the Fortran, the C/C++ can be far worse... and I wish that I could retroactively restrict some of my colleagues and customers to sticking with Fortran.

Wow, look at the downvotes! I've been saying this to people for 20+ years, in conference talks and face-to-face conversations, and I've never had a negative reaction. Happy to learn if folks think I'm wrong about my personal experiences. If you've had a different experience, happy to learn about that too (although I'd be surprised if you'd downvote me for that.)

(edit: thanks for the upvotes! and discussion!)

Personally I'd rather see academic c than academic fortran; academic c++ is in the 'oh please no' category.

Yes, I've seen too many unfortunate template libraries, which never work in combination with a different template library.

I've used some quite reasonable academic C++ for geometry processing and meshing.

Note that I'm not saying that science academic C++ is never great; I'm saying that it's often unfortunate. Is the code you're talking about code that's progressed through several maintainers without becoming sad? That's often where things go wrong, not the initial code, but keeping it going, cleanly, through multiple collaborators / post-docs / grad students.

I agree that meshing is definitely a difficult situation; I've never seen a solution (for unusual geometry) that I liked much. Things that can pretend to be rectangular? No problem. Very adaptive meshes? Argh.

Hard to say, but I'm looking at 39 contributors, 539 stars, 3,005 commits [1]. Nicely documented, pretty extensive list of users, CI on Linux and Windows. But with that said, yes, I've also dealt with absolutely horrendous academic C++.

[1] https://github.com/libigl/libigl/

Are there any users that are science academics? I wasn't talking about computer scientists / computer graphics / visualization researchers at all in this thread; they do great with C++. I clicked on a bunch of links of the users of this library and didn't see any physics, astronomy, chemistry, fluid dynamics, etc.

Not sure! We are using it in (non-academic) CFD work.

So it sounds like you're somewhat like the folks I'm talking about, even if the other users of this code are not... but have you had a lot of turnover of your coders, like grad students and postdocs that turn over every 2-4 years? Because that was the full nature of my complaint: science not-computer-scientists over several maintainers doing better in Fortran.

You can also quite simply use local variables, and unroll loops manually:

   for (i ...) {
     double a0 = a[i];
     double a1 = a[i+1];
     double b0 = b[i];
     double b1 = b[1+1];
     c[i] = a0 + b0;
     c[1+1] = a1 + b1;
   }

The hand-unrolling is predicated on the knowledge that the arrays do not overlap. It is obvious to the compiler that the locals do not alias with the arrays, or with each other.

The stores to c[] could alias with a[] or b[], but that hardly matters; we assume that all the array loads and stores take place. We have done ourselves the unrolling that was hindered by the suspicion that the arrays overlap, and could do more of it.

When coding in C for tightness, it helps to pretend that locals (at least basic type scalar ones) are registers and references to arrays and structs are memory loads and stores. Then think like an assembly language programmer, somewhat. It's better to use a few more locals than a few more loads and stores. Try to load once, process with locals, store once.

yes, but restrict is limited in that it can only be used to tell two pointed memories could never alias. Finer grain aliasing can be done using unions and structs, but it's not as flexible as Fortran and the possibility of non-aliasing is easily broken just by allowing address to be taken.

As an aside, how to you accurately measure memory bandwidth? I've been looking for a good solution but couldn't find one.

Some processors expose performance counters but you don't need to accurately measure memory bandwidth to know you are maxing out. Suppose you know that the memory bandwidth is 1Gb/s and you program reads 100Gb of data to produce the result. If the execution takes ~100s then you know you hit the maximum memory bandwidth.

Were you using C99 during your HPC programming in academia?

I used C99 for the computing kernel of one standalone project (simulation of gravitational lensing). Apart from this, it was mostly C++.

What made your work so dependent upon memory bandwidth?

As an introduction, you might want to look into John McCalpin's history as a shallow-water oceanographer and his subsequent fixation on the STREAM benchmark. TL;DR: some algorithms aren't friendly to caching.

I have no idea if this applies to smarnach, but it definitely applies to other people.

John McCalpin did a nice presentation at SC16 giving an overview of memory BW/latency issues in HPC, slides: http://sites.utexas.edu/jdm4372/2016/11/22/sc16-invited-talk...

Your point is true but I thought you Fortran users might find the OS below interesting since they wrote it in Fortran. Only one I know of.

https://en.m.wikipedia.org/wiki/PRIMOS

Nice bit of trivia! I was an admin of a Primos machine briefly while also doing scientific programming.

Primos wasn't much good compared to the commercial Unixes that were appearing. It did have one nice-ish feature late in it's history; the CPL scripting language.

Wow, big trip down memory lane here because my first job was on a Prime mini-computer! Wrote a complete system using CPL scripts (with the occasional bit of Fortran thrown in) on Primos.

Shame to hear that Prime / Primos is now defunct but that job was back in the 80's!

What does non-aliasing mean in this context?

It means that the compiler knows that pointers/arrays don't overlap in memory. Data can therefore be safely stored in registers, instead of continually re-read from memory, to account for possible alterations between computations [1].

[1]: http://stackoverflow.com/questions/146159/is-fortran-faster-...

I wonder if Rust is fundamentally able to take advantage of its borrow checker and outperform C in this regard.

As a counterpoint, it is worth noting that essentially no high-performance BLAS or LAPACK implementations (e.g., GotoBLAS/OpenBLAS, ATLAS, BLIS, etc.) are implemented in Fortran and instead tend to use a combination of C (in order to use vector intrinsics) and assembly.

To be fair, there are GCC built-ins that allow you to tell the compiler that arrays are nonaliasing, but you're right in that it's more difficult in general to get a compiler to spit out fast numeric code with C than it is with Fortran.

I feel like it's a "right tool for the job" situation. Fortran was literally built for the purpose of numerical computation. Multi-dimensional arrays, non-aliasing arrays and pointers, multi-dimensional array intrinsics, intrinsic array binary operators... all come with Fortran right out of the box. Of course, C/C++ can do all of these things, but all of the proposed methods, at least that I have seen, feel like they are trying to hammer a screw.

I have one question - if I compile and link a program written in C and Fortran, can the linker inline Fortran functions in C? I know it can do it if I link object files produced from C source, but can it do it if I link object files from separate languages?

Most Fortran compilers are actually Fortran/C/C++ compilers, and cross-link and cross-inline with no problems.

Having never done it with Fortran, I don't know. It should be possible, because e.g. on Linux, everything adheres to the ABI. However, the ABI doesn't specify how to interpret pointers AFAIK. Array ordering [0] and indexing is different in Fortran compared to C, so I'm sure there's pitfalls you'd have to account for.

[0] https://en.wikipedia.org/wiki/Row-_and_column-major_order

Fortran arrays can be indexed however you like. Code that is designed to be called from C can be written with 0 indexed arrays.

Sure, but you still have to remember to do so :)

It's not the link editor's function to inline; an optimizing compiler does that. A linker can, eventually and depending on the platform, perform some binary optimization, as is for example the case with Sun Studio compilers on SPARC (-xlinkopt=2), or with GCC (-flto), but that's not inlining.

To nitpick, compilers nor linkers typically are not able to do "binary optimization", if by that you mean analyzing the asm and optimizing that. What happens when you compile with -flto or equivalent is that the compiler hides some internal representation in a special section in the object files, then at the link stage (if you also use -flto when linking) is that the linker invokes the compiler again which deserializes those trees from the object files, and does the inter-procedural optimizations on those.