Degenerate Conic

Refactoring (Part 2)

curve

I'm working on a modern Fortran version of the simulated annealing optimization algorithm, using this code as a starting point. There is a Fortran 90 version here, but this seems to be mostly just a straightforward conversion of the original code to free-form source formatting. I have in mind a more thorough modernization and the addition of some new features that I need.

Refactoring old Fortran code is always fun. Consider the following snippet:

C  Check termination criteria.
      QUIT = .FALSE.
      FSTAR(1) = F
      IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE.
      DO 410, I = 1, NEPS
        IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE.
410   CONTINUE

This can be greatly simplified to:

! check termination criteria.
fstar(1) = f
quit = ((fopt-f)<=eps) .and. (all(abs(f-fstar)<=eps)

Note that we are using some of the vector features of modern Fortran to remove the loop. Consider also this function in the original code:

      FUNCTION  EXPREP(RDUM)
C  This function replaces exp to avoid under- and overflows and is
C  designed for IBM 370 type machines. It may be necessary to modify
C  it for other machines. Note that the maximum and minimum values of
C  EXPREP are such that they has no effect on the algorithm.

      DOUBLE PRECISION  RDUM, EXPREP

      IF (RDUM .GT. 174.) THEN
         EXPREP = 3.69D+75
      ELSE IF (RDUM .LT. -180.) THEN
         EXPREP = 0.0
      ELSE
         EXPREP = EXP(RDUM)
      END IF

      RETURN
      END

It is somewhat disturbing that the comments mention IBM 370 machines. This routine is unchanged in the f90 version. However, these numeric values are no longer correct for modern hardware. Using modern Fortran, we can write a totally portable version of this routine like so:

pure function exprep(x) result(f)

use, intrinsic :: ieee_exceptions

implicit none

real(dp), intent(in) :: x
real(dp) :: f

logical,dimension(2) :: flags
type(ieee_flag_type),parameter,dimension(2) :: out_of_range = [ieee_overflow,ieee_underflow]

call ieee_set_halting_mode(out_of_range,.false.)

f = exp(x)

call ieee_get_flag(out_of_range,flags)
if (any(flags)) then
  call ieee_set_flag(out_of_range,.false.)
  if (flags(1)) then
    f = huge(1.0_dp)
  else
    f = 0.0_dp
  end if
end if

end function exprep

This new version is entirely portable and works for any real kind. It contains no magic numbers and uses the intrinsic IEEE exceptions module of modern Fortran to protect against overflow and underflow, and the huge() intrinsic to return the largest real(dp) number.

So, stay tuned...

End of Life

python2clock

It amuses me somewhat to see the push to get people to stop using Python 2. Python 2 was sort of replaced with Python 3 in 2008. However, Python 3 broke backward compatibility, and Python 2 continued to be supported and even developed (some of the Python 3 features have even been backported to Python 2 which is somewhat bizarre). However, the official Python maintainers have declared January 1st, 2020 to be the end of life for Python 2. However, many Linux distos (as well as MacOS) still include Python 2 as the default.

The situation with Fortran is significantly different. You can still get (even purchase) a Fortran compiler that can compile obsolescent FORTRAN 77 written 40 years ago. A lot of people haven't gotten the memo that Fortran has moved on. Not a day goes by where there isn't a StackOverflow question from somebody about some godawful FORTRAN 77 code that they can't figure out how to get working. So, I wonder how successful this effort will be with Python. Will there still be Python 2 holdouts 30 years from now? Of course, the Python community is significantly different from the Fortran community (to whatever extent Fortran could be said to even have a community). The Python implementation is open source so people are free in theory to just fork it and continue to update it. But, as the official source dries up I suspect eventually people will just move on.

There is basically no one centralized place for Fortran users to learn about Fortran, download Fortran, work on or even comment on changes to Fortran, or anything else really. Backward compatibility is actually one of the major strengths of Fortran, but there just isn't anyone to tell you to move on. Many of the major Fortran libraries freely available on the internet are still, in 2018, written in Fortran 77 (see NetLib, where Fortran codes go to die). There is also SOFA and SPICELIB (two libraries still being developed in Fortran 77 for some reason). There is LAPACK, probably the most visible FORTRAN 77 library. It turns out, you can't even link LAPACK and SPICELIB in the same program anymore, because both now have a routine called DPSTRF! If these libraries weren't using this depreciated source form and were actually using Fortran 90 modules (or god forbid, Fortran 2008 submodules) we wouldn't have this problem. Nobody should be publishing a Fortran library at this point that is just a pile of subroutines. And yet, they do.

Sure, there is a Fortran standards committee. But they write the standard and don't do much else. There isn't any official website that serves as any kind of central location to learn about the language or download any compilers. There is no nonprofit "Fortran Foundation". There isn't even a Fortran logo. What we basically have is a few compiler vendors such as NAG, Intel, PGI, as well as the open source gfortran developers, each with their own websites and their own schedules. Even the C++ language (another ISO standard language) has a GitHub site and even an actual homepage. The Fortran committee website is extremely underwhelming, and is basically just a list of links to plain text files that have meeting minutes and some other documents. The odds of a new user even finding this site are pretty slim (and of course, there really isn't anything useful here anyway). The Fortran standards process seems very opaque (I guess it involves infrequent meetings and typing up plain text files).

