Degenerate Conic

Polynomial Roots

Consider the problem of finding all the roots of a polynomial. This involves solving the following equation:

This problem occurs in many different applications, such as computing atmospheric density in a spacecraft orbit simulation. In general, both the coefficients and roots can be complex numbers. A special case involves only real coefficients. General methods exist for any degree polynomial, and custom methods also exist for polynomials of specific degree. The quadratic formula is the most famous analytical solution for polynomials of degree 2. In general, numerical methods are necessary to solve polynomials of degree 5 and higher.

Newton's method is frequently used, which requires an initial guess and an iterative procedure. Another elegant way to solve for all the roots at once is to find the eigenvalues of the companion matrix.

Fortran library

My polyroots-fortran library includes various methods to solve this problem. Many of the implementations have legacy going back decades (for example, from SLATEC). It's easy to incorporate into other codes using the Fortran Package Manager. The available methods in the library are:

Example

As an example, let's generate a plot of all the roots for all \(11^{th}\) degree polynomials with coefficients either +1 or -1. That is:

We'll use the polyroots method, which solves for the roots of a real polynomial equation by computing the eigenvalues of the companion matrix (using the LAPACK eigen solver DGEEV). To get all the coefficient permutations, we will use a recursive algorithm, and pyplot-fortran to generate the plot. Here is the code:

program main

use polyroots_module, only: polyroots, wp => polyroots_module_rk
use pyplot_module,    only: pyplot

implicit none

integer,parameter :: degree = 11 !! polynomial degree
integer,dimension(2) :: icoeffs = [-1,1] !! set of coefficients
integer,dimension(degree+1) :: a !! coefficients of polynomial
type(pyplot) :: plt  !! for making the plot
integer      :: ierr !! error code from [[polyroots]]
integer      :: j    !! counter

call plt%initialize(grid=.true.,xlabel='$\\Re(z)$',ylabel='$\\Im(z)$',&
                    title='Degree 11 Polynomial Roots',usetex=.true.,&
                    font_size=30, axes_labelsize=30, &
                    xtick_labelsize=30, ytick_labelsize=30, &
                    figsize=[20,10])
call generate(1)
call plt%showfig()

contains

    recursive subroutine generate(i)
    integer, intent(in) :: i
    integer :: ix
    real(wp),dimension(degree) :: zr, zi !! roots
    if (i > degree+1) then
        ! solve for this root and add it to the plot:
        call polyroots(degree,real(a,wp),zr,zi,ierr)
        call plt%add_plot(zr,zi,label='',linestyle='bo',markersize=2)
    else
        do ix = 1,size(icoeffs)
            a(i) = icoeffs(ix)
            call generate(i+1)
        end do
    end if
    end subroutine generate

end program main

The first set of roots, for \(a = [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]\), are: [\(-1\), \(\pm i\), \(\pm \sqrt3/2 \pm 0.5 i\), \(\pm 0.5 \pm \sqrt3/2\)].

The resultant plot with all the roots is shown here:

References

• A. Edelman and H. Murakami, "Polynomial roots from companion matrix eigenvalues", Math. Comp. 64 (1995), 763-776
• C. Moler, "ROOTS - Of Polynomials, That Is", MathWorks, 1991
• Abel–Ruffini theorem [wikipedia]
• Less verbose combinations, fortran-lang Discourse, Nov 2022

Converting Legacy C++ Code To Modern Fortran

"I am bound Upon a wheel of fire, that mine own tears Do scald like molten lead." - Shakespeare, "King Lear", Act IV, Scene 7

Typically, we see examples of poor saps who think that legacy Fortran codes must be converted to "modern" languages such as C++. Well, let's try to do the reverse! Let's try converting some "legacy" C++ code to modern Fortran. The example I will show is the Jacchia-Roberts atmosphere model from GMAT (General Mission Analysis Tool). This model is widely used for satellite orbit simulations to model satellite drag. There is no modern Fortran implementation publicly available as far as I can tell, so let's create one. And just to see what happens, I'll do this with assistance from AI, particularly Claude Sonnet (initially 4.5, but also trying 4.6) in GitHub Copilot. What could go wrong?

All About the Vibes

First, we start with the C++ files we want to convert. These files are at least 20 years old, with heritage going back even further. The main ones are:

• JacchiaRobertsAtmosphere.cpp - the main implementation of the Jacchia-Roberts density model.
• SolarFluxReader.cpp - the algorithm to read in a space weather file. There are a couple different file formats supported here, but I'm only going to worry about the CSSI one.
• A1Mjd.cpp - Some time conversion routines. This will be important later when trying to diagnose some validation issues.

Using AI is a very interesting and also frustrating experience. It's like working with someone who is a genius, but is also an idiot and sometimes a compulsive liar. The initial conversion got me 80% of the way there but was subtly wrong in various ways. I started with a basic prompt ("Convert this C++ file into modern Fortran") and went from there. There was a lot of back and forth and manual interventions to fix little issues (syntax errors, etc.) In one instance, I tried to get it to identify the source of a bug, and it "thought" for 15 minutes until I finally just gave up and looked at the code myself and fixed it instantly. Sometimes, it would get fixated on one approach that wasn't correct, and I'd have to keep telling it to ignore that approach and try something else (e.g., I had to figure out the correct way to do something, and once I gave it more directions, it was fine and could do it). One thing that was annoying was that it kept trying to generate code from first principles rather than performing a straight up conversion of the C++ code to Fortran. I had to keep reminding it that what I wanted was a faithful translation. At one point, I switched to Sonnet 4.6 and it also found a couple of code typos that Sonnet 4.5 had generated.

