Weather and climate modeling codes from Fortran to C++

42

A more user-friendly interface for LAPACK/BLAS/etc libraries, I only use it for solving linear systems (for now at least) and when starting to learn with inherited codes I couldn’t find much sense to that “weird” dgsev call until reading further. Now my students have the wrapper function solve_linear_system

So the call

call dgesv(n, nrhs, A, lda, ipiv, b, ldb, info)

becomes

X = solve_linear_system(A, b)

This is a really simple/basic example, but I think that correctly wrapping this routines (being more similar to numpy in some way) is a good approach to get a better Fortran ecosystem

3 Likes
43

I see that @FedericoPerini wrapper implementation is already doing this, makes me really happy!

2 Likes
44

@fedebenelli that is exactly right. The high level design:

  • if you just want to solve a system (and stop the program if it can’t be solved), x = solve(A, b) will work
  • if you want to handle errors yourself, you can do x = solve(A, b, err) and then you can handle err yourself and have full flexibility.
1 Like
45

I think the 2nd case should be written

call solve_system(a, b, x, err)

although it’s not a big deal, since one can easily write a subroutine wrapper to do this.

Richard Maine (editor of the Fortran 2003 standard) explained in comp.lang.fortran 16 years ago why he avoids having functions return error flags.

I’m even less a fan of having functions return values in arguments.
I’d probably recast the function as a subroutine if an error return was
a possibility. After all, to me most of the point of functions is to be
able to use them in the middle of expressions, and if your function is
instead invoked in the middle of an expression such as

something = other_stuff + whatever*evaluate(xin)

then you’ve “lost” your opportunity to test for the error. It is both
too late in that you have already used the “bad” result, plus the
appropriate test on something is now non-trivial. That’s probably a
bigger caveat than the others mentioned above, but I also presume it is
obvious. In my view, if a function can’t be used in an expression like
this, then there isn’t much point in it being a function. But that’s
more my personal style choice than a real caveat.

2 Likes
46

It will be impossible to make everybody happy, but we’ve been discussing interfaces at length with @jeremie.vandenplas, @everythingfunctional, @hkvzjal, @gnikit et al. during our past Monthly calls.

So far we’ve come up with some interesting things:

  • All procedures that could be made pure were made pure (first time ever?)

  • Generalised KIND-dependent functions and 128-bit linear algebra support (first time ever?)

  • 2/3 interfaces for each operation, ideally:

    • subroutine interface that minimizes allocations
    • function interface returning optional error flag
    • .operator. interface with no error flag, pure if applicable

We also do not want to deviate too much from common interfaces such as NumPy and SciPy, because that’s what most people are used to nowadays.

See for example the current matrix inversion interface:

Of course all comments and ideas are very welcome and we’ll do our best to optimize across all the requests!

2 Likes
47

I fully agree with that. The other drawback of a return flags is that the function cannot be pure (or simple in F2023), which potentially prevent some compiler optimisations. In the current design choice purity is apparently reserved for the operator flavor, but not all operations can have an operator version.

48

Nice to hear about the progress you’re making with user-friendly linear algebra.

Concerning high-precision LAPACK and BLAS, there is at least one implementation from a Japanese group(GitHub - nakatamaho/mplapack: The MPLAPACK: multiple precision version of BLAS and LAPACK). I think I’ve seen other high-precision implementations before but can’t find the link now.

49

It looks like @FedericoPerini has the subroutine version, as requested. However, it is not pure. Why not?

50

Some of the high-level APIs have not been pureified yet because we’re still focussing on the BLAS/LAPACK backends, but the will come soon!

However, some functions in LAPACK contain intent(out) arguments, which prevents them, and all the upstream procedures using them, from being pure. So, some operations i.e. SVD cannot be made pure unless we change the LAPACK interface (which would break compatibility with the external libraries, or require large workarounds). One random example:

52

Yes, but it’s a separate issue. Maybe in some future this will be fixed, by changing the interfaces of these functions (and of the underlying implementation). In contrast, there will be no way to make the higher level functions pure in the future if they have an intent(out) return flag.

53

for those relying on pure BLAS/LAPACK procedures 2 out of 3 interfaces can be pure:

  • subroutine interface
  • .operator. interface
54

I’ve hoped for a long time that the standards committee would allow programmers to use any ASCII symbol they wanted in a user defined operator so I could mimic MATLAB and use a backslash () as an operator to solve linear systems. ie

interface operator(.\.)
   function solve (a, b) result(x)
   ..
   end function
end interface

x = a.\.b

I know it’s never going to happen but its something I would probably use a lot if available.

2 Likes
55

I’m not opposing your idea, but we could define a .div. operator and write x = b.div.a

56

Great idea! Let’s find one that would be suitable for LU solve in stdlib.
Some further ideas:

  • .backslash. (I would rather avoid .bs. … lol)
  • .ld. (left divide? like Matlab’s mldivide)
  • .inv. that would support both .inv.A meaning A^{-1} and A.inv.B meaning A^{-1}B
  • // (yes, we may overload the intrinsic concatenation, looks close to a division but may be confusing?)
57

My motivation for the backslash is I think trying to attract students whose only programming experience is MATLAB to Fortran would be easier if there was a way to make Fortran a little closer to what they already know from MATLAB. The interactive capabilities of LFORTRAN will go a long way in that direction. Enabling syntax that approximates MATLAB would also help.

4 Likes
58

I feel that inverse(...) would be less ambiguous, so that you could form, for example:

C=matmul(conjg(transpose(A)),inverse(B))
59

Here is how I see it… As mentioned above, one advantage of functions is that they can appear in expressions. Non-trivial expressions involving arrays can be tricky to optimize by the compilers (particularly to avoid creating too many temporaries…), and my understanding is that the pure or simple properties can help the compilers in these cases. You say that operators can do in this case: that’s correct, but style preferences set apart (I prefer functions to custom operators), not all functions will have an operator version I guess (returning the k first eigenvectors for instance require an additional argument)… maybe “most of them” is enough.

Having the return flags only in the “expert” subroutine versions would look fine to me.

How does Eigen manages errors?

60

Although the operations A.inv.B and matmul(inverse(A),B) are the same in abstract mathematics, they are very different in numerical analysis. Both the operation counts and the numerical properties differ.

61

Maybe, but it’s a bit ambiguous. If I saw A.inv.B I would expect that A is inverted and multiplied with B rather than being the solution to AX = B. Perhaps .div. would be better.

Also, are there any provisions for declaring the matrices as symmetric or Hermitian?

Finally, I think if the ‘info’ from LAPACK is not going to be returned, then if the matrix A is singular the entirety of B should be set to NaN. This way it fails safely.

(How do you include LaTeX formulae on these forums?)

1 Like
62

By enclosing them in a pair of dollar signs ($), looking at what others have done. The syntax is described at LaTeX - Mathematical Python.