Earth to Mars Trajectory

From: "Jonathan Thornburg [remove -animal to reply]"
Subject: Earth to Mars Trajectory
Newsgroups: sci.astro.research
References:


Kent Paul Dolan wrote:
[[about an orbital mechanics simulation]]
unless you are tied to BLAS or
something similar in the way of a numerical software
library, why not a more modern programming language
that would make your software accessible to a larger
audience?

Fortran is dragging a _lot_ of baggage these days.

This does depend rather strongly on where in the scientific-computing
ecosystem you look. For example, in the Cactus toolkit for the parallel
solution of hyperbolic PDEs (http://www.cactuscode.org), newly-written
physics code is something like
50% Fortran 90
30% C
15% C++
5% Fortran 77
[I emphasize that this is "physics" code; the "infrastructure" code for
MPI, parallel I/O, memory management, etc, is mainly C/C++.]

Partly this is historical inertia, but partly it's that Fortran (still)
has much better support than C/C++/Java for certain kinds of numerical
computation. For example, consider the following trivial Fortran 77
subroutine to multiply two n*n matrices:

subroutine matrix_multiply(n, A, B, AB)
integer n
double precision A(n,n), B(n,n), AB(n,n)
integer i,j,k
double precision sum
do 300 j = 1, n
do 200 i = 1, n
sum = 0.0d0
do 100 k = 1, n
sum = sum + A(i,k)*B(k,j)
100 continue
AB(i,j) = sum
200 continue
300 continue
return
end

Until the C99 standard (which isn't yet widely used), you couldn't
write equivalent code in C without significant uglyness (measured in
programmer (in)convenience, (in)efficiency, or both).
[N.b. "equivalent" means n isn't known at compile time,
and the matrices are stored contiguously in memory. We
can ignore row- versus column-major storage.]
C99 has features ("variable length arrays") which let you write this
subroutine in a natural & efficient manner, but alas these are not
described in most C textbooks.

And as for C++, well, it's considerably *worse*: if you want a useful
matrix class, you have to roll your own. And C99 variable length arrays
aren't in the current ANSI/ISO C++ standard.

For these and othe reasons, the vast majority of high-quality numerical
software libraries (eg., the repository at http://www.netlib.org) are
still written in Fortran, most commonly Fortran 77. A classic example
is ODE integration. IMHO the three best general-purpose ODE integration
codes around are
* Shampine & Gordon's ODE/STEP (written in Fortran 66)
* Hindmarsh's ODEPACK/LSODE (written in Fortran 77)
* Brankin, Gladwell, & Shampine's RKSUITE (written in Fortran 77)

Some newer packages now provide C or C++ wrappers or translations
(eg RKSUITE now includes a 3rd-party translation to C++), but the
primary documentation is still usually written using Fortran
terminology, and "advanced features" are often only available through
the Fortran API.

I conclude that at least a reading knowledge of Fortran remains very
useful for doing scientific computation.

Fortunately, interlanguage calling is fairly easy. In particular, it's
easy (albeit usually slightly operating-system--dependent) to call
Fortran subroutines from a {C,C++,Java,Ada,Pascal,Modula-n,PL/1,...}
program. This is what I do in my own research codes (which are almost
entirely C++ or Perl.)

--
-- From: "Jonathan Thornburg [remove -animal to reply]"
Dept of Astronomy, Indiana University, Bloomington, Indiana, USA
"Washing one's hands of the conflict between the powerful and the
powerless means to side with the powerful, not to be neutral."
-- quote by Freire / poster by Oxfam

In article , Jonathan
Thornburg writes:

Fortran is dragging a _lot_ of baggage these days.

This does depend rather strongly on where in the scientific-computing
ecosystem you look. For example, in the Cactus toolkit for the parallel
solution of hyperbolic PDEs (http://www.cactuscode.org), newly-written
physics code is something like
50% Fortran 90
30% C
15% C++
5% Fortran 77

Indeed. Weather/climate modelling is another area where Fortran is
still big. I worked in this field a bit several years ago. When I
heard that some special users have more disk space if they need more
than 100 MB, I thought "that's not so much" (even though disks were much
more expensive 15 years ago) until I realised they were talking about
100 MB of source code. Also, since I knew that these folks had a lot of
computing power, I didn't think that a couple of days on a Cray was that
much computing time, until I realised it was just the time needed for
compiling the code.

In short, if you are doing serious number crunching, you need the
efficiency of Fortran.

Note that the standard is still evolving. The last update was HUGE.

Partly this is historical inertia, but partly it's that Fortran (still)
has much better support than C/C++/Java for certain kinds of numerical
computation. For example, consider the following trivial Fortran 77
subroutine to multiply two n*n matrices:

[code deleted]

Consider that Fortran 90 has the built-in MATMUL intrinsic function:

FORTRAN

Intrinsic_Procedures

MATMUL

MATMUL (matrix-a, matrix-b)

Class: Transformational function - Generic

Performs matrix multiplication of numeric or logical matrices.

The "matrix"s can be arrays of rank one or two. At least one
argument must be rank two. The size of the first (or only)
dimension of "matrix-b" must equal the last (or only) dimension of
"matrix-a".

The type of the resulting array depends on the data types of the
arguments. The rank and shape of the result follows:

o If "matrix-a" has shape (n,m) and "matrix-b" has shape (m,k),
the result is a rank-two array with shape (n,k).

o If "matrix-a" has shape (m) and "matrix-b" has shape (m,k), the
result is a rank-one array with shape (k).

o If "matrix-a" has shape (n,m) and "matrix-b" has shape (m), the
result is a rank-one array with shape (n).


Examples:

Consider the following:

A is the matrix |2 3 4|, B is the matrix |2 3|,
|3 4 5| |3 4|
|4 5|

X is vector (1, 2), and Y is vector (1, 2, 3).


The result of MATMUL (A, B) is the matrix-matrix product AB with
the value

|29 38|
|38 50|

The result of MATMUL (X, A) is the vector-matrix product XA with
the value (8, 11, 14).

The result of MATMUL (A, Y) is the matrix-vector product AY with
the value (20, 26).

Note that all you need to code is "MATMUL(A,B)".

Built-in functions like SIN, SQRT etc will take an array as an argument
and perform the funtion on every element. I could go on and on. For
serious heavy-duty numerical computation, Fortran is the only serious
choice, both for efficiency of the code and for producing readable code.
(The latter was a disadvantage of Fortran77 and even older versions, but
even Fortran77 was more readable than C.)

Getting back to astronomy, IIRC a lot of code for processing data from
artificial satellites, space probes etc is in Fortran, as is much code
for modelling galaxy formation etc.