The Result

See jacchia-roberts-fortran for the end result. The library provides a class that can be used to initialize the model with a given space weather file (or constant values can be specified). Then a method is provided for computing the density. An example use case is:

program example

use jacchia_roberts_module, only: jacchia_roberts_type
use jacchia_roberts_kinds,  only: ip, dp

implicit none

type(jacchia_roberts_type) :: jr
integer(ip) :: status
real(dp) :: density

! example inputs:
real(dp),parameter :: rad_earth = 6356.766_dp ! Earth polar radius (km)
character(len=*),parameter :: sw_file = 'data/SpaceWeather-All-v1.2.txt' ! space weather file to load
real(dp),parameter :: utc_mjd = 59215.5_dp ! MJD: Jan 1, 2021 12:00 UTC
real(dp),parameter :: alt_km = 200.0_dp ! geodetic altitude (km)
real(dp),parameter :: lat = 20.0_dp ! geodetic latitude (deg)
real(dp),dimension(3),parameter :: r = [7000.0_dp, 0.0_dp, 0.0_dp] ! spacecraft position vector (km)
real(dp),dimension(3),parameter :: sun = [1.0_dp, 0.0_dp, 0.0_dp] ! sun direction unit vector

! initialize the model:
call jr%initialize(rad_earth, filename=sw_file, status=status)

! compute the density:
density = jr%density(alt_km, r, sun_vector, lat, utc_mjd)

end program example

The idea is, first you would initialize the class, and then during your simulation, call the model to get the atmospheric density using the current epoch and state of the spacecraft. A flow chart of the process would look a little like this:

The resultant Fortran code is very clean compared to the C++ code. Consider this example (the rho_cor function). The original C++ code:

Real JacchiaRobertsAtmosphere::rho_cor(Real height, Real a1_time, Real geo_lat,
                                       GEOPARMS *geo)
{
   Real geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha;
   Real sin_lat, eta_lat;

   // Compute geomagnetic activity correction
   if (height < 200.0)
   {
      geo_cor = 0.012 * geo->tkp + 0.000012 * exp(geo->tkp);
   }
   else
   {
      geo_cor = 0.0;
   }

   // Compute semiannual variation correction
   f = (5.876e-7 * pow(height, 2.331) + 0.06328) *
      exp(-0.002868 * height);
   day_58 = (a1_time - 6204.5)/365.2422; // @todo - should this use GmatConstants value?
   tausa = day_58 + 0.09544 * (
         pow( 0.5*(1.0 + sin(2*GmatMathConstants::PI*day_58 +6.035)), 1.65 )  - 0.5);
   alpha = sin(4.0*GmatMathConstants::PI*tausa + 4.259);
   g = 0.02835 + (0.3817 + 0.17829 * sin(2*GmatMathConstants::PI*tausa + 4.137)) *
         alpha;
   semian_cor = f * g;

   // Compute seasonal latitudinal variation
   sin_lat = sin(geo_lat);
   eta_lat = sin(2.0*GmatMathConstants::PI*day_58 + 1.72) * sin_lat * fabs(sin_lat);
   slat_cor = 0.014 * (height - 90.0) * eta_lat *
      exp(-0.0013 * (height - 90.0) * (height - 90.0));

   return pow(10.0, geo_cor + semian_cor + slat_cor);
}

And the Fortran equivalent:

pure function rho_cor(height, utc_mjd, geo_lat, geo) result(correction)
    real(dp), intent(in) :: height !! Spacecraft altitude (km)
    real(dp), intent(in) :: utc_mjd !! UTC Modified Julian Date
    real(dp), intent(in) :: geo_lat !! Geodetic latitude of spacecraft (radians)
    type(geoparms_type), intent(in) :: geo !! Geomagnetic parameters
    real(dp) :: correction !! Correction factor (multiplicative)

    real(dp) :: geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha, sin_lat, eta_lat

    ! Compute geomagnetic activity correction
    if (height < 200.0_dp) then
        geo_cor = 0.012_dp * geo%tkp + 0.000012_dp * exp(geo%tkp)
    else
        geo_cor = 0.0_dp
    end if

    ! Compute semiannual variation correction
    f      = (5.876e-7_dp * height**2.331_dp + 0.06328_dp) * exp(-0.002868_dp * height)
    day_58 = (utc_mjd - 36204.0_dp) / 365.2422_dp ! years since Jan 1, 1958 midnight
    tausa  = day_58 + 0.09544_dp * &
            ((0.5_dp * (1.0_dp + sin(2.0_dp * PI * day_58 + 6.035_dp)))**1.65_dp - 0.5_dp)
    alpha  = sin(4.0_dp * PI * tausa + 4.259_dp)
    g      = 0.02835_dp + (0.3817_dp + 0.17829_dp * sin(2.0_dp * PI * tausa + 4.137_dp)) * alpha
    semian_cor = f * g

    ! Compute seasonal latitudinal variation
    sin_lat  = sin(geo_lat)
    eta_lat  = sin(2.0_dp * PI * day_58 + 1.72_dp) * sin_lat * abs(sin_lat)
    slat_cor = 0.014_dp * (height - 90.0_dp) * eta_lat * &
                exp(-0.0013_dp * (height - 90.0_dp)**2)

    correction = 10.0_dp**(geo_cor + semian_cor + slat_cor)

