Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Feb 10, 2018

Another One Bites the Dust: GRAM Atmosphere Model

It looks like the Fortran community has lost another venerable NASA Fortran library to C++. I recently noticed that the latest release of NASA's Global Reference Atmosphere Model (GRAM) is now C++. From the release link:

Earth-GRAM 2016 is now available as an open-source C++ computer code that can run on a variety of platforms including PCs and UNIX stations. The software provides a model that offers values for atmospheric parameters such as density, temperature, winds, and constituents for any month and at any altitude and location within the Earth's atmosphere. Earth-GRAM 2010 is available in FORTRAN. Similar computer models of Mars, Venus, Titan, and Neptune are available for diverse mission applications through NASA’s Space Environments and Effects (SEE) Program, a partnership with industry, academia, and other government agencies that seeks to develop more reliable, more effective spacecraft.

globe_west

The Earth-GRAM model, originally written in Fortran 77, is produced by MSFC and has been around for a while, with many revisions (i.e., GRAM-86, GRAM-88, GRAM-90, GRAM-95, GRAM-99, GRAM 2007, and GRAM 2010). In the 2010 revision, they converted the old Fortran 77 code to Fortran 90, which I thought was a step in the right direction (although maybe 20 years too late). It seems like they could have used Fortran 2003 and provided a C interface using the standard Fortran/C Interoperability feature, which would have enabled calling it from a variety of languages (including C++).

I remain utterly unconvinced that C++ is the programming language that engineers and scientists should be using. Are CS people taking over programming jobs once done by engineers? Is this a failure of university engineering departments? Although it seems like even CS people don't want to be using C++ anymore. Every time I turn around a new programming language (Java, Rust, Go, D, C#, Swift, ...) is created by CS people fed up with the C++ dumpster fire. So I don't know who is telling engineers that C++ is the future. And the perception that Fortran is an obsolete programming language is very strong, even within the engineering/technical disciplines that it was designed for. Sometimes this is due to sheer ignorance, but sometimes it is due to ancient Fortran 77 legacy code that people have to deal with. I sympathize with them, Fortran 77 code is terrible and needs to go away, but it should be modernized, not thrown out and replaced with C++. Think of the children! (and the memory leaks!)

Parting words of wisdom

Security tips when programming in C (2017 edition):

1) Stop typing

2) Delete what you've already typed

— ryan huber (@ryanhuber) June 21, 2017

See also

Dec 10, 2017

Fortran 2018

fortran2018

The upcoming Fortran standard formerly known as Fortran 2015 has a new name: Fortran 2018. It was decided to change it in order to match the expected year of publication. This makes sense. The previous standard (Fortran 2008) was published in 2010.

