Fortran

Although in popular lore there are various Fortran dialects, such as the ancient Fortran-IV , which was then followed by Fortran-66  and Fortran-77 , to be eventually superseded by Fortran-8x  (which never materialised), and, at long last by Fortran-90  and, more recently, by Fortran-95 , officially there is only one Fortran, the ANSI Fortran . The language, as currently defined by the ANSI Committee, is closest to what people on the streets call Fortran-95. What many good housewives call Fortran-77, is ANSI Fortran's official subset with a number of deprecated features. These, it is said, may even become extirpated one day, but, until that day of deliverance, which inspired prophets of old foretold and many a prophet of today foretell still, they shall continue to be with us, and you are welcome to use them in your programs, as you are welcome to drink alcohol, smoke tobacco, chew gum (though not in Singapore ), and make indecent proposals to consenting adults.

Let us now return to the problem at hand, that is to $\chi ^2$-fitting and equations 2.1 in section 2.1.8 on page . The equations take as their input arrays x_i, y_i and $\sigma _i$, all of length N, and return a, $\sigma_a$, b, $\sigma_b$, and the goodness of fit, Q.

In order to convert the equations into a Fortran program we must begin by declaring the corresponding Fortran variables in the preamble of the program. The arrays  can be statically allocated or dynamically allocated. In the latter case, we don't have to specify their length in advance. In the former, we do. Fortran also allows us to specify any precision we may wish for, within certain limits. The precision  can be specified in many ways. For example we can use the selected_real_kind function, which takes two arguments: number of significant digits we wish for, and the highest expected exponent (the power of 10), and returns the smallest size of a real on a given architecture that will satisfy our requirements.

integer, parameter :: long = selected_real_kind(9,99)
integer, parameter :: n = 100
real(kind=long), dimension(n) :: x, y, sigma
real(kind=long) :: a, b, sigma_a, sigma_b, Q

Having defined the variables we can now commence our computation.

First observe that all the $\chi ^2$-fit quantities can be expressed in terms of

$$\begin{eqnarray*}S_\sigma &=& \sum_{i=1}^N \frac{1}{\sigma_i^2}, \\ S_x &=& \su... ...ma_i^2}, \qquad S_{xy} = \sum_{i=1}^N \frac{x_i y_i}{\sigma_i^2} \end{eqnarray*}$$

so we should evaluate these first. We begin by evaluating $S_\sigma$:

one_by_sigma_2 = 1 / sigma**2
S_sigma = sum(one_by_sigma_2)

In the code above we have used two additional variables: one_by_sigma_2, which is an array that stores $1/\sigma_i^2$, and S_sigma, which corresponds to $\sum_i 1/\sigma_i^2$.

On seeing this part of the code Fortran will first recognise that sigma is an array. 1/sigma**2 means: take each element of the array and square it, and then take its inverse. The resulting floating point number, of type real(kind=long), will be assigned to the corresponding entry in the one_by_sigma_2 array, that is to a slot with the same index number.

On a parallel  or on a vector machine (but not necessarily on a multicomputer) this computation should be performed in parallel without any special requests from your side.

The second line of the code takes the one_by_sigma_2 array, sums all its component, and writes the result on S_sigma. Again, special hardware, which some supercomputers or even more advanced microprocessors may have, may be employed in order to evaluate this sum very fast, possibly in parallel too.

The definitions of the new variables introduced in this part of the code should be added to the preamble of the program. We will make use of this opportunity and also define the remaining S parameters:

real(kind=long), dimension(n) :: one_by_sigma_2
real(kind=long) :: S_sigma, S_x, S_y, S_xx, S_xy

Having evaluated one_by_sigma_2, we can easily compute the remaining parameters:

S_x = sum(x * one_by_sigma_2)
S_y = sum(y * one_by_sigma_2)
S_xx = sum(x**2 * one_by_sigma_2)
S_xy = sum(x * y * one_by_sigma_2)

An alternative solution would be to do this as follows:

S_sigma = sum (1 / sigma**2)
S_x = sum(x / sigma)
S_y = sum(y / sigma)
...

and so on. But the resulting program would be a lot slower, even though, it would be somewhat more concise. The reason why it would be so much slower  is because divisions, and especially floating point divisions, are extremely expensive on just about any hardware that I know of. So, if in your program you have to divide by the same quantity very frequently, it pays to evaluate its inverse beforehand and then replace all divisions with multiplications, as we have done in our first version of the code.

The next thing to observe, while scrutinising the $\chi ^2$-fit equations is that all final expressions are divided by the same quantity, which is often referred to as $\Delta$:

$$\Delta = S_\sigma S_{xx} - \left( S_x \right)^2$$ (2.9)

Because we will have to divide by this quantity many times, again, we are more interested in its inverse $1/\Delta$:

real(kind=long) :: one_by_delta
...
one_by_delta = 1 / (S_sigma * S_xx - S_x**2)

Before going any further let us still do some pencil (or, say, Maxima) work on our expressions for $\sigma_a$ and $\sigma_b$, in case something may simplify in the process. Indeed it does:

$$\sigma_a^2 \sum_{i=1}^N \sigma_i^2 \left(\frac{S_{xx} - x_i S_x} {\sigma_i^2... ...a^2} \sum_{i=1}^N \frac{S_{xx}^2 - 2 S_{xx} x_i S_x + x_i^2 S_x^2} {\sigma_i^2}$$ =

$$\frac{1}{\Delta^2} \left( S_{xx}^2\sum_{i=1}^N \frac{1}{\sigma_i^... ...^N \frac{x_i}{\sigma_i^2} + S_x^2 \sum_{i=1}^N \frac{x_i^2}{\sigma_i^2} \right)$$ =

$$\frac{1}{\Delta^2} \left( S_{xx}^2 S_\sigma - 2 S_{xx} S_x S_x + S_x^2 S_{xx} \right) = \frac{S_{xx}}{\Delta^2}\left(S_{xx}S_\sigma - S_x^2\right)$$ =

$$\frac{S_{xx}}{\Delta^2}\Delta = \frac{S_{xx}}{\Delta}$$ = (2.10)

and similarly

$$\sigma_b^2 \sum_{i=1}^N \sigma_i^2 \left(\frac{x_i S_\sigma - S_x} {\sigma_i... ... \sum_{i=1}^N \frac{S_\sigma^2 x_i^2 - 2 S_\sigma x_i S_x + S_x^2} {\sigma_i^2}$$ =

$$\frac{1}{\Delta^2} \left(S_\sigma^2 S_{xx} - 2 S_\sigma S_x^2 + S_x^2 S_\sigma \right)$$ =

$$\frac{S_\sigma}{\Delta^2} \left(S_\sigma S_{xx} - S_x^2\right) = \frac{S_\sigma}{\Delta}$$ = (2.11)

Now we can easily finish the computation of our $\chi ^2$-fit parameters:

a = one_by_delta * (S_xx * S_y - S_x * S_xy)
b = one_by_delta * (S_sigma * S_xy - S_x * S_y)
sigma_a = sqrt(S_xx * one_by_delta)
sigma_b = sqrt(S_sigma * one_by_delta)

This is not the end yet. We must evaluate $\chi ^2$ and Q. In this section we will do $\chi ^2$ only, because Q needs a special section for itself.

real(kind=long) :: chi_2
...
chi_2 = sum(one_by_sigma_2 * (y - a - b * x)**2)

Whereas it is pretty obvious what multiplication of an array by a scalar should yield, it is less obvious what should be an outcome of adding a scalar to an array (at least less obvious to a mathematician, because such an operation is not defined in mathematics). Fortran in this context will simply perform the addition on every element of the array. In other words, it is as if the scalar has been upgraded to an array.

It is rather polite  to end your Fortran program with the word

end

so that both the run time system and the compiler know that this is it^2.6. It is also polite to begin the program with a declarative statement such as

program chi

so that you and your colleagues know how to refer to this program.

Figure 2.4 presents our first Fortran program in full.

  

Figure: Our first attempt at a $\chi ^2$ fitting program

Here you have encountered for the first time  the Fortran line continuation character &. In Fortran newlines are statement terminators. If a statement doesn't fit on a single line, the line must be explicitly continued, and that's done by placing & at the end of the line.

On a Sun workstation running Solaris 2.6, you would compile  this program as follows:

gustav@blanc:../src 16:16:13 !550 $ f90 -o chi chi.f90
gustav@blanc:../src 16:16:22 !551 $

whereas on Silicon Graphics workstation you could try:

gustav@fesarius:../src 16:32:15 !86 $ f90 -n32 -o chi chi.f90
gustav@fesarius:../src 16:32:28 !87 $

The -n32 option means that we are compiling for a 32-bit architecture. This is the only architecture supported on the Ships cluster. The other architecture supported by the MIPS f90 compiler is -64:

gustav@fesarius:../src 16:33:39 !89 $ f90 -64 -o chi chi.f90
ld64: FATAL 9: I/O error (/usr/lib64/crt1.o): No such file or directory
f90 ERROR:  /usr/lib64/cmplrs/ld64 returned non-zero status 32
gustav@fesarius:../src 16:33:51 !90 $

Well, if it's not there, it's not there.

The file name extension must be .f90. Otherwise the f90 compiler will assume that this is a Fortran-77 file and will spew great many error messages.

I don't advise running this program at this stage, because we haven't initialised the arrays x_i, y_i and, even more importantly, $\sigma _i$ yet. If the system kindly initialises the latter to zero, as they often do, the program will crash on the first attempt to divide $1/\sigma_1 = 1/0$, sic!

But even at this stage we can draw a number of important lessons from this exercise.

1.

The body of the program is only one line longer than the number of equations that it implements.

2.

The equations map almost verbatim onto the body of the program. We have, in fact, named our variables in a way that makes that mapping more obvious.

3.

Array operations indicate explicit parallelism.

4.

The program isn't cluttered by DO or for loops, as it would invariably be in C, C++, Fortran-77, Common Lisp and Scheme. Array indices aren't explicitly invoked anywhere in the program either.

5.

Only one auxiliary variable, one_by_sigma_2 was needed to map our $\chi ^2$-fit equations onto Fortran. We could even do without it, but the program would run slower.

6.

Fortran notation is austere: there are no baroque decorations such as braces, begin .. end brackets, etc. Even the word program is not mandatory (but the terminating command end is!).

7.

The program is entirely portable. It does not depend on includes or any system specific features. There is nothing in it that would relate explicitly to an operating system or to a hardware architecture. Even the precision request is architecture independent.

In short, Fortran is a language specifically designed for science. It is a very elegant language, which is a pleasure to use on large scientific codes. Its evolution was long and arduous, and its beginnings modest. Present day Fortran offers nearly every features found in other languages, e.g., C and C++, and more. High Performance Fortran, which is closely related to ANSI Fortran provides the easiest way to write parallel programs. This is due, in part, to the fact that plain ANSI Fortran is already an explicitly parallel language. It is therefore relatively easy to distribute Fortran array operations over a parallel computer.

These comments, do not apply to Fortran's older brother, Fortran-77. More often than not, a specific advice about parallelism must be passed to a Fortran-77 compiler by a programmer.