end function rho_cor

The main differences in the Fortran are:

• The real kinds use the dp kind variable, allowing us to easily compile the code with a specified real kind (e.g., single, double, or quad precision). This is one of the neat features of Fortran, and how I write all my codes.
• The crazy GSFC MJD (a1_time) is replaced with real MJD, so we had to modify that one constant (the epoch of Jan 1, 1958).
• Normal syntactical differences: abs instead of fabs, no semicolons, ** operator instead of the pow function, then/endif instead of {} symbols, etc. The C++ code using (height - 90.0) * (height - 90.0) instead of (height - 90.0_dp)**2 is also kind of funny.

I find the Fortran version easier to read without the extra C++ clutter, but maybe it's just because i'm more used to it? This aspect of Fortran (it being easy to read) isn't always appreciated or discussed in these "language war" type discussions, but it's very real in my opinion, especially for code is that mostly math.

Here's an example output density plot (generated using my pyplot-fortran library) for three epochs near the max/min/average parts of the solar cycle:

Validation

One of the most interesting parts of the process was using AI to help validate the result. I had no unit tests to start with, but I had two sources of "truth":

• GMAT itself, which I could run to generate density data for different altitudes, epochs, etc. and compare it to the Fortran results.
• An older Fortran implementation of Jacchia-Roberts. This was a Fortran 77 implementation by Valdemir Carrara (Instituto Nacional de Pesquisas Espaciais, INPE) from the 80s, which I obtained from his website, which is no longer there (I archived the codes on GitHub in 2019 thinking I would do something with them one day, and now that day has arrived). This was also another AI opportunity, since the space weather component of this code was not preserved, so I had Copilot generate a wrapper so that it could call the modern Fortran module we had produced instead.