Waiting for an updated Fortran standard is an exercise in Zen-like patience. Almost a decade after Fortran 2008, we'll get a fairly minor update to the core language. And it will be years after that before it's fully supported by any compiler that most users will have available (gfortran still doesn't have a bug-free implementation of all of Fortran 2008 or even 2003). Fortran was essentially reborn in the Fortran 2003 standard, which was an amazing update that brought Fortran into the modern world. It's a terrific programming language for scientific and technical computing. However, the limitations are all too clear and have been for a long time:

  • We need better facilities for generic programming. It's impossible to do some things without having to duplicate code or use "tricks" like preprocessing or include files (see JSON-Fortran for examples).
  • We need some kind of exception handling. Fortran does have a floating-point exception handling feature, but honestly, it's somewhat half-baked.
  • We need a better implementation of strings. Allocatable strings (introduced in Fortran 2003) are great, but not enough, since they can't be used in all instances where strings are needed.
  • We need the language to be generally less verbose. I'm tired to having to type multiple nested SELECT TYPE statements to do something that is a one liner in Python (sure I know Fortran will never be as succinct as Python, but some of the verbosity required for object-oriented Fortran is just perverse).
  • We need any number of new features to make it easier to extend the language with third-party libraries (so we don't have to wait two decades for a feature we want).
  • We also need the language development process to embrace a more open collaborative model using modern tools (Usenet is not the future). I guess the recent survey was unprecedented, but it's not enough.

fortran

Fortran is a programming language that needs a better PR department. Legacy Fortran codes are being rewritten in C++, Python, or even Julia. NASA frequently throws massive Fortran 77 libraries with decades of heritage (e.g., DPTRAJ/ODP, GTDS, SPICELIB) into the trash in order to rewrite it all from the ground up in C++, without ever considering Fortran 2003+ (or maybe not realizing it exists?). The information about modern Fortran on the internet is spotty at best, outdated, or downright wrong (what is the deal with REAL*8?). In popular consciousness Fortran is mostly a punchline (usually something to do with punchcards and your granddad). A language like Python (which was never designed for technical computing) is now seen by many as a superior solution for technical computing. Matlab refers to itself as "the only top programming language dedicated to mathematical and technical computing"! The Julia website lists somewhat misleading benchmarks than implies that C, Julia, and even Lua are faster than Fortran.

Now, get off my lawn, you kids!

db3996347690687770d0e98039208811

See also

Aug 11, 2017

Another One Bites the Dust

jpl

JPL recently released an update to their awesome SPICE Toolkit (it is now at version N66). The major new feature in this release is the Digital Shape Kernel (DSK) capability to define the shapes of bodies (such as asteroids) via tessellated plate models.

Unfortunately for Fortran users, they also announced that they have decided to reimplement the entire library in C++. SPICELIB is currently written in Fortran 77, which they f2c to provide a C version (which is also callable from IDL, Matlab, and Python, among others). Their reason for this "upgrade" is to provide thread safety and object oriented features. Of course, modern Fortran can be thread safe and object oriented, and upgrading the code to modern standards could be done in a fraction of the time it will take to rewrite everything from scratch in C++. SPICELIB is extremely well-written Fortran 77 code, and is not infested with COMMON blocks, EQUIVALENCE statements, etc. I actually don't think it would take much effort to modernize it. In addition, Fortran/C interoperability could be employed to easily provide an interface that is callable from C without source transformation.

However, I guess it isn't meant to be, and the science/engineering community will lose another Fortran code to C++ like many times before, in spite of C++ being a terrible language for scientists and engineers.

See also

Oct 29, 2016

GOTO Still Considered Harmful

For a number of years I have been familiar with the observation that the quality of programmers is a decreasing function of the density of go to statements in the programs they produce. -- Edsger W. Dijkstra

goto

One of the classics of computer science is Edsger Dijkstra's "Go To Statement Considered Harmful", written in 1968. This missive argued that the GOTO statement (present in several languages at the time, including Fortran) was too primitive for high-level programming languages, and should be avoided.

Most people now agree with this, although some even today think that GOTOs are fine under some circumstances. They are present in C and C++, and are apparently used extensively in the Linux kernel. A recent study of C code in GitHub concluded that GOTO use is not that bad and is mainly used for error-handling and cleanup tasks, for situations where there are no better alternatives in C. However, C being a low-level language, we should not expect much from it. For modern Fortran users however, there are better ways to do these things. For example, unlike in C, breaking out of multiple loops is possible without GOTOs by using named DO loops with an EXIT statement like so:

a_loop : do a=1,n
    b_loop: do b=1,m
        !do some stuff ...
        if (done) exit a_loop ! break out of the outer loop
        !...
    end do b_loop
end do a_loop

In old-school Fortran (or C) this would be something like this:

do a=1,n
    do b=1,m
        ! do some stuff ...
        if (done) goto 10 ! break out of the outer loop
        ! ...
    end do
end do
10 continue

Of course, these two simple examples are both functionally equivalent, but the first one uses a much more structured approach. It's also a lot easier to follow what is going on. Once a line number is declared, there's nothing to stop a GOTO statement from anywhere in the code from jumping there (see spaghetti code). In my view, it's best to avoid this possibility entirely. In modern Fortran, DO loops (with CYCLE and EXIT), SELECT CASE statements, and other language constructs have obviated the need for GOTO for quite some time. Fortran 2008 added the BLOCK construct which was probably the final nail in the GOTO coffin, since it allows for the most common use cases (exception handing and cleanup) to be easily done without GOTOs. For example, in this code snippet, the main algorithm is contained within a BLOCK, and the exception handling code is outside:

main: block
    ! do some stuff ...
    if (error) exit main ! if there is a problem anywhere in this block,
    ! then exit to the exception handling code.
    ! ...
    return ! if everything is OK, then return
end block main

! exception handling code here

The cleanup case is similar (which is code that is always called):

main: block
    ! do some stuff ...
    if (need_to_cleanup) exit main ! for cleanup
    ! ...
end block main

! cleanup code here

I don't believe any of the new programming languages that have cropped up in the past couple of decades has included a GOTO statement (although someone did create a goto statement for Python as an April Fool's joke in 2004). Of course, the presence of GOTO's doesn't mean the programmer is bad or that the code isn't going to work well. There is a ton of legacy Fortran 77 code out there that is rock solid, but unfortunately littered with GOTOs. An example is the DIVA integrator from the JPL MATH77 library (the screenshot above is from this code). First written in 1987, it is code of the highest quality, and has been used for decades in many spacecraft applications. However, it is also spaghetti code of the highest order, and seems like it would be unimaginably hard to maintain or modify at this point.

Source: XKCD

Source: XKCD

See also

Apr 24, 2016

Fortran & C Interoperability (Part 2)

C

Here's another example using the C interoperability features of modern Fortran. First introduced in Fortran 2003, this allows for easily calling C routines from Fortran (and vice versa) in a standard and portable way. Further interoperability features will also be added in the next edition of the standard.

For a Fortran user, strings in C are pretty awful (like most things in C). This example shows how to call a C function that returns a string (in this case, the dlerror function).

module get_error_module

use, intrinsic :: iso_c_binding

implicit none

private

public :: get_error

interface
    !interfaces to C functions
    function dlerror() result(error) &
    bind(C, name='dlerror')
    import
    type(c_ptr) :: error
    end function dlerror

    function strlen(str) result(isize) &
    bind(C, name='strlen')
    import
    type(c_ptr),value :: str
    integer(c_int) :: isize
    end function strlen
end interface

contains

function get_error() result(error_message)
!! wrapper to C function char *dlerror(void);

character(len=:),allocatable :: error_message

type(c_ptr) :: cstr
integer(c_int) :: n

cstr = dlerror() ! pointer to C string

if (c_associated(cstr)) then

    n = strlen(cstr) ! get string length

    block
        !convert the C string to a Fortran string
        character(kind=c_char,len=n+1),pointer :: s
        call c_f_pointer(cptr=cstr,fptr=s)
        error_message = s(1:n)
        nullify(s)
    end block

else
    error_message = ''
end if

end function get_error

end module get_error_module

First we define the bindings to two C routines so that we can call them from Fortran. This is done using the INTERFACE block. The main one is the dlerror function itself, and we will also use the C strlen function for getting the length of a C string. The bind(C) attribute indicates that they are C functions. The get_error function first calls dlerror, which returns a variable of type(c_ptr), which is a C pointer. We use c_f_pointer to cast the C pointer into a Fortran character string (a CHARACTER pointer variable with the same length). Note that, after we know the string length, the block construct allows us to declare a new variable s of the length we need (this is a Fortran 2008 feature). Then we can use it like any other Fortran string (in this case, we assign it to error_message, the deferred-length string returned by the function).

Of course, Fortran strings are way easier to deal with, especially deferred-length (allocatable) strings, and don't require screwing around with pointers or '\0' characters. A few examples are given below:

subroutine string_examples()

implicit none

!declare some strings:
character(len=:),allocatable :: s1,s2

!string assignment:
s1 = 'hello world' !set the string value
s2 = s1 !set one string equal to another

!string slice:
s1 = s1(1:5) !now, s1 is 'hello'

!string concatenation:
s2 = s1//' '//s1 ! now, s2 is 'hello hello'

!string length:
write(*,*) len(s2) ! print length of s2 (which is now 11)

!and there are no memory leaks,
!since the allocatables are automatically
!deallocated when they go out of scope.

end subroutine string_examples

See also

Apr 09, 2016

C++ vs Fortran (Part 3)

matrix

Built-in high-level support for arrays (vectors and matrices) using a very clean syntax is one of the areas where Fortran really shines as a programming language for engineering and scientific simulations. As an example, consider matrix multiplication. Say we have an M x K matrix A, a K x N matrix B, and we want to multiply them to get the M x N matrix AB. Here's the code to do this in Fortran:

AB = matmul(A,B)

Now, here's how to do the same thing in C or C++ (without resorting to library calls) [1]:

for (int i = 0; i < M; i++) {
    for (int j = 0; j < N; j++) {
        AB[i][j]=0;
        for (int k = 0; k < K; k++) {
            AB[i][j] += A[i][k] * B[k][j];
        }
    }
}

MATMUL is, of course, Fortran's built-in matrix multiplication function. It was added to the language as part of the Fortran 90/95 standards which dramatically upgraded the language's array-handling facilities. Later standards also added additional array features, and Fortran currently contains a wide range of intrinsic array functions:

Routine Description
ALL True if all values are true
ALLOCATED Array allocation status
ANY True if any value is true
COUNT Number of true elements in an array
CSHIFT Circular shift
DOT_PRODUCT Dot product of two rank-one arrays
EOSHIFT End-off shift
FINDLOC Location of a specified value
IPARITY Exclusive or of array elements
IS_CONTIGUOUS Test contiguity of an array
LBOUND Lower dimension bounds of an array
MATMUL Matrix multiplication
MAXLOC Location of a maximum value in an array
MAXVAL Maximum value in an array
MERGE Merge under mask
MINLOC Location of a minimum value in an array
MINVAL Minimum value in an array
NORM2 L2 norm of an array
PACK Pack an array into an array of rank one under a mask
PARITY True if number of elements in odd
PRODUCT Product of array elements
RESHAPE Reshape an array
SHAPE Shape of an array or scalar
SIZE Total number of elements in an array
SPREAD Replicates array by adding a dimension
SUM Sum of array elements
TRANSPOSE Transpose of an array of rank two
UBOUND Upper dimension bounds of an array
UNPACK Unpack an array of rank one into an array under a mask

As an example, here is how to sum all the positive values in a matrix:

b = sum(A, mask=(A>0.0))

In addition to the array-handling functions, we can also use various other language features to make it easier to work with arrays. For example, say we wanted to convert all the negative elements of a matrix to zero:

where (A<0.0) A = 0.0

Or, say we want to compute the dot product of the 1st and 3rd columns of a matrix:

d = dot_product(A(:,1),A(:,3))

Or, replace every element of a matrix with its sin():

A = sin(A)

That last one is an example of an intrinsic ELEMENTAL function, which is one that can operate on scalars or arrays of any rank. The assignment is done for each element of the array (and can potentially be vectored by the compiler). We can create our own ELEMENTAL functions like so:

elemental function cubic(a,b,c) result(d)
real(wp),intent(in) :: a,b,c
real(wp) :: d
d = a + b**2 + c**3
end function cubic

Note that the inputs and output of this function are scalars. However, since it is ELEMENTAL, it can also be used for arrays like so:

real(wp) :: aa,bb,cc,dd
real(wp),dimension(10,10,10) :: a,b,c,d
!...
dd = cubic(aa,bb,cc) ! works on these
d = cubic(a,b,c) ! also works on these

The array-handling features of Fortran make scientific and engineering programming simpler and less error prone, since the code can be closer to the mathematical expressions, as well as making vectorization by the compiler easier [2]. Here is an example from the Fortran Astrodynamics Toolkit, a vector projection function:

pure function vector_projection(a,b) result(c)

implicit none

real(wp),dimension(:),intent(in) :: a
!! the original vector
real(wp),dimension(size(a)),intent(in) :: b
!! the vector to project on to
real(wp),dimension(size(a)) :: c
!! the projection of a onto b

if (all(a==0.0_wp)) then
    c = 0.0_wp
else
    c = a * dot_product(a,b) / dot_product(a,a)
end if

end function vector_projection

Line 15 is the vector projection formula, written very close to the mathematical notation. The function also works for any size vector. A similar function in C++ would be something like this:

inline std::array<double,3>
vector_projection(const std::array<double,3>& a,
const std::array<double,3>& b)
{
    if (a[0]==0.0 & a[1]==0.0 & a[2]==0.0){
        std::array<double,3> p = {0.0, 0.0, 0.0};
        return p;
    }
    else{
        double aa = a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
        double ab = a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
        double abaa = ab / aa;
        std::array<double,3>
        p = {a[0]*abaa,a[1]*abaa,a[2]*abaa};
        return p;
    }
}

Of course, like anything, there are many other ways to do this in C++ (I'm using C++11 standard library arrays just for fun). Note that the above function only works for vectors with 3 elements. It's left as an exercise for the reader to write one that works for any size vectors like the Fortran one (full disclosure: I mainly just Stack Overflow my way through writing C++ code!) C++'s version of dot_product seems to have a lot more inputs for you to figure out though:

double ab = std::inner_product(a.begin(),a.end(),b.begin(),0.0);

In any event, any sane person using C++ for scientific and engineering work should be using a third-party array library, of which there are many (e.g., Boost, EIGEN, Armadillo, Blitz++), rather than rolling your own like we do at NASA.

References

  1. Walkthrough: Matrix Multiplication [MSDN]
  2. M. Metcalf, The Seven Ages of Fortran, Journal of Computer Science and Technology, Vol. 11 No. 1, April 2011.
  3. My Corner of the World: C++ vs Fortran, September 07, 2011

Apr 03, 2016

C++ vs Fortran (Part 2)

Here is some C++ code (from actual NASA software):

aTilde[ 0] = aTilde[ 1] = aTilde[ 2] =
aTilde[ 3] = aTilde[ 4] = aTilde[ 5] =
aTilde[ 6] = aTilde[ 7] = aTilde[ 8] =
aTilde[ 9] = aTilde[10] = aTilde[11] =
aTilde[12] = aTilde[13] = aTilde[14] =
aTilde[15] = aTilde[16] = aTilde[17] =
aTilde[18] = aTilde[19] = aTilde[20] =
aTilde[21] = aTilde[22] = aTilde[23] =
aTilde[24] = aTilde[25] = aTilde[26] =
aTilde[27] = aTilde[28] = aTilde[29] =
aTilde[30] = aTilde[31] = aTilde[32] =
aTilde[33] = aTilde[34] = aTilde[35] = 0.0;

Here is the same code translated to Fortran:

aTilde = 0.0

😀

References

Jul 17, 2015

MATH77 Library

jpl

It looks like JPL just opensourced their MATH77 Fortran library. MATH77 is a library of Fortran 77 subroutines implementing various numerical algorithms. It was developed over decades at JPL, and contains some very high-quality and time-tested code. The code is released under a BSD-type license. There is also a C version for people who love semicolons.

See also

  1. "MATH77/mathc90 Libraries Available at netlib", July 17, 2015. [comp.lang.fortran]

Jun 21, 2015

Refactoring

refactoring

SLSQP [1-2] is a sequential least squares constrained optimization algorithm written by Dieter Kraft in the 1980s. Today, it forms part of the Python pyOpt package for solving nonlinear constrained optimization problems. It was written in FORTRAN 77, and is filled with such ancient horrors like arithmetic IF statements, GOTO statements (even computed and assigned GOTO statements), statement functions, etc. It also assumes the implicit saving of variables in the LINMIN subroutine (which is totally non-standard, but was done by some early compilers), and in some instances assumes that variables are initialized with a zero value (also non-standard, but also done by some compilers). It has an old-fashioned "reverse communication" style interface, with generous use of "WORK" arrays passed around using assumed-size dummy arguments. The pyOpt users don't seem to mind all this, since they are calling it from a Python wrapper, and I guess it works fine (for now).

I don't know if anyone has ever tried to refactor this into modern Fortran with a good object-oriented interface. It seems like it would be a good idea. I've been toying with it for a while. It is always an interesting process to refactor old Fortran code. The fixed-form to free-form conversion can usually be done fairly easily with an automated script, but untangling the spaghetti can be a bit harder. Consider the following code snipet from the subroutine HFTI:

C DETERMINE PSEUDORANK
    DO 90 j=1,ldiag
90  IF(ABS(a(j,j)).LE.tau) GOTO 100
    k=ldiag
    GOTO 110
100 k=j-1
110 kp1=k+1

A modern version that performs the exact same calculation is:

!determine pseudorank:
do j=1,ldiag
    if (abs(a(j,j))<=tau) exit
end do
k=j-1
kp1=j

As you can see, no line numbers or GOTO statements are required. The new version even eliminates one assignment and one addition operation. Note that the refactored version depends on the fact that the index of the DO loop, if it completes the loop (i.e., if the exit statement is not triggered), will have a value of ldiag+1. This behavior has been part of the Fortran standard for some time (perhaps it was not when this routine was written?).

There is also a C translation of SLSQP available, which fixes some of the issues in the original code. Their version of the above code snippet is a pretty straightforward translation of the FORTRAN 77 version, except with all the ugliness of C on full display (a[j + j * a_dim1], seriously?):

/* DETERMINE PSEUDORANK */
i__2 = ldiag;
for (j = 1; j <= i__2; ++j) {
    /* L90: */
    if ((d__1 = a[j + j * a_dim1], fabs(d__1)) <= *tau) {
        goto L100;
    }
}
k = ldiag;
goto L110;
L100:
k = j - 1;
L110:
kp1 = k + 1;

See Also

  1. Dieter Kraft, "A Software Package for Sequential Quadratic Programming", DFVLR-FB 88-28, 1988.
  2. Dieter Kraft, "Algorithm 733: TOMP–Fortran Modules for Optimal Control Calculations," ACM Transactions on Mathematical Software, Vol. 20, No. 3, p. 262-281 (1994).
  3. John Burkardt, LAWSON Least Squares Routines, 21 October 2008. [Fortran 90 versions of the least squares routines used in SLSQP].

Mar 21, 2015

Too Much Confusion

fortran-wheel

There is a lot of confusion and misinformation about the Fortran programming language on the internet, and a general ignorance about it among programmers. Most younger programmers who use languages invented five minutes ago probably have never seen it, and may only be dimly aware of it as some obsolete language that nobody uses anymore.

You may be surprised to learn that Fortran is a modern, object-oriented, general-purpose programming language. Fortran is not a programming language created by computer scientists to write operating systems, nor does it include every single programming concept anyone's ever heard of. It is a language designed for computational efficiency and efficient array manipulation, and has a clear uncluttered syntax that can be read and understood by non-experts with only a little effort and training. It is particularly suited for numerical and scientific programming by non-expert programmers such as engineers and scientists. Learning all of Fortran is considerably easier than learning all of C++ (as an example). Sure, you can't do template metaprogramming, but few engineer/scientist types would ever want to do that anyway (besides, it's bad for you and could make you go blind).

By "Fortran", I mean modern Fortran (i.e., Fortran 2003 and 2008, which is the latest standard). Yes, the roots of Fortran go way back to the 1950s. Sure, early Fortran programs were written on punched cards. So what? Latin was once scratched into wax tablets, but that isn't really relevant to modern Italian and French speakers. In fact, the Fortran language has evolved considerably since it was first standardized in 1966. It generally has followed a cycle where a major update is followed by a minor update (1977=minor, 1990=major, 1995=minor, 2003=major, 2008=minor). It has been said that the 2003 update was as big an update to Fortran 95 as C++ was to C! Ten years later, the GNU Fortran compiler is still not fully F2003 compliment (the Intel compiler only recently became so).

People who attempt to enter the world of Fortran programming are easily corrupted and discouraged by misinformation. Even a lot of old-school Fortran users are unaware of the later standards. This is too bad, because modern Fortran is actually quite a respectable programming language for a lot of technical applications. This article is a pretty good overview of Fortran for C/C++ programmers. However, it is outdated, since it is confined to Fortran 95. Most of the limitations it mentions (no procedure pointers, clunky character strings, lack of an intent attribute for pointer dummy arguments, the nonstandardness of the ; character) have been rectified in subsequent standards.

The fact is the internet is not really the best source of information for modern Fortran. One day, maybe, there will be vibrant community of Fortran users on the internet, extensive online documentation, open source projects, and all your questions will simply be a web search away (cf., Python). But for now, you'll probably have to buy some books. If a book has the numbers 77, 90, or 95 in the title, don't open it, it will only confuse you. This is not to say that there aren't friendly Fortran folk on the internet who will also help you out. Two of the best places to go with questions are the Intel Fortran forum and the comp.lang.fortran newsgroup (yes, apparently, Usenet still exists).

References

Next → Page 1 of 2