In many ways, the standards committee has failed the user base. The slowness of the language development process, the refusal to adopt modern practices of collaboration, and the refusal to address major shortcomings of the language has allowed other languages to overtake Fortran in the fields it was once dominant. In fields such as machine learning, one of the major computational activities of modern times, there is no significant presence of Fortran. The committee has given us incremental, somewhat half-baked features that don't really solve real-world problems in any kind of satisfactory manner (like ISO_VARYING_STRING, parametrized derived types, user-defined IO, floating point exceptions). While shortcomings from decades ago are left unaddressed (lack of a useful mechanism for generic programming, lack of access to the file system, no intrinsic high-level classes for dynamic structures such as stacks and queues, no useful varying string class, a non-standard module file format that is a constant source of annoyance for both users and library developers, even the lack of a standardized file extension for Fortran, which leads to nonsense like .f90, .f03, .f08, .f18? being used). It's incredible but true that the Fortran standard doesn't actually define a file extension to be used for Fortran. One of the respondents to a recent survey on Fortran language development had this to say:

The standard committee is too inbred with compiler developers who only see the language from the inside out and lacking in users who know what features they need for their particular application space.

Amazing libraries in the scientific/technical/numerical domain are being written in other programming languages. DifferentialEquations.jl is the type of high-quality math-oriented library that no one is writing anymore in Fortran (indeed it makes use of features of Julia that aren't even possible in Fortran). This article about using Julia as a differentiable language is the sort of thinking we desperately need in the Fortran world. Not more excuses about why fundamental changes can't be made to Fortran for fear of breaking somebody's 60 year old spaghetti code.

GFortran 9.1

gccegg-65

Gfortran 9.1 (part of GCC) has been released. Apparently this is a significant GCC release with a "huge number of improvements" including a new D language component. Of course, the Fortran updates are a little more modest. According to the release notes, the updates are:

Asynchronous I/O is now fully supported [Fortran 2003].
The BACK argument for MINLOC and MAXLOC has been implemented [Fortran 2008].
The FINDLOC intrinsic function has been implemented [Fortran 2008].
The IS_CONTIGUOUS intrinsic function has been implemented [Fortran 2008].
Direct access to the real and imaginary parts of a complex variable via c%re and c%im has been implemented [Fortran 2008].
Type parameter inquiry via str%len and a%kind has been implemented [Fortran 2008].
C descriptors and the ISO_Fortran_binding.h source file have been implemented [Fortran 2018].
The MAX and MIN intrinsics are no longer guaranteed to return any particular value in case one of the arguments is NaN.
A new command-line option -fdec-include, has been added as an extension for compatibility with legacy code using some non-standard behavior from the old DEC compiler.
A new BUILTIN directive, has been added. The purpose of the directive is to provide an API between the GCC compiler and the GNU C Library which would define vector implementations of math routines.

In addition, the release includes a bunch of bug fixes. Gfortran has more-or-less complete support for Fortran 2003, and only a couple things missing from Fortran 2008. There is a ways to go for full Fortran 2018 support. Gfortran is maintained by a very small number of volunteers, and their hard work is greatly appreciated!

Fortran + LLVM Update

There are now at least five open source Fortran compilers (in various stages of completion) that are based on LLVM:

Flang -- Original attempt (possibly defunct?) by the LLVM Team at the University of Illinois at Urbana-Champaign.
Flang -- The first attempt by NVIDIA/PGI. It's some kind of open-sourced version of their commercial compiler being funded by Lawrence Livermore, Sandia and Los Alamos National Laboratories. See previous post when this was first announced in 2015.
f18 -- The second (modernized) attempt by NVIDIA/PGI, built from the ground up. Intended to be a replacement for the previous one.
DragonEgg -- This one uses LLVM as a GCC backend. Also seems to be defunct, the website says it only works with the very old GCC 4.6.
lfortran -- This one is very interesting, since it isn't just a run of the mill Fortran compiler, it extends the language a little bit to include a REPL and some other great ideas. This one seems to have some connection to Los Alamos National Laboratory as well, but doesn't appear to be related to the NVIDIA Flang one.

None of these really seem to be finished yet. Hopefully, one or more will achieve full Fortran 2018 compliance and be good enough for production work. I'm particular interested to see how lfortran matures. In recent years, LLVM has taken the compiler world by storm, and it will be nice to see Fortran get in on the action.

References

The LLVM Compiler Infrastructure Project
M. Larabel, NVIDIA Has Been Working On A New Fortran "f18" Compiler It Wants To Contribute To LLVM, 1 March 2019 [Phoronix]
NNSA, national labs team with Nvidia to develop open-source Fortran compiler technology, Nov. 13, 2015 [Lawrence Livermore National Laboratory]

posted at 21:29 · Programming · Flang Fortran LFortran LLVM

Old School

I came across this old NASA document from 1963, where a cross product subroutine is defined in Fortran (based on some of the other routines given, it looks like they were using FORTRAN II):

cross

A cross product function is always a necessary procedure in any orbital mechanics simulation (for example when computing the angular momentum vector \(\mathbf{h} = \mathbf{r} \times \mathbf{v}\)). The interesting thing is that this subroutine will still compile today with any Fortran compiler. Of course, there are a few features used here that are deemed obsolescent in modern Fortran (fixed-form source, implicit typing, and the DIMENSION statement). Also it is using single precision reals (which no one uses anymore in this field). A modernized version would look something like this:

subroutine cross(a,b,c)
implicit none
real(wp),dimension(3),intent(in) :: a,b
real(wp),dimension(3),intent(out) :: c
c(1) = a(2)*b(3)-a(3)*b(2)
c(2) = a(3)*b(1)-a(1)*b(3)
c(3) = a(1)*b(2)-a(2)*b(1)
end subroutine cross

But, basically, except for declaring the real kind, it does exactly the same thing as the one from 1963. The modern routine would presumably be put in a module (in which the WP kind parameter is accessible), for example, a vector utilities module. While a subroutine is perfectly acceptable for this case, my preference would be to use a function like so:

pure function cross(a,b) result(c)
implicit none
real(wp),dimension(3),intent(in) :: a,b
real(wp),dimension(3) :: c
c(1) = a(2)*b(3)-a(3)*b(2)
c(2) = a(3)*b(1)-a(1)*b(3)
c(3) = a(1)*b(2)-a(2)*b(1)
end function cross

This one is a function that return a \(\mathrm{3} \times \mathrm{1}\) vector, and is explicitly declared to be PURE (which can allow for more aggressive code optimizations from the compiler).

This document, which is the programmer's manual for an interplanetary error propagation program, contains a lot of other Fortran gems. It also has some awesome old school flow charts:

Fortran 2018 Unleashed

fortran-2018

Fortran 2018 (ISO/IEC 1539:2018) has been officially published by the ISO, so it is now the official Fortran standard, replacing the nearly decade-old Fortran 2008 (which was actually published in 2010). Since this is a copyrighted ISO document, you can't actually view it online for free, but you can view the final draft which is basically the same thing.

The final draft document for Fortran 2018 is 630 pages. By contrast, the original Fortran Programmer's Reference Manual, produced by IBM in 1956, was only 51 pages. The language has certainly changed a lot in 62 years. Modern Fortran is still the programming language of choice for many of us solving very large and complex computational problems in science and engineering. It is the only scientific programming language that does not compromise on speed in order to provide dynamic typing, a REPL, JIT, or whatever the latest computing fad is. It is the only scientific programming language that is an international standard. It is the only scientific programming language that has multiple complete implementations from different vendors. It has a far superior array and matrix syntax than C-based languages. Its syntax is straightforward and can be mastered by non-professional programmers. It's hard to write slow code in Fortran, even if you barely know what you're doing. Modern Fortran is object-oriented and also provides a standardized interoperability with C-based languages (nowadays this is quite useful for interfacing with Python). Fortran-based tools are used by NASA to design and optimize the trajectories of interplanetary spacecraft. The next generation of crewed spacecraft missions are being designed with Fortran-based tools.

Now, Fortran is not without its flaws (perhaps I will write a post one day on the things I don't like about Fortran), but it is still very good at doing what it was designed to do. Improvements added in each new revision continue to breathe new life into this venerable programming language.

Modern Fortran Explained

The new edition of "Modern Fortran Explained" is out. This is the latest in the venerable series that provides an informal but mostly-definitive description of the Fortran standard (of course, the actual ISO standard is not that useful to anybody trying to learn the language, since it's very technical and almost entirely unreadable). The latest edition of Modern Fortran Explained incorporates the updates that will be included in the upcoming Fortran 2018 standard. I highly recommend this book. However, the authors chose to present the new stuff at the end of the book, rather than in the main text. This is somewhat unfortunate. For example, the do concurrent construct description begins on p 137. However, additional features added in Fortran 2008 and 2018 are not described or even mentioned in this section, for that you have to flip to p 435. So, all the information is there, but the organization could be improved to be of more use for beginners.

Intel Fortran Compiler 19.0

Intel has just released v19 of the Intel Fortran Compiler (part of Intel Parallel Studio XE 2019). The new release adds some new features from the upcoming Fortran 2018 standard:

Coarray events
Intrinsic function coshape
Default accessibility for entities accessed from a module
Import Enhancements
All standard procedures in ISO_C_BINDING other than C_F_POINTER are now PURE

Presumably, it also include some bug fixes. Intel used to provide a list of bug fixes, but they seem to have stopping doing that for some reason.

Circular Restricted Three-Body Problem

*Libration Points in the Earth-Moon System. From [Farquhar, 1971](https://www.lpi.usra.edu/lunar/documents/NASA%20TN%20D-6365.pdf).*

The Circular Restricted Three-Body Problem (CR3BP) is a continuing source of delight for astrodynamicists. If you think that, after 200 years, everything that could ever be written about this subject has already been written, you would be wrong. The basic premise is to describe the motion of a spacecraft in the vicinity of two celestial bodies that are in circular orbits about their common barycenter. The "restricted" part means that the spacecraft does not affect the motion of the other two bodies (i.e., its mass is negligible).

The normalized equations of motion for the CR3BP are very simple (see the Fortran Astrodynamics Toolkit):

subroutine crtbp_derivs(mu,x,dx)

implicit none

real(wp),intent(in) :: mu
!! CRTBP parameter
real(wp),dimension(6),intent(in) :: x
!! normalized state [r,v]
real(wp),dimension(6),intent(out) :: dx
!! normalized state derivative [rdot,vdot]

real(wp),dimension(3) :: r1,r2,rb1,rb2,r,v,g
real(wp) :: r13,r23,omm,c1,c2

r = x(1:3) ! position
v = x(4:6) ! velocity
omm = one - mu
rb1 = [-mu,zero,zero] ! location of body 1
rb2 = [omm,zero,zero] ! location of body 2
r1 = r - rb1 ! body1 -> sc vector
r2 = r - rb2 ! body2 -> sc vector
r13 = norm2(r1)**3
r23 = norm2(r2)**3
c1 = omm/r13
c2 = mu/r23

!normalized gravity from both bodies:
g(1) = -c1*(r(1) + mu) - c2*(r(1)-one+mu)
g(2) = -c1*r(2) - c2*r(2)
g(3) = -c1*r(3) - c2*r(3)

! derivative of x:
dx(1:3) = v ! rdot
dx(4) = two*v(2) + r(1) + g(1) ! vdot
dx(5) = -two*v(1) + r(2) + g(2) !
dx(6) = g(3) !

end subroutine crtbp_derivs

There are many types of interesting periodic orbits that exist in this system. I've mentioned Distant Retrograde Orbits (DROs) here before. Probably the most well-known CR3BP periodic orbits are called "halo orbits". Halo orbits can be computed in a similar manner as DROs (see previous post). For halos, we can get a good initial guess using an analytic approximation, and then use this guess to target the real thing. See the halo orbit module in the Fortran Astrodynamics Toolkit for an example of this approach. We can also use a continuation type method to generate families of orbits (say, for the L2 libration point, we have a range of halo orbits from near-planars that "orbit" L2, to near-rectilinears with close approaches to the secondary body). In a realistic force model (with the real ephemeris of the Earth, Moon, and Sun), the orbits generated using CR3BP assumptions are not periodic, but we can use the CR3BP solution as an initial guess in a multiple-shooting type approach with a numerical nonlinear equation solver to generate quasi-periodic multi-rev solutions. See the references below for details.

I generated the following image of a bunch of planar periodic CR3BP orbits using a program incorporating the Fortran Astrodynamics Toolkit and OpenFrames (an interface to the open source OpenSceneGraph library that can be used to visualize spacecraft trajectories). The initial states for the various orbits were taken from this database. Eventually, I will make this program open source when it is finished. There's also an online CR3BP propagator here.

screenshot5

CMLIB

nbs

The National Bureau of Standards (NBS) Core Math Library (CMLIB) is a "collection of high-quality, easily transportable Fortran subroutine sublibraries solving standard problems in many areas of mathematics and statistics". It was written in FORTRAN 77 and is available from a NIST FTP link. The first version (1.0) was released in March 1986, and the last update (3.0) was in 1988. The software is no longer maintained, and I assume one day the link will probably disappear. It was compiled mostly from other externally-available libraries available at the time (including BLAS, EISPACK, FISHPAK, FNLIB, FFTPACK, LINPACK, and QUADPACK) but there are also some original codes. My Bspline-Fortran package is a modernized update and extension of the DTENSBS routines from CMLIB. The core spline interpolation routines (slightly updated) still work great after more than 30 years.

The National Bureau of Standards (NBS) was renamed the National Institute of Standards and Technology (NIST) in August 1988 (a few months after CMLIB was last updated). The nice thing about these old government-produced codes is that they are public domain in the United States. Indeed any software written by a U.S. government employee (but not a government contractor) as part of their official duties is automatically public domain. I generally prefer permissive software licenses, and you can't get more permissive than that.

Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Sep 15, 2019

Refactoring (Part 2)

See also

Aug 03, 2019

End of Life

See also

May 04, 2019

GFortran 9.1

See also

Mar 23, 2019

Fortran + LLVM Update

References

Dec 10, 2018

Old School

See also

Dec 03, 2018

Fortran 2018 Unleashed

See also

Nov 10, 2018

Modern Fortran Explained

See also

Sep 13, 2018

Intel Fortran Compiler 19.0

See also

Aug 05, 2018

Circular Restricted Three-Body Problem

See also

Jul 14, 2018

CMLIB

See also