AI seems to get "unit test happy" sometimes, spewing out numerous long unit tests that do something simple and then prints up elaborate tables with unicode icons and paragraphs of text that tells you what it just tested. That seemed a bit excessive, and I ended up deleting many of the tests it generated, but kept a good number of them. Overall it was quite useful for helping to generate and evaluate tests. When tests failed (e.g., the new results didn't match the other two data sets) I would tell it to investigate to locate the source of the discrepancy, which it could sometimes do successfully, and sometimes not. The AI went on a flight of fancy for a while trying to figure out a timing discrepancy. It did help me debug it, but it couldn't quite figure it out. At one point I asked if it could go to GitHub and look at the full GMAT code that was there, it told me it could and provided an unconvincing reply suspiciously quickly that it had done so. I guess it was lying to me (Copilot I thought we were friends). I ended up downloading the extra file manually, putting it in the workspace, and directing the AI's attention to it, which seemed to help.

During the validation process, the AI did help me identify three possible issues with the original C++ code. I link the three bug reports I submitted here:

1. PrepareKpData uses wrong row index for F10.7a - This is one that the AI detected that I would never have noticed in a million years. It looks like a copy/paste type bug but I'm not 100% sure if it's a real issue or not.
2. Legacy UTC calculation used by JacchiaRobertsAtmosphere - This is one that I eventually noticed while trying to debug timing discrepancies during the validation process comparing the results from GMAT to the new version. GMAT (and other GSFC tools) use this crazy alternate definition of modified Julian date (MJD). So while trying to figure out how to convert what it was outputting in a real MJD, I realized they were using an approximation of UTC (it turns out the correct UTC value output from GMAT is slightly different from what is being used internally by the atmosphere model). It looks like something that was a holdover from the earlier tools, but I'm not sure if there's a good reason to keep using it. According to the git history when this file was first committed, the file was originally created for the GSS project in 1995 (I presume this is referring to the Generalized Support Software project at GSFC).
3. Low-altitude Jacchia-Roberts bugs - This one is very interesting. There appear to be bugs in the GMAT code for low (<125 km) altitudes. Note that this part of the code is there but not used, since an exception is raised for this case, even though the Jacchia model should support the low altitude ranges (in Roberts' paper, he says his results "are identical" to Jacchia's original equations in the 90 - 125 km altitude regime). According to the comments and git history on SourceForge, this file was committed in 2004 and the source was the Windows Swingby code (Swingby was another legacy NASA tool that was one of the progenitors of STK/Astrogator). The roots function comments says it was converted from existing Fortran code (likely from GTDS, which included the Jacchia-Roberts model). However, the polynomial deflation technique was written in 1993 based on the "Numerical Recipes in C" book (note that Numerical Recipes talks about the need to "polish" a root found by this method, something not implemented in this code. Maybe that's where they went wrong?) Is it possible that Swingby had this same bug all those years ago? Or maybe my Fortran conversion is the buggy one? In any event, I replaced this algorithm in the Fortran version with an alternate root finding approach (based on the one in the INPE code) which seems to work fine. I also tested routines from polyroots-fortran and those gave the same results (I might end up just replacing this custom routine with the external dependency).

Legacy Schmegacy

The earlier Fortran tool that some of this code was based on was the Goddard Trajectory Determination System (GTDS). This tool, like many Fortran tools, had a long heritage, but was thrown in the trash in the end, never modernized but replaced with C++ rewrites. Incidentally, GTDS and Swingby were developed by Computer Sciences Corporation (CSC), a company that existed from 1959 to 2017, which was cofounded by Roy Nutt, one of the members of the original team that created Fortran. So, now we have come full circle on this wheel of fire with a new modern Fortran port.

Conclusions

The end result of our AI-assisted C++ to Fortran conversion produced something that I think is likely quite useable. Overall I believe it did save me time compared to how long it would have taken me to do manually. But, getting this over the finish line definitely required manual intervention and domain knowledge. The process found a few potential bugs in the original code that need to be looked at by actual experts, since I don't trust the AI not a gaslight me with plausible-sounding reasons for them. Having the AI write the unit tests was also a huge time saver and help with debugging the implementation. In the end, for a task of this magnitude, I think AI can get you something like 95% of the way there, but then that last 5% needs to be completed and checked very carefully by somebody who more or less knows what they are doing. Blindly "vibe coding" is dangerous for technical applications such as this, since it's very easy for an AI to provide plausible looking procedures that are totally (or even worse, subtly) wrong.

Note: No AI was used to write the text of this article.

References

1. Jacchia, L. G., "New Static Models of the Thermosphere and Exosphere with Empirical Temperature Profiles", Smithsonian Astrophysical Observatory Special Report No. 313, 1970.
2. Roberts, C. E., Jr., "An Analytic Model for Upper Atmosphere Densities Based Upon Jacchia's 1970 Models", Celestial Mechanics, Vol. 4, pp. 368-377, 1971.
3. GMAT: General Mission Analysis Tool
4. Archive of legacy INPE atmosphere models github/@jacobwilliams
5. Atmosphere Models [degenerateconic.com] 2021-10-24
6. "Does the Programming Language Still Matter When AI Writes Your Code? I Tried Python and Came Back to C#" medium.com/@binnmti, Mar 16, 2026.
7. J. Carrico & E. Fletcher, "Software Architecture and Use of Satellite Tool Kit’s Astrogator Module for Libration Point Orbit Missions", Libration Point Orbits and Applications, pp. 471-487 (2003)
8. A. C. Long, et al., "Goddard Trajectory Determination System (GTDS) Mathematical Theory Revision 1", FDD/552-89/001, Computer Sciences Corporation, July 1989.

Fortran Compiler Updates (2026)

Fortran compiler development continues, with regular updates across a range of compilers. Some of the recent updates are summarized below:

Intel

Intel has announced an update to their Fortran compiler (ifx), which is now at version 2026.0. ifx replaced the venerable ifort compiler a few years ago, and fully supports the Fortran 2018 standard and parts of the Fortran 2023 standard. The Fortran compiler is bundled as part of the Intel's oneAPI Toolkit. According to the release notes, there are only a few updates in this release:

• Support for the leading zero edit mode, which allows a user to control the output of optional leading zeros in numeric outputs.
• A coarray update related to NOTIFY_TYPE, which I don't understand. Maybe this is useful for somebody. Coarrays are the built-in MPI type feature of Fortran.
• Support for C-style (a ? b : c) conditional expressions. This was a feature added to Fortran 2023, which is a profoundly weird fit for Fortran syntax, but it does allow a one-line short-circuited logical expression that wasn't otherwise possible before (the merge statement was the closest we had for that before but isn't quite the same).
• Various new OpenMP 6.0 features and a new -qopenmp-threadprivate option.
• On Windows, support for Microsoft Visual Studio 2026.

Presumably bugs have also been fixed, since the road to ifx adoption has been a little bumpy, with all new bugs that weren't present in ifort. Intel used to call out what bugs were fixed in their compilers, but they haven't done that in years, and users are left to figure out for themselves what might have been fixed.

GFortran

The GNU Project also just released GCC 16.1. According to the release notes the changes in GFortran are:

• Coarrays using native shared memory multithreading on single node machines and handling Fortran 2018's TEAM feature.
• Support the Fortran 2018 extensions to the import statement, the reduce intrinsic and the new generic statement.
• The Fortran 2023 additions to the trigonometric functions are now supported (such as the sinpi intrinsic).
• The Fortran 2023 split intrinsic subroutine is now supported and c_f_pointer now accepts an optional lower bound as an argument.
• The -fexternal-blas64 option has been added to call external BLAS routines with 64-bit integer arguments for matmul.
• Parameterized Derived Types (PDT) support is improved. GFortran's implementation of PDT's has been buggy and incomplete for years. PDT was a feature added in Fortran 2003 that looks like it might be useful for writing generic code, but in practice it is not very useful. Fortran 202y, the next major revision of the Fortran standard is slated to include better generics capabilities.

Honorable mentions

• LFortran also continues to make steady progress towards its first beta release.
• LLVM Flang development also continues, the latest release was 22.1.0. We still don't have Mac builds of this compiler available through conda, so I've never tested this one before. I've gotten this far without having to manually compile a Fortran compiler, and I'm not about to start now.

See also

• Intel® Fortran Compiler 2026.0 is available, Intel Fortran Compiler Forum, 05-04-2026
• Gcc/gfortran 16 released Fortran-lang forum, 2026-Apr-23
• GCC 16 Release Series Changes, New Features, and Fixes
• LFortran Blog
• LLVM Download Page
• Fortran Compilers [fortran-lang.org]

Root Finding

Finding the \(x\) where a scalar function \(f(x) = 0\) is a common problem in numerical programming. There are numerous applications to this. In orbital mechanics, for example, we may want to find the time at which a spacecraft reaches a certain altitude, or enters the Earth's shadow. For these cases, the independent variable is time, and the function \(f\) requires integration of the equations of motion of the spacecraft (for example, by using an IVP solver such as DDEABM). Integrating to an event is often referred to as a "G-stop" feature in these types of applications.

One way to solve this problem is using a derivative-free bracketing method. For these, only function evaluations are used, and the algorithm begins within a set of bounds (\([x_a, x_b]\)) that bracket the root (i.e., \(sign~f(x_a) \ne sign~f(x_b)\)). Over at the FORTRAN 77 Netlib graveyard, the classic zeroin root finding method can be found that implements such an algorithm, namely Brent's method by Australian mathematician Richard Brent from 1971.

In my roots-fortran library, I have modern Fortran implementations of this method and many others, as well as a lot of benchmark test cases. The library can be compiled for single (real32), double (real64), or quad (real128) precision reals. For real64 numbers and convergence tolerance \(\epsilon\) = 1.0e-15, the following table lists the number of function evaluations for various methods to converge to the root for various test functions:

The methods are: BS = Bisection, PE = Pegasus, ITP, TM = TOMS 748, BR = Brent, AB = Anderson Bjorck, MU = Muller, and CH = Chandrupatla.

There seems to be cottage industry of creating slight variations of this type of method. There are any number of recent papers on this topic, all claiming some improvement. In reality, a lot of the newer methods really aren't great. An obvious trick is a choose functions where your method behaves well. Another is to only compare against very simple methods (such as bisection or regula falsi, which are known to be inefficient). Finally, the most devious trick is to list the number of iterations of the method, rather than the number of function evaluations, which is basically a useless comparison that obscures the real cost for methods that have multiple function evaluations per iteration.

The Brent method is still one of the best performers. The Chandrupatla method from 1997 is also quite good. For the full set of test cases, see the roots-fortran git repo.

Minimization methods

A related type of algorithm are derivative-free minimization methods. The FMIN function at Netlib is based on Brent's original localmin algorithm (modern version here). The PRIMA library implements modern versions of M.J.D. Powell's derivative-free minimizers for multidimensional functions.

See also

• R. P. Brent, "An algorithm with guaranteed convergence for finding a zero of a function", The Computer Journal, Vol 14, No. 4., 1971.
• R. P. Brent, "Algorithms for Minimization Without Derivatives", Prentice-Hall, 1973.
• Gregory W. Chicares, Brent's method outperforms TOMS 748 [Let me illustrate...]
• T. R. Chandrupatla, "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without derivatives", Advances in Engineering Software, Vol 28, 1997, pp. 145-149.
• M. J. D. Powell, "The BOBYQA algorithm for bound constrained optimization without derivatives," Department of Applied Mathematics and Theoretical Physics, Cambridge England, NA2009/06, 2009.
• polyroots-fortran -- Another modern Fortran library, specifically for finding the roots of polynomials.

Fortran 2023 Unleashed

Fortran 2023 (ISO/IEC 1539-1:2023) has been officially published. This replaces the previous standard (Fortran 2018). In a previous post I listed the updates that are in this release (then called "Fortran 202x"). It will be a while before we have fully-compliant Fortran 2023 compilers. In the latest update to the Intel Fortran Compiler, they have started to add a few minor things.

Fortran is the oldest surviving programming language still in common use. Originally developed in 1957, first standardized in 1966, and continuously updated since that time, Fortran remains at the heart of scientific and high-performance computing.

The next standard (unofficially called Fortran 202y) is also already in work.

See also

• Programming languages: Fortran, ISO, Nov. 17, 2023.
• S. Lionel, Fortran 2023 has been Published!, Nov. 23, 2023.
• John Reid, The new features of Fortran 2023, Mar. 13, 2023.
• J. Williams, The New Features of Fortran 202x, Mar 27, 2022.
• Fortran: High-performance parallel programming language

Fortran Filesystem Blues

One of my major pet peeves about Fortran is that it contains virtually no high-level access to the file system. The file system is one of those things that the Fortran standard pretends doesn't exist. It's one of those idiosyncratic things about the Fortran standard, like how it uses the word "processor" instead of "compiler", and doesn't acknowledge the existence of source code files so it doesn't bother to recommend a file extension for Fortran (gotta preserve that backward compatibility just in case punch cards come back).

There are various hacky workarounds to do various things that poor Fortran programmers have had to use for decades. Basically, all you have is OPEN, CLOSE, READ, WRITE, and INQUIRE. Here are a few examples:

Deleting a file

Here's a standard Fortran routine that can be used to delete a file:

function delete_file(name) result(istat)

  implicit none

  character(len=*),intent(in) :: name
  integer :: istat

  integer :: iunit

  open(name=name,newunit=iunit,status='OLD',iostat=istat)
  if (istat==0) close(iunit,status='DELETE',iostat=istat)

end function delete_file

As you can see, it's just a little trick, using a feature of CLOSE to delete a file after closing it (but first, we had to open it). Of course, note that the error codes are non-standard (the Fortran standard also doesn't consider this information important and lets the different compilers do whatever they want -- so of course, they all do something different). And of course, there is no exception handling in Fortran either (a post for another time).

