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 -fitting
and equations 2.1 in section 2.1.8
on page . The equations take as their
input arrays *x*_{i}, *y*_{i} and ,
all of length *N*, and return
*a*, ,
*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, QHaving defined the variables we can now commence our computation.

First observe that *all* the -fit quantities
can be expressed in terms of

so we should evaluate these first. We begin by evaluating :

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
,
and `S_sigma`

, which corresponds to
.
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

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

(2.9) |

Because we will have to

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
and ,
in case
something may simplify in the process. Indeed it does:

= | |||

= | |||

= | |||

= | (2.10) |

and similarly

= | |||

= | |||

= | (2.11) |

Now we can easily finish the computation of our -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
and *Q*.
In this section we will do
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

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

endso that both the run time system and the compiler know that this is it

program chiso that

Figure 2.4 presents our first Fortran program in full.

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,
yet.
If the system kindly initialises the latter to zero, as they
often do, the program will crash on the first attempt to
divide
,
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 -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.

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.