I used Fortran 2008 a lot some years ago; initially I wanted to write a C++ wrapper around a Fortran library that seemed useful, but ended up liking the language better than C++ and just did the whole project in it.
I think a lot of people don’t realize that modern Fortran actually supports object-oriented programming and Matlab-style array slicing, and comes with built-in support for e.g. matrices and complex numbers. It even has some features inspired by functional languages, like being able to mark functions as pure (side-effect free). Combined with its general speed, this makes it a pretty nice language for numerical code.
Although for most new users, I think Julia might be a better option than Fortran these days.
This reads like the maintainer, mind you a maintainer of scipy, the premier foundational library in numerical python, does not even understand what PRIMA is or what modern fortran is like. I get people are busy, but this shakes my opinion of scipy a little bit. It doesn't take that long to read the readme... It does not reflect well on them.
I've been trying to debug some of scipy's statistics and it's... incomprehensible. I am not impressed. I'll happily take fortran over its growing morass of wierdo boost wrappers.
I thought that the PRIMA developer’s comment in this thread was really mature and thoughtful. Reproduced below (minus some formatting):
> My goal is to make Professor Powell 12’s solvers as accessible as possible to scientists, engineers, and algorithm researchers. I am not particularly in favor of or against any language. I hope that everyone can easily use Powell’s solvers in her/his favorite languages.
> The first implementation of PRIMA is in modern Fortran simply because Powell’s implementation was in Fortran 77. Using Fortran, I can systematically verify the bit-to-bit faithfulness of the modernized implementation (not only “faithful up to an epsilon”). In addition, the intrinsic support for matrix-vector calculations is a strong advantage when developing reference implementations (or templates) of numerical solvers — most numerical algorithms are combinations of such calculations anyway.
> The major motivation for developing the modern Fortran version is to provide a reference for the implementation in other languages, namely Python, MATLAB, C++, Julia, and R. A reference implementation must be structured, modularized, readable, understandable, and extendable. The original Fortran 77 code is a true masterpiece, but it is not proper at all for being used as a reference implementation. You do not want to use a spaghetti-style codebase with 244 GOTOs as a reference, or your implementation will be of the same style.
> Putting it more straightforwardly, I implemented the modern-Fortran version of PRIMA in order to develop versions that are entirely Fortran-free 3.
> Coming back to the point, I fully understand why the SciPy community has “little appetite for taking on more Fortran code”. The reputation of Fortran has been damaged over the years. I do not agree with the damaged reputation, and I feel sorry for those who do not have a chance to know (or refuse to know) modern Fortran due to this false reputation, but I do not blame them. It is the responsibility of the Fortran community to re-establish the reputation. I do not regard it as a bully to request non-Fortran implementations. It is not a question of surrendering or not.
> I do consider myself a member of the Fortran community. Taking my share of the aforementioned responsibilty, I will try to promote the usage of modern Fortran via the PRIMA project. There is nothing more convincing than a successful real-life project.
> For the inclusion of PRIMA in SciPy, I will keep communicating with both the SciPy and the Fortran community (e.g., those on this discourse), trying to find the best route. As pointed out by others, f2c is not an option due to its incapability of handling modern Fortran. Official C++ and Python implementations are being planned, but they will not be delivered in the near future. The most probable and practical solution, as suggested here and under the other thread, is to wrap the modern Fortran implementation of PRIMA using iso_c_bindings + ctypes or similar facilities. I hope the SciPy maintainers will accept this solution.
Yeah, I agree 100% - I really support the approach of the developer here and am totally aligned with all the reasoning.
The libprima/prima codebase is very readable, even if you're not accustomed to modern Fortran (let alone Fortran period). It has fantastic comments throughout the numerical algorithms, and even though there are a lot of lines of code, in most places it really seems to be a minimally-complex implementation with very little magic. I haven't built anything against the codebase yet but based on the examples I feel it will be far easier than many other libraries out there.
A motivated developer could very quickly port this code to their preferred numerical programming language.
Thank you so much for your very warm and encouraging comments!
PRIMA has been a black hole that absorbs all my time and energy in the past three years, which even puts my career (as a junor professor) in real danger. The positive feedback like yours is vital for me. Without it, I would not have the energe or courage to continue. Many thanks!
> A motivated developer could very quickly port this code to their preferred numerical programming language.
I am much glad to hear a person other than myself saying this. It is the very reason why I develop this reference implementation. PRIMA achieves its success if others can implement Powell's solvers to high quality using PRIMA as a referene, without the genius like Powell and without the long-term experience and strugling like me.
Hey, all the thanks go to you! I really wish we could do more to reward contributions like this when they're coming from folks in the academic world... the impact of this work is easily on par with anything else you could do as a junior faculty.
Have you considered setting up a regular search on GitHub to try to find public codebases that are likely using your library? That might be a great complementary set of statistics that you can use to promote the impact of this work.
I wonder why not just try to use more Pythonic solutions such as numba within scipy. I realize that maybe not in all cases can the same performance as Fortran be achieved but it strikes me as much more hackable for Python programmers than having to learn an new (old) language just for certain core algorithms. One of the principal reasons for switching to C implementations is because more people understand C than Fortran. This seems doubly true for Python-based jit solutions like numba or jax, so why not just skip C altogether and go for a more native port?
The last turtle is probably assembly, or at least intrinsics, in most good BLAS implementations. Probably called from C, because that’s easier in C.
The nice thing about Fortran is that you get something almost as good, but a normal engineer can write it. There are more C programmers than Fortran programmers, but there are more Fortran programmers than there are C programmers who can write really good assembly kernels. And all the C programmers in that last group are already working on vendor BLAS implementations.
That's a 6 year old documentation page, modern numpy does not use fortran in any way except that maybe you have a particular BLAS reference implementation.
Numpy is basically a wrapper to BLAS and LAPACK though, so it seems a bit odd to ignore the fact that those system libraries are usually written in Fortran.
The reference implementation is written in Fortran, but I'm not sure about the implementations actually used. OpenBLAS has fortran components, but I believe cuBLAS and MKL are vast majority C/C++.
This is a very fair challenge. My understanding was that the main ones in use were fortran and it is very interesting to see that may not be the case. Love to see any benchmark or other comparisons on the options, if you know of any good links.
cuBLAS is not Fortran because it targets NVIDIA GPUs. Despite what the grandparent says, MKL does include some Fortran in its codebase. Most other vendor-provided BLAS and LAPACK libraries are based on the reference Fortran implementation. Even OpenBLAS, the open-source reimplementation, has large parts written in Fortran.
Most BLAS and LAPACK implementations are written in Fortran. There are some C/C++ alternatives, but those are the exceptions (even if some are widely used in the FOSS community) and they still don’t benchmark as well even in 2023.
The good BLAS libraries are MKL and BLIS nowadays, I’m pretty sure.
BLIS is C and assembly if I remember correctly.
MKL is a big library, Wikipedia lists the languages used as C, C++, Fortran, and I think it is also pretty well known that they’ve sprinkled some magic assembly in there.
Intel hired Mr. Goto a while ago, he wrote gotoBLAS (from which openBlas is derived, so this one dude is responsible for a ton of FLOPs).
The thing that made gotoBLAS good was the hand-tuned assembly kernels, so we can be reasonably sure that MKL has hand tuned assembly kernels at this point.
I think the tongue-in-cheek comment is actually really good. It is a reminder that BLAS is more like a linear algebra API than a particular library.
Only someone at Intel could tell us for sure, but my understanding is that MKL is a C/C++ interface to a rather bog-standard BLAS/LAPACK library written in Fortran, but with various inner loops replaced or rewritten at runtime with a selection of assembly-optimized x86 SIMD codes chosen based on the capabilities of the CPU.
Reference BLAS is of course Fortran but nobody would actually use it. Some other good BLAS are (in order of historical to cutting edge), ATLAS, OpenBLAS, and then BLIS. They all mostly use C and assembly.
BLIS makes a good BLAS. It almost keeps up with MKL, according to the nice plots on their GitHub.
They use C mostly, because the magic is in the GEMM (and possibly TRSM) kernels, which you are going to write in assembly anyway.
Intel also released oneDNN open source, which has an awful lot of BLAS-like functionality. No Fortran in there.
MKL is older than all these libraries of course (maybe not ATLAS, I’m not sure actually), and has a lot of other functionality (sparse stuff, FFTs), so I’m sure it has more Fortran! But it would be atypical compared to the rest of the BLAS implementations, if, like, their DGEMM was Fortran.
> While a FORTRAN 77 compiler is not necessary for building NumPy, it is needed to run the numpy.f2py tests. These tests are skipped if the compiler is not auto-detected.
Numpy is most importantly BLAS and LAPACK, which can be implemented in any language. The best CPU version is probably still MKL, which is a mix of Fortran, C, and assembly. C is used for these kinds of projects because it is easy to sprinkle assembly into C.
Fortran solves a particular problem: produce a pretty good binary from readable code that a normal PHD student or engineer can write. This is not really the problem that a tuned BLAS/LAPACK library has to face, vendors throw very experienced folks at these libraries to chase after those last 5% performance boosts that Fortran leaves on the table.
A native port is indeed planned. However, since we are talking about a project of about 10K lines of code (though many of the lines are comments), such a port will not be delivered very soon.
In fact, native implementations of PRIMA in Python, MATLAB, C++, Julia, and R will all be done in the future. See https://github.com/libprima/prima#other-languages . But they will take time. PRIMA has been a one-man project since it started three years ago. Community help is greatly needed.
This is crazy. Get over yourself! Whether you like it or not, Fortran is still THE industry standard for scientific codes, and the fastest compiled language on many architectures.
In the defense of the SciPy response, while there may be some programming language zealotry behind the position, it seems to be from a lack of expertise in SciPy to maintain its Fortran code that's pushing SciPy to avoid new Fortran additions. Hopefully OP gets PRIMA through! Been admiring this effort, Dr. Zhang.
F2py, which is already a standard part of numpy, does the heaving lifting and generates idiomatic python wrappers for Fortan code.
It's just zealotry of them going "eeew fortran that's like really really old". Fortran got multidimensional dense numeric arrays right. Modern Fortran allows you to do elemental/vectorized operations on them. C does not even come close, and requires tons of preprocessor macros if you try. Every other 'modern' language introduces a bunch of tooling and dependencies just to get where Fortran is, because they just don't have a focus on scientific computing.
> It's just zealotry of them going "eeew fortran that's like really really old".
The developers of SciPy are maintaining a codebase that is 18% fortran, according to [1]. The only language that makes up a larger part of the codebase is python.
Don't you think it's possible their opinion on fortran is actually informed by quite a bit of experience?
> The central issue is that their experience with the old Fortran doesn't translate to modern Fortran; they are two different beasts.
This is so true! I hope anyone could spend a few minutes checking what modern Fortran is before blaming it. A good reference is https://fortran-lang.org/en/
I guess one of the major issues is the compilation and packaging of the Fortran code. It is also my experience when developing PDFO ( https://www.pdfo.net/ ), the predecessor of PRIMA.
As long as PRIMA is maintained (which seems likely independently), it's only the bindings that need to be maintained. And it's trivial to write a Fortran interface in C or C++ (perhaps that can be included in PRIMA anyway) and then SciPy can pretend like it's not Fortran...
Looking at the original issue, scipy maintainers don't want to use LGPL code (why?!?) but they use code that appears in Numerical Recipes (which has super onerous licensing)?!?
First, if you read through that scipy issue (https://github.com/scipy/scipy/issues/18118 ) the author was willing and able to relicense PRIMA under a 3-clause BSD license which is perfectly acceptable for scipy.
For the numerical recipes reference, there is a mention that scipy uses a slightly improved version of Powell's algorithm that is originally due to Forman Acton and presumably published in his popular book on numerical analysis, and that also happens to be described & included in numerical recipes. That is, unless the code scipy uses is copied from numerical recipes, which I presume it isn't, NR having the same algorithm doesn't mean that every other independent implementation of that algorithm falls under NR copyright.
thanks, I misunderstood the NR reference. (In my field... it is common for people to copy and paste terrible NR code without understanding how bad it is or the licensing terms...).
Any resources for learning Fortran in 202X that don't bother going over the old cruft?
Last time I looked (admittedly, years ago), most resources seemed to be in the vein of "here's Fortran 1990, and once you've learned that here are a few chapters on F03 and F08, treated separately".
It'd be great to learn there exists a "Fortran: the good parts" sort of thing that just completely elides the outdated stuff.
where you can find all the gurus of modern Fortran. They are very willing to help if you post a question. I have asked all my Fortran questions over there during the development of PRIMA.
Even though it's not great but I understand why this maybe the direction SciPy maintainers want to go in. They have limited capacity and I think focusing on a single language makes sense especially if there are way more developers in that language, C.
One way to solve this would be of course, to have more people from the Fortran community being part of the SciPy community. But I don't know how big the Fortran community really is to be able to do that.
There just isn't much interest. HPC focuses on distributed and product doesn't care about open source. That mostly leaves just professors and oddballs.
Interesting response. I develop the Julia SciML organization https://sciml.ai/ and we'd be more than happy to work with you to get wrappers for PRIMA into Optimization.jl's general interface (https://docs.sciml.ai/Optimization/stable/). Please get in touch and we can figure out how to set this all up. I personally would be curious to try this out and do some benchmarks against nlopt methods.
72 comments
[ 2.8 ms ] story [ 129 ms ] threadOptimization Without Derivatives: PRIMA Fortran Version and Inclusion in SciPy, https://news.ycombinator.com/item?id=35959991
SciPy enhancement: The Fortran 77 implementation of COBYLA is buggy and challenging to maintain. Switch to the PRIMA implementation? #18118, https://github.com/scipy/scipy/issues/18118#issuecomment-155...
You only need Python's ctypes module to load the compiled Fortran library thanks to Fortran's iso_c_binding module.
Using f2c as suggested seems like a worse alternative, not a solution.
essentially madness.
Modern fortran is quite nice, much nicer than C from a Python programmer's perspective.
Totally agree. I hope everyone making comments on Fortran can spend a few minitues reading what modern Fortran is, e.g., by checking
https://fortran-lang.org/en/
I think a lot of people don’t realize that modern Fortran actually supports object-oriented programming and Matlab-style array slicing, and comes with built-in support for e.g. matrices and complex numbers. It even has some features inspired by functional languages, like being able to mark functions as pure (side-effect free). Combined with its general speed, this makes it a pretty nice language for numerical code.
Although for most new users, I think Julia might be a better option than Fortran these days.
> My goal is to make Professor Powell 12’s solvers as accessible as possible to scientists, engineers, and algorithm researchers. I am not particularly in favor of or against any language. I hope that everyone can easily use Powell’s solvers in her/his favorite languages.
> The first implementation of PRIMA is in modern Fortran simply because Powell’s implementation was in Fortran 77. Using Fortran, I can systematically verify the bit-to-bit faithfulness of the modernized implementation (not only “faithful up to an epsilon”). In addition, the intrinsic support for matrix-vector calculations is a strong advantage when developing reference implementations (or templates) of numerical solvers — most numerical algorithms are combinations of such calculations anyway.
> The major motivation for developing the modern Fortran version is to provide a reference for the implementation in other languages, namely Python, MATLAB, C++, Julia, and R. A reference implementation must be structured, modularized, readable, understandable, and extendable. The original Fortran 77 code is a true masterpiece, but it is not proper at all for being used as a reference implementation. You do not want to use a spaghetti-style codebase with 244 GOTOs as a reference, or your implementation will be of the same style.
> Putting it more straightforwardly, I implemented the modern-Fortran version of PRIMA in order to develop versions that are entirely Fortran-free 3.
> Coming back to the point, I fully understand why the SciPy community has “little appetite for taking on more Fortran code”. The reputation of Fortran has been damaged over the years. I do not agree with the damaged reputation, and I feel sorry for those who do not have a chance to know (or refuse to know) modern Fortran due to this false reputation, but I do not blame them. It is the responsibility of the Fortran community to re-establish the reputation. I do not regard it as a bully to request non-Fortran implementations. It is not a question of surrendering or not.
> I do consider myself a member of the Fortran community. Taking my share of the aforementioned responsibilty, I will try to promote the usage of modern Fortran via the PRIMA project. There is nothing more convincing than a successful real-life project.
> For the inclusion of PRIMA in SciPy, I will keep communicating with both the SciPy and the Fortran community (e.g., those on this discourse), trying to find the best route. As pointed out by others, f2c is not an option due to its incapability of handling modern Fortran. Official C++ and Python implementations are being planned, but they will not be delivered in the near future. The most probable and practical solution, as suggested here and under the other thread, is to wrap the modern Fortran implementation of PRIMA using iso_c_bindings + ctypes or similar facilities. I hope the SciPy maintainers will accept this solution.
The libprima/prima codebase is very readable, even if you're not accustomed to modern Fortran (let alone Fortran period). It has fantastic comments throughout the numerical algorithms, and even though there are a lot of lines of code, in most places it really seems to be a minimally-complex implementation with very little magic. I haven't built anything against the codebase yet but based on the examples I feel it will be far easier than many other libraries out there.
A motivated developer could very quickly port this code to their preferred numerical programming language.
PRIMA has been a black hole that absorbs all my time and energy in the past three years, which even puts my career (as a junor professor) in real danger. The positive feedback like yours is vital for me. Without it, I would not have the energe or courage to continue. Many thanks!
> A motivated developer could very quickly port this code to their preferred numerical programming language.
I am much glad to hear a person other than myself saying this. It is the very reason why I develop this reference implementation. PRIMA achieves its success if others can implement Powell's solvers to high quality using PRIMA as a referene, without the genius like Powell and without the long-term experience and strugling like me.
Many thanks! --- Zaikun
Have you considered setting up a regular search on GitHub to try to find public codebases that are likely using your library? That might be a great complementary set of statistics that you can use to promote the impact of this work.
> The core of NumPy is well-optimized C code.
https://numpy.org/doc/1.13/user/building.html
The nice thing about Fortran is that you get something almost as good, but a normal engineer can write it. There are more C programmers than Fortran programmers, but there are more Fortran programmers than there are C programmers who can write really good assembly kernels. And all the C programmers in that last group are already working on vendor BLAS implementations.
[0]: https://github.com/search?q=repo%3Anumpy%2Fnumpy+path%3A**%2...
BLIS is C and assembly if I remember correctly.
MKL is a big library, Wikipedia lists the languages used as C, C++, Fortran, and I think it is also pretty well known that they’ve sprinkled some magic assembly in there.
The reply about cuBLAS was tongue in cheek as it is a GPU BLAS library
The thing that made gotoBLAS good was the hand-tuned assembly kernels, so we can be reasonably sure that MKL has hand tuned assembly kernels at this point.
I think the tongue-in-cheek comment is actually really good. It is a reminder that BLAS is more like a linear algebra API than a particular library.
Reference BLAS is of course Fortran but nobody would actually use it. Some other good BLAS are (in order of historical to cutting edge), ATLAS, OpenBLAS, and then BLIS. They all mostly use C and assembly.
BLIS makes a good BLAS. It almost keeps up with MKL, according to the nice plots on their GitHub.
https://github.com/flame/blis/blob/master/docs/Performance.m...
They use C mostly, because the magic is in the GEMM (and possibly TRSM) kernels, which you are going to write in assembly anyway.
Intel also released oneDNN open source, which has an awful lot of BLAS-like functionality. No Fortran in there.
MKL is older than all these libraries of course (maybe not ATLAS, I’m not sure actually), and has a lot of other functionality (sparse stuff, FFTs), so I’m sure it has more Fortran! But it would be atypical compared to the rest of the BLAS implementations, if, like, their DGEMM was Fortran.
Fortran solves a particular problem: produce a pretty good binary from readable code that a normal PHD student or engineer can write. This is not really the problem that a tuned BLAS/LAPACK library has to face, vendors throw very experienced folks at these libraries to chase after those last 5% performance boosts that Fortran leaves on the table.
In fact, native implementations of PRIMA in Python, MATLAB, C++, Julia, and R will all be done in the future. See https://github.com/libprima/prima#other-languages . But they will take time. PRIMA has been a one-man project since it started three years ago. Community help is greatly needed.
Thanks.
It's just zealotry of them going "eeew fortran that's like really really old". Fortran got multidimensional dense numeric arrays right. Modern Fortran allows you to do elemental/vectorized operations on them. C does not even come close, and requires tons of preprocessor macros if you try. Every other 'modern' language introduces a bunch of tooling and dependencies just to get where Fortran is, because they just don't have a focus on scientific computing.
The developers of SciPy are maintaining a codebase that is 18% fortran, according to [1]. The only language that makes up a larger part of the codebase is python.
Don't you think it's possible their opinion on fortran is actually informed by quite a bit of experience?
[1] https://github.com/scipy/scipy/
The central issue is that their experience with the old Fortran doesn't translate to modern Fortran; they are two different beasts.
This is so true! I hope anyone could spend a few minutes checking what modern Fortran is before blaming it. A good reference is https://fortran-lang.org/en/
This is also reflected by the following comment at https://github.com/scipy/scipy/issues/18118#issuecomment-155...
> The real maintenance burden is the compilation and the packaging not the code itself that we are carrying around.
For the numerical recipes reference, there is a mention that scipy uses a slightly improved version of Powell's algorithm that is originally due to Forman Acton and presumably published in his popular book on numerical analysis, and that also happens to be described & included in numerical recipes. That is, unless the code scipy uses is copied from numerical recipes, which I presume it isn't, NR having the same algorithm doesn't mean that every other independent implementation of that algorithm falls under NR copyright.
Modern Fortran is quite different from Fortran 77, while being as powerful, if not more.
In addition, there has been a significant community effort on improving and modernising the legacy packages, the ecosystem, and the language itself.
With projects like LFortran (https://lfortran.org/), fpm (https://github.com/fortran-lang/fpm), and stdlib (https://github.com/fortran-lang/stdlib), I believe that Fortran will enjoy prosperity again.
Last time I looked (admittedly, years ago), most resources seemed to be in the vein of "here's Fortran 1990, and once you've learned that here are a few chapters on F03 and F08, treated separately".
It'd be great to learn there exists a "Fortran: the good parts" sort of thing that just completely elides the outdated stuff.
Also https://fortran-lang.org/en/ is a good hub to find more information, packages, etc.
https://fortran-lang.org/en/learn/
BTW, the Fortran Lang website (https://fortran-lang.org/) is THE official website of modern Fortran, maintained by the modern Fortran community.
You may also take a look at the modern Fortran discourse at
https://fortran-lang.discourse.group/
where you can find all the gurus of modern Fortran. They are very willing to help if you post a question. I have asked all my Fortran questions over there during the development of PRIMA.
One way to solve this would be of course, to have more people from the Fortran community being part of the SciPy community. But I don't know how big the Fortran community really is to be able to do that.