Creating a directory

There is absolutely no standard ability in Fortran to create or delete a directory. Again, Fortran basically pretends that directories don't exist. Note that the Intel Fortran Compiler provides a super useful non-standard extension of the INQUIRE statement to allow it to be used to get information about directories. It's a pretty ridiculous state of affairs (Fortran added a bazillion IEEE routines in Fortran 2018 for some reason that nobody needed, but it still doesn't have something like this that everybody has needed for decades). Intel's IFPORT portability module provides many useful (non-standard) routines for accessing the file system (for example DELFILESQQ for deleting files and DELDIRQQ for deleting a directory). I use these all the time and the fact that they are non-standard and not present in other compilers is a major source of annoyance.

For some file or directory operations, you can always resort to using system calls, and thus have to provide different methods for different platforms (of course Fortran provides no standard way to get any information about the platform, so you have to resort to non-standard, possibly compiler-specific, preprocessor directives to do that). For example, to create a directory you could use:

subroutine make_directory(name)

  implicit none

  character(len=*),intent(in) :: name

  call execute_command_line ('mkdir '//name)

end subroutine make_directory

That one was easy because "mkdir" works on macOS, Window, and Linux. See previous post about getting command line calls into strings, which can also be useful for some things like this.

Getting a list of files

How about getting a list of all the files that match a specific pattern? Again, this can be done if you are using the Intel compiler using the IFPORT module GETFILEINFOQQ routine:

function get_filenames(pattern, maxlen) result(list)

  use ifport

  implicit none

  character(len=*),intent(in) :: pattern
  integer,intent(in) :: maxlen
  character(len=maxlen),dimension(:),allocatable :: list

  integer(kind=int_ptr_kind()) :: handle
  type(file$infoi8) :: info
  integer :: length

  allocate(list(0))
  handle = file$first
  do
    length = GETFILEINFOQQ(trim(pattern), info, handle)
    if ((handle==file$last) .or. (handle==file$error)) exit
    list = [list, info%name]
  end do

end function get_filenames

The file$infoi8 type is some super-weird non-standard Intel thing. Note that this example only returns the file names, not the full paths. To get the full path is left as an exercise to the reader. Also note that Fortran doesn't provide a good string class, so we just set a maximum string length as an argument to this function (but note that Intel is only returning a 255-character string anyway).

Final thoughts

There are various libraries out there that can do some of this. For example, M_system (a Fortran module interface for calling POSIX C system routines). Hopefully, a comprehensive set of filsystem routines will eventually make it into the Fortran Standard Library.

See also

• Portability Library Routines [Intel]
• getfileinfoqq -- Intel Portability Function that returns information about the specified file. [Intel]
• Interface to POSIX I/O API [Fortran stdlib issue]

Visual Studio Code + Fortran

With the right plugins, Microsoft Visual Studio Code can be turned into a very good modern Fortran IDE. Here are some of my favorites:

• FORTRAN IntelliSense (hansec.fortran-ls). To use this, you must pip install fortran-language-server.
• vscode-modern-fortran-formatter (yukiuuh.vscode-modern-fortran-formatter). To use this, you must pip install fprettify.
• Modern Fortran (krvajalm.linter-gfortran)
• Fortran Breakpoint Support (ekibun.fortranbreaker)

Historically, good full-featured IDEs for modern Fortran have been hard to find. Intel Fortran integrates with the full MS Visual Studio, and I used that for many years. But, except for syntax highlighting and automatic determination of compile order, it is basically terrible. I always have to disable the "IntelliSense" features since it brings the interface to a crawl (it's been like this for years and they don't seem to care enough to fix it). Unfortunately, Intel's excellent Fortran debugger also only works from within Visual Studio, so if you need to use it, you are stuck in that environment.

I've pretty much converted to use VSCode now for all Fortran coding (on both Mac and Windows). In the past, I have also used TextWrangler, Sublime Text, and Atom.

See also:

• vscode-icons -- A nice set of icons for various file formats, although the Fortran icon is terrible. I submitted a pull request to use the new one from fortran-lang.org. As of Oct 2021, it has yet to be merged.
• "Parse File…" with lagging visual studio [Intel Fortran Forum]. Note: I believe one of my comments on this post was censored by the admins. I distinctly remember asking if we were a bunch of suckers for continuing to give money to a multi-billion dollar corporation (Intel) when they can't be bothered to fix their defective product.
• Fortran IDEs at fortran-lang.org.
• Is Visual Studio Code Really The Best Code Editor?, February 28, 2021. [tabnine.com]

The Prehistory of Fortran-Lang

At the recent FortranCon2021, I was amused to see in the presentation "The State of Fortran 2021" that the first bullet point on the slide "Formation of Fortran-Lang" was:

August 2019 Conversations on Twitter between Ondřej Čertík, Milan Curcic, and Jacob Williams bring out common perceived shortcomings in the Fortran ecosystem

Laurence Kedward, et. al., "The State of Fortran 2021", Fortran2021, Sept. 2021

So, now, over two years later, I think it's finally time to tell my side of the story. 😀

I was first introduced to Fortran 25 years ago, in 1997 (my only previous programming experience was with BASIC on Apple ]['s). We used the WATFOR FORTRAN 77 compiler on DOS machines in an undergraduate programming class for engineers. Little did I know (or little did the university know) that Fortran 90 had been published years earlier (and indeed Fortran 95 was published at the end of that very semester).

Over the course of my academic career, I encountered and used Fortran a few times, mostly FORTRAN 77, and mostly to solve homework problems or do class projects. I came to use IMSL, Numerical Recipes in FORTRAN 77, Harwell Subroutine Library, and had some minimal exposure to Fortran 90 one semester when asked to modernize some old code by one of the professors for some project that I've forgotten about (something to do with gravity anomalies). But mostly I used a lot of Matlab, which is very popular among engineering students (due to cheap student licenses). I also used a little MathCad, some Mathematica, REALBasic (since renamed Xojo), and Visual Basic (classic). There was then no exposure in the aerospace engineering curriculum of two major US universities to any sort of formal software development practices (for all I know this is still the case). I never once encountered version control, unit testing, continuous integration, etc. Maybe those weren't as common then? (Subversion didn't come out until 2000, and Git until 2005).

When I entered professional life in late 2006, it wasn't long before I became the lead developer of a Fortran tool at NASA called Copernicus. Originally created by one of my former professors, Copernicus is a tool for spacecraft trajectory optimization and design. While working on Copernicus, I was introduced to Fortran 90/95 in earnest. Copernicus was originally developed on Windows with Compaq Visual Fortran (which, by the time I came along, had recently become defunct) and one of my first projects was to convert the code over to use the Intel Fortran compiler (which we have used ever since). During this period I first became aware of modern-ish (Fortran 95) concepts such as user derived types, operator overloading, function overloading, etc. Eventually Fortran 2003 came along, which was quite a revelation, and as soon as non-buggy implementations of the new features appeared in the Intel compiler (it took a while) I was using it like crazy. It's hard to imagine programming Fortran now without modern features such as type bound procedures. FORTRAN 77 (and indeed, even Fortran 95), to me, now seems like something out of the dark ages. And so I came to really appreciate the Modern Fortran programming language, and wondered why more people didn't use it (or seem to even know about it).

The fact was, all was not well in the larger Fortran community. Really, there was very little community. Most of the code I encountered on the internet was still FORTRAN 77, like the last 25 years never happened. There were also few if any Fortran blogs (there were a few of exceptions: Fortran Dev [now defunct], Fortran in a C World [now also defunct], and Steve Lionel's Doctor Fortran column at Intel [still updated infrequently]). One of the bright spots was the Intel Fortran Compiler forum, which was always full of helpful folks (Steve, of course, the most helpful of all). There were very few websites to speak of that had much useful Modern Fortran content. Fortran.com is a mostly-useless site, and fortran.org has been held hostage by a domain name squatter for years. Ondřej Čertík's site fortran90.org is good (but of course it is folly to put the standard name in the URL since that will eventually be out of date). Fortranwiki.org, administrated by Jason Blevins is also a good site.

On June 26, 2014, I started this blog, with the intent to focus on modern Fortran and astrodynamics algorithms. Earlier that year, I had also discovered GitHub and started to work on some projects there as well. At work, we were still using Subversion for version control at this time, but I wanted to learn more about git. The blog and my presence on GitHub was an experiment to attempt to publicly develop an ecosystem of Modern Fortran libraries (either by creating new ones, or refactoring old codes that hadn't been updated in decades that nobody else seemed to want to update), as well as to advocate for Fortran and show what it could do. On GitHub I encountered others who were also working on Fortran tools, namely Stefano Zaghi (author of lots of great libraries and tools, including FoBiS, the Fortran Building System for poor people, which is very dear to my heart), Zaak Beekman (who helped me greatly with JSON-Fortran and introduced me to continuous integration), Chris MacMackin (author of the ford automatic documentation generator for Fortran, which I contributed a very very small amount to), and a few others. Zaak also started a GitHub group called the Fortran F/OSS Programmers Group, which I also participated in and view as sort of a precursor to the fortran-lang group. Zaak and I coauthored a paper in 2018 on Modern Fortran applications in my field of spacecraft trajectory optimization, which I consider to be my magnum opus.

In April 2017, I was approached by Manning Publications, who wanted to discuss the possibility of writing a Fortran book. I didn't really have time to do that so I declined. But luckily, they also approached Milan Curcic, who did end up writing it (Modern Fortran, 2020) and it turned out excellent.

I currently have 67 repositories on GitHub, the vast majority of them being Fortran libraries, and they are used by lots of people all over the world in various fields of science and engineering. But one thing I didn't really have the aptitude or know-how to do was to form an online community. Years ago I considered buying the fortran.io domain when I noticed it was still up for grabs, but it was eventually bought by Nick Doiron around 2016 (see FORTRAN.io -- finally, a Fortran Web Framework) which made the rounds somewhat on the tech internet, basically just as a novelty. I guess the .io fad for tech domain names is over now anyway?

There just wasn't a website that really served as a focal point for the entire Fortran community. It was a huge image problem. This was one of the things I lamented in my (infamous?) August 2019 post. Luckily, Milan had purchased fortran-lang.org sometime in 2018, and that is what he activated for the new site in April 2020. Milan and Ondřej have already written about their amazing efforts to rejuvenate the Fortran community and ecosystem (I first became aware of Ondřej in April 2019 when he announced LFortran, which totally blew my mind). The fortran-lang website, the Fortran Standard Library, the Fortran Package Manager, and everything else that has come out these activities has been spectacular. I have actually contributed very little to all of this, although I did help with the logo (see this previous post for the story on that). I also continue to toil away on my libraries like before (recently getting a good number of them set up to compile with the new Fortran Package Manager).

So, now you know my side of the story. The future of Fortran is looking bright, and I'm pleased to have contributed in some small way. Now get off my lawn.

See also

• J. Williams, R. D. Falck, and I. B. Beekman, "Application of Modern Fortran to Spacecraft Trajectory Design and Optimization", 2018 Space Flight Mechanics Meeting, 8–12 January 2018, AIAA 2018-145.
• J. Williams, End of Life [degenerateconic.com] August 3, 2019
• J. Williams, New Blood [degenerateconic.com] April 26, 2020
• M. Curcic, First year of Fortran-lang [medium.com] Dec 21, 2020
• O. Čertík, Resurrecting Fortran [ondrejcertik.com] March 12, 2021

JSON

There are a lot a JSON libraries available for Fortran nowadays. The complete list, as far as I know is:

• fson Fortran 95 JSON Parser. JSON-Fortran is a fork of this library. (started 2012).
• YAJL-Fort A Modern Fortran Interface to YAJL. This one is an interface to a C library, and not pure Fortran. See also petaca. (started 2013).
• JSON-Fortran A Fortran 2008 JSON API. This is my library. As far as I know, this was the first production-ready JSON parser written in modern Fortran. (started 2014).
• fortjson JSON library written in Fortran 2003. Designed with portability across HPC architectures in mind. (started 2018).
• jsonff JSON for Fortran. (started 2019).

Example

Here's an example of using JSON-Fortran to read in some data from a JSON string:

program test

 use json_module, wp => json_RK, ip => json_IK

 implicit none

 type(json_file) :: json
 integer(ip) :: t
 real(wp),dimension(:),allocatable :: x

 call json%deserialize('{"t": 1, "x": [2.0, 3.0, 4.0]}')

 call json%get('t', t); write(*,*) t
 call json%get('x', x); write(*,*) x

end program test

This prints:

           1
   2.00000000000000        3.00000000000000        4.00000000000000

Note that JSON-Fortran has all sorts of options for:

• Using different real (real32, real64, and real128) and integer kinds (int8, int16, int32, and int64).
• Controlling the JSON format for printing, including number of spaces for indenting, printing vectors on one line only, or full minification with no extra whitespace or line breaks.
• Support for comments in a JSON file.
• Multiple ways to get data from a JSON structure, such as RFC 6901 "JSON Pointer" paths or JSONPath) "bracket-notation".
• Graceful handing of unusual inputs such as NaN or Infinity.
• Thread safe error handling.

See also

• json.org
• Robinson, T. E., Replacing Fortran Namelists with JSON, American Geophysical Union, Fall Meeting 2017.
• JSON + Fortran [degenerateconic.com], June 28, 2014.
• Some projects using JSON-Fortran [JSON-Fortran wiki]

oneAPI

Intel recently released its new oneAPI Toolkit, which had been in beta for a while. According to their documentation, oneAPI includes various products:

• Intel® oneAPI Base Toolkit
  • Includes the DPC++/C++ Compiler, the Intel Math Kernel Library, and the Intel Distribution for Python
• Intel® oneAPI HPC Toolkit
  • Includes the new version of the Intel Fortran Compiler (now known as the Classic version) as well as a new LLVM-based Fortran compiler that is currently in beta. Also includes Intel MPI and Intel Inspector.
• Intel® oneAPI IoT Toolkit
• Intel® oneAPI Rendering Toolkit
• Intel® AI Analytics Toolkit
• Intel® Distribution of OpenVINO™ Toolkit

The big news here is that this is all now entirely free to use. They are all available for Windows, Linux, and macOS. The Fortran compiler also now includes all of the Fortran 2018 standard. Apparently, Intel will still sell you a license, but it seems the only reason to buy one would be to get the Premier Support.

See also

• Intel oneAPI Toolkits [Intel]
• Intel's Fortran 2018 now released [comp.lang.fortran]
• Intel releases oneAPI Toolkit, free, Fortran 2018 [fortran-lang.discourse.group]