SciPy: "there is interest in adopting this (PRIMA) implementation, but little appetite for taking on more Fortran code"

Thank you @conradoat for your comment and suggestion.

My goal is to make Professor Powell’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.

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.

9 Likes