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 xi, yi 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
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
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
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
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
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 -fit equations is
that all final expressions are divided by the same quantity,
which is often referred to as :
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
something may simplify in the process. Indeed it does:
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 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
endso that both the run time system and the compiler know that this is it2.6. It is also polite to begin the program with a declarative statement such as
program chiso that you and your colleagues know how to refer to this program.
Figure 2.4 presents our first Fortran program in full.
Here you have encountered for the first time
the Fortran line continuation
&. 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
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
-n32option 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
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 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 xi, yi 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.
forloops, as it would invariably be in
C++, Fortran-77, Common Lisp and Scheme. Array indices aren't explicitly invoked anywhere in the program either.
one_by_sigma_2was needed to map our -fit equations onto Fortran. We could even do without it, but the program would run slower.
begin .. endbrackets, etc. Even the word
programis not mandatory (but the terminating command
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.