Hacker News Reader

Comment by sampo | Hacker News Reader

After Fortran got its array syntax 29 years ago, there aren't many languages that you can use to write array-manipulating code like in Fortran. Maybe just Julia and NumPy?

View on HN

All array languages (apl and its derivatives).

Raku (perl6).

R, Matlab/Octave, Mathematica

What are the highlights of Fortran's array manipulation capabilities?

Well, things like... you have 2, 2 dimensional arrays of the same size, A and B.

You want to make a new array, where each cell in A is added to the same cell in B. So you are going to do a loop over 2 variables, right?

in FORTRAN the code is...

C = A + B

in Julia the code is C = A .+ B

Interesting. Does it result in a compile error if you try arrays of different sizes? (If it's just a run time error, I don't see how it's different to implementing this with operator overloading in any language that supports it).

Sure, if you try with non conformable arrays, gfortran for example would tell you something like this:

  Different shape for array assignment at (1) on dimension 1

Note that a scalar is conformable to any array, so that A = 2 for example would work for any shape of A, and would fill A with 2's.

That's what I do in C++ with Eigen or boost::ublas or any other library. Operator overloading is not a feature limited to Fortran.

Fortran gets its performance from prohibition of aliasing, which is something that C++ admits as per the standard. One could enable the same performance boosts with a non standard compiler flag, but this breaks a lot of codebases, including the Linux kernel AFAIK.

Aliasing is a pox upon the world. There's a lot of work gone into stuff like UBSan for C++. I think aliasing should be added to stuff that it checks for. Dunno how much overhead that would be though.

Taking rows and columns is still nice in Eigen. For example taking the second column of a matrix is m.col(1) (Eigen, 0-based indexing) and m(:,2) (Fortran, 1-based indexing unless specified something else).

But more generic array slicing with strides seems to get nasty in Eigen. For example if

    v2 = v(1:2:20)

    Map<RowVectorXf,0,InnerStride<2> > v2(v.data(), v.size()/2);

in Eigen, I don't want to read the docs for how to do the same for 2- or 3-dimensional arrays.

https://eigen.tuxfamily.org/dox/group__TutorialReshapeSlicin...

In Julia, it is C = A + B too. (You don't need the elementwise operator ('.') if the operation is a well defined expression in linear algebra.)

Reply on news.ycombinator.com