Chapter 1 Programming in a high-level language – Fortran 9095

Document Sample
Chapter 1 Programming in a high-level language – Fortran 9095 Powered By Docstoc
					Chapter 1

Programming in a high-level
language – Fortran 90/95

Programming Class 1 (Week 1)
Why learn a high-level language?
The cutting edge in scientific computing is efficient implementation of float-
ing point arithmetic, i.e. arithmetic involving numbers expressed in decimal
expansions. Although integrated systems such as Matlab are highly convenient
you pay a penalty for this in terms of efficiency and memory requirements.
Matlab on its own is not suitable for very large problems (or even medium-sized
ones without careful coding). High-level compiled languages (Fortran, C, C++,
Pascal) are more efficient.

Which high-level language to learn?
Fortran77. “Old fashioned”, for example it does not allow abstract data types
     or recursion. However, since a great deal of well-engineered scientific software
     is in Fortran77, it remains important today. The simplicity of its syntax al-
     lows compilers to highly optimise the resulting machine code making full use
     of the cache, the pipeline and other hardware features in modern computers.
Fortran90/95. modern upgrade of Fortran77 which supports many important fea-
     tures missing from Fortran77. It is backward compatible with Fortran77.
     (Fortran2003 standard coming!)

C and C++. Many scientific programmers also use languages like C (especially)

     or C++. In terms of efficiency comparable to Fortran.
Pascal. “Old fashioned” like Fortran77.
Java. Not as efficient for scientific applications.
There is a growing interest in object-oriented programming languages (like C++
or Java) for scientific applications as opposed to procedural languages (like Fortran).
Modern C++ compilers allow for similar efficiency as Fortran77, while allowing
for data hiding, code abstraction and other modern developments in computer
science. (Fortran2003 standard includes more object-orientation!)
    Many workers in scientific computing will need to work with multi-language
codes (and indeed it is nowadays quite easy to couple Fortran and C or C++ code).
However, for us it is more important to learn the principles of software design,
then to learn a particular language. In this course we will use Fortran95 and we
will spend the first 2 weeks learning the basics of this language.

Example 1
Assemble the n × n Hilbert matrix
                       Hi,j = 1/(i + j − 1) , i, j = 1, . . . n .
Programs given:
   • Fortran95 program in week1/hilbert.f90
   • Three Matlab programs in week1: hilbert.m, hilbert1.m and hilbert2.m

 Lab Exercise 1. Connect to the computers in the MSc room and edit these files.

    1. Start → Programs → Unix Applications → X-Windows Manager

    2. On grey area right click, select any BUCS machine. Logon.

    3. Type
                      ssh -l<userid> mapc-desknode-0**.maths

    4. Copy my files over to your file space:

                        cp ˜ jvl20/ma50177/week1/hilbert* .

    5. Type emacs& to start the Emacs-editor. Use the file button to open

The Fortran95 source code in week1/hilbert.f90:

program hilbert                  ! optional program label

! f95 program to assemble the Hilbert matrix
! We begin by declaring our variables

real, allocatable, dimension(:,:) :: H    ! declares a real matrix
                                          ! H with two dimensions
                                   ! (i.e. rank 2). The extent
                                   ! of each dimension is declared
                                   ! in the allocate command below

real :: t1,t2                      ! real numbers used for timing
                                   ! this program
integer :: i,j,n                   ! integers, used for counting

 print*,’Type n’                   !   message for the user
 read*, n                          !   extent of the dimensions of H
 allocate(H(n,n))                  !   allocate storage
 call cpu_time(t1)                 !   look at the clock

 do i = 1,n                        ! nested loop for assembling H
   do j = 1,n
    H(i,j) = 1.0/(i+j-1.0)         ! mixed mode expression for H(i,j)
   end do
 end do

 call cpu_time(t2)

 print*,’n=’,n,’cputime taken = ’,t2-t1 ! prints out timing results

 print*, H(1,2)                    ! prints out (1,2)th element of H

end program hilbert              ! mandatory end statement

To compile this program type ifort hilbert.f90 to the $ prompt. You will
get an executable a.out .

Type a.out to run this executable.

Notes on hilbert.f90:

  1. “Dynamic” allocation of array storage for H.

  2. All variables must be declared!

  3. Study the syntax of the do loop.

  4. Real and integer data types appear – different in Fortran95 !

  5. “Mixed mode” arithmetic in the formula for H(i,j) ([MeRe p.36], [Sm p.16]).

  6. Use of intrinsic function cpu_time.

  7. Only H(1,2) printed to screen, as example.

 Lab Exercise 2. What is the effect of the following changes to hilbert.f90?

    1. Change H to h in the last print statement.

    2. Comment out the line: real,allocatable...
       (use ! to comment out a line).

    3. Comment out the line: integer ::      i,j,n

    4. Change the line integer ::       i,j,n to integer ::      i,j,n,k

 Make each of the changes, recompile and rerun to see the result. What do you
 learn from this?

 Now, repeat these experiments using the compile command: ifort -u .

 Try some other changes.

To illustrate the efficiency of a compiled language, compare the program hilbert.f90
with two Matlab programs which carry out the same task:

% program hilbert.m : script to assemble the n x n Hilbert matrix
% using a double loop

n = input(’type n:   ’)
H = zeros(n,n);             % Initialise H as an n x n matrix

tic                         % look at clock

for j = 1:n                 % Use a double do loop
  for i = 1:n
    H(i,j) = 1/(i+j-1);         % doesn’t matter whether we
                                % write 1 as 1.0 or not
   end % end of j loop
end % end of i loop



% end of program hilbert.m

% program hilbert2.m : script to assemble the n x n Hilbert matrix
% avoiding loops by using array manipulation facilities of Matlab

n = input(’type n:    ’)
H = zeros(n,n);    % Initialise H

H = ones(n,n)./((1:n)’*ones(1,n) + ones(n,1)*(1:n) - ones(n,n));


% end of program hilbert2.m

Efficiency: CPUtime in seconds (2.1GHz Athlon, 3.7 GByte RAM)

        Matlab 6 (Release 13)         Matlab 7.3 (Release 2006b)
 n      hilbert.m   hilbert2.m    hilbert.m    hilbert1.m    hilbert2.m    hilbert.f90
1000       5.4         0.29          0.038        0.067         0.14          0.018
2000      21.6         1.13           0.15         0.61         0.53          0.075
4000      88.6          4.5           0.60         4.2           2.1          0.30
8000      ****         ****           2.4          17.1          8.4           1.2
16000     ****         ****          ****         ****          ****           4.9

(**** = out of memory (3.7 GByte RAM))

In Matlab 6 (and in any older version) hilbert.m is almost 300 times slower
than the compiled code hilbert.f90 and hilbert2.m is still 15 times slower. In
Matlab 7 the times are closer, but the compiled code is still more than twice as

    The high cost of hilbert.m (in Matlab 6) is due to the fact that Matlab is an
interpreted language: each line of code is translated to machine instructions each
time it is reached: here the same line is retranslated n2 times. We can see that in
the later versions of Matlab (since Release 14, September 2004) the interpreter is
able to identify the loop and to avoid retranslating the same line n2 times. However,
this will not always be possible when a for-loop becomes more complicated. The
faster code hilbert2.m relies on Matlab built-in matrix manipulation functions
to avoid loops. This built-in functions are themselves Fortran routines from fast
libraries like BLAS and LAPACK which we will come across later in this course.

   A further problem with interpreted languages like Matlab is memory. Even
on my computer, which has 3.7 GByte memory, I can only go with n up to 4000 or
8000 for Matlab 6 and Matlab 7.3 respectively, whereas in hilbert.f90 I can go
up to 20000!

   The bottom line is that Matlab has problems to compete for large problems.
This will become even more apparent for more complicated problems.

Lab Exercise 3. Complete the two tables on any of the linux machines in the
MSc room (3.0Ghz Intel, 1 GByte memory, Matlab 7 (Release 14)).

          n     hilbert.m   hilbert1.m hilbert2.m     hilbert.f90

Programming Class 2 (Week 1)
IMPORTANT: You should read this already before your Thursday class. The
exercises will be done in the class. Please come prepared to ask questions.

This class will take place in the MSc room. The purpose is:
   • to learn some Emacs skills
   • to learn about elementary arithmetic, data types, and intrinsic functions in
   • to learn about control of flow in Fortran95 programs
   • to translate an elementary MATLAB program into Fortran95
By the end of Friday’s class you should have done all the exercises in this section.

The Emacs editor
Emacs is the recommended editor for developing Fortran95 code. To start Emacs
type emacs& to the $ prompt. You can use the menu bar at the top to load, insert
and save files.

You can use Emacs as a file manager. Under the ‘file” button, select “Open direc-
tory”. You can then navigate around your files and if you like you can use Emacs
to create and delete directories and files.

To get pretty colours for your Fortran code, go to the Options menu and select
Syntax Highlighting. Then click on Save Options.

Compiling and running programs
Exercise 1. A simple program for numerical integration.
Suppose f is a given function on an interval [α, β]. Then Simpson’s rule for com-
puting an approximation to the exact integral α f (x)dx is:
                       Iα,β (f ) =     {f (α) + 4f (γ) + f (β)} ,
where h = (β − α) and γ = (α + β)/2 .
A program to implement this in a very simple way is at


Copy this program to your filespace. Compile and run it. What does the result tell
you? Replace f (x) = x3 by f (x) = x5 . What do you find?

Intrinsic functions and data types
Exercise 2. Intrinsic functions.
Look up the list of intrinsic functions in the Fortran95 book of Smith [Sm] or Met-
calf and Reid [MR]. Find the entry for atan. Using the fact that tan(π/4) = 1, write
a line of code to compute π and then compute an approximation to 0 π cos(πx)dx
using Simpson’s rule. What is the error in the result?

Exercise 3. Data types.
Fortran95 allows reals of both single and double precision. You can specify the
precision by specifying a kind parameter:

   • real (kind=4)      ... a single precision real
   • real (kind=8)      ... a double precision real

Note that if you do not specify the kind parameter, single precision is used by
Also, if a program requires the real constant 10.0 in double precision, it should be
typed as 10.0_8. Otherwise it will only be of single precision and might flaw any
double precision calculations.
Edit your s_simpson.f90 so that the computations are done in double precision.
Compile and run again. How can you be sure that the double precision is working?

Hint: To avoid having to edit the declaration of every variable when you want to
change the precision of a program, you can simply declare the kind as a parameter
and then you only need to change it once. So in the declaration of variables you add
something like:
                          integer, parameter :: rk=1

Then when declaring a real variable integral you say

                          real (kind=rk) :: integral

Also, if a program requires the real constant 10.0, it can be typed as 10.0_rk to
obtain the required precision. Then when editing you only need to change the value
in the parameter declaration once.

Warning: Fortran95 is a weakly typed language. Variables that start with
i,j,k,l,m or n are understood implictly to be integers unless declared otherwise.
This is a rather dangerous Fortran77 legacy. To override it you can put

                                 implicit none

at the top of each program unit. Alternatively, compile with the -u option. Recom-
mendation: always compile with the -u option.

Writing your first program
Exercise 4. Control of Flow.
As you will remember in Matlab , the most common constructions for controlling
the flow in a program are the

  (i) for loop,

 (ii) while loop,

(iii) if statement.

In Fortran95 the counterpart of (i) is a do loop. You have seen an example of
this already in hilbert.f90 . For (iii) there is an if statement which is similar
to Matlab but with slightly different syntax. For (ii), there is no explicit while
loop. However, the required action can be achieved with a do loop combined with
the if and exit statements. For example, to loop over all values of i from 1 to n,
stopping when i=j we can write

       do i = 1,n
         if (i==j) exit
       end do

The exit command causes exit from the do loop and the program continues with
the command after the end do statement. Look in the Fortran95 books for more

What to do: Copy over the program ~jvl20/ma50177/week1/macheps.m. This
is a Matlab program to compute the machine epsilon for a given machine. You
should have seen this before in MA50174. The machine epsilon is the smallest
positive number ǫ such that the machine representation of 1 + ǫ differs from the
machine representation of 1. Translate this program from Matlab to Fortran95.
Compile and run your Fortran95 code and compare its output to the Matlab code.

IMPORTANT: Please try to complete these exercises before Friday’s class.

Programming Class 3 (Week 1)


 Main Program and Subprograms
 Subprograms (Procedures):

               function and subroutine

 The argument list, intent of an argument
 Interface blocks
 The Makefile
 Examples: various quadrature procedures

Example: Numerical Integration
Simpson’s rule for approximating      α
                                          f (x)dx is
                      Iα,β (f ) =     {f (α) + 4f (γ) + f (β)}

where h = (β − α) and γ = (α + β)/2 .
Task 1. Open the file s_simpson.f90 in emacs (see also below).

program s_simpson

! simple program to approximate the integral of a function f over the
! interval [alpha, beta] using Simpson’s rule. Here f(x) = x**3

real :: alpha, beta, gamma
real :: f_alpha, f_beta, f_gamma, h
real :: integral

print*, "Type alpha,beta"
read*, alpha,beta

gamma = (alpha + beta)/2

h = (beta-alpha)

f_alpha = alpha**3
f_beta = beta**3
f_gamma = gamma**3

integral = (h/6)*(f_alpha + 4*f_gamma + f_beta)

print*, integral

end program s_simpson

It works but is not very convenient/flexible:
  1. To change the function f we have to edit 3 lines and recompile.
  2. We often will want to apply the rule over many small subintervals of a
     larger interval. This is the “composite Simpson’s rule”. To do this we would
     have to rerun the current program many times, type in new data each
     time, and add up the answers afterwards.

To fix Problem 1 we put the definition of f in a function subprogram:

  function func (x) result (func result)
  real, intent(in) :: x
  real :: func result
  func result = x**3
  end function func

This is stored as func.f90 (in week1/simpson).

The intent statement is optional here. It is good programming practice to put
it in. x is a dummy argument in func. If a dummy argument is specified as
intent(in) then it cannot be changed inside the function.

Now in s_simpson we can write
  f_alpha = func(alpha)
Then alpha is the actual argument, which is passed to func.

To address Problem 2, it would be better to rewrite the whole Simpson procedure
as a subprogram, which takes as arguments the function func and the end points
alpha, beta and produces the approximate integral. This could then be called
many times from a main program, with different values of the arguments, if
required. The best way to do this is using a subroutine.

Unlike a function, a subroutine often has both input and output arguments,
or even arguments which have intent(inout). It has no “result” as such. All data
transfer with the calling program is via the arguments. Also, unlike a function, a
subroutine is invoked with the call command.
subroutine simpson(integral,func,alpha,beta)

! simple program to approximate the integral of a function func
! over the interval [alpha, beta] using Simpson’s rule

 real, intent(out) :: integral           ! The result of Simpson’s rule
 real, intent(in) :: alpha, beta         ! The inputs
 real :: gamma, h                        ! variables needed internally

! Interface block to ensure consistency and detect errors in passing
! arguments to subprograms essentially to let the subroutine know
! what type of thing func is

   function func(x) result(func_result)
   real, intent(in) :: x
   real :: func_result
   end function func
 end interface

 gamma = (alpha + beta)/2
 h = (beta-alpha)
 integral = (h/6)*(func(alpha) + 4*func(gamma) + func(beta))

end subroutine simpson

This subroutine is called from a main program, in our case integration.f90

All these files are located in week1/simpson.

Task 2. Copy the entire directory week1/simpson by typing:

           cp -R ~jvl20/ma50177/week1/simpson .          (‘‘-R’’ for recursive)

Task 3. Compile with the command

              ifort -u func.f90 simpson.f90 integration.f90

Code development: The Makefile
Since we now have three program units which can be developed separately it is best
to control their development using a Makefile. This contains all the commands
which are needed to compile the program units and then make an executable. When
you have a Makefile, type the command make to run it. (Only program units which
have been edited/changed since the last time you ran make are recompiled.)

A suitable Makefile to compile and link all the components of our code in the files
func.f90 simpson.f90 integration.f90 is in week1/simpson/Makefile:

# This is a comment line in the Makefile

OPT = -u      # all ifort calls use the -u option
              # (suppresses implicit typing)

all: integration           # forces all the commands needed to create
                           # ‘‘integration’’ to be executed,
                           # unless there is nothing to do

func.o: func.f90                     #   compiles func.f90. The -c option
      ifort $(OPT) -c func.f90       #   causes the production of an object
                                     #   file only which will be linked
                                     #   with other object files below

# then does the same for simpson.f90 and integration.f90

simpson.o: simpson.f90
      ifort $(OPT) -c simpson.f90

integration.o: integration.f90
      ifort $(OPT) -c integration.f90

# now link the three object files to get an executable
# -- in this case called integration

integration: integration.o simpson.o func.o
      ifort -o integration $(OPT) integration.o simpson.o func.o

         rm -rf *.o    core    # creates a command to remove old .o and
                               # core files

# type ‘‘make clean’’ to run this

Task 4. Go to the directory simpson in your filespace and type
                           make clean and then make

Note. Using ifort -o testx test.f90 creates and executable called testx
(rather than a.out).

Passing parameters
Let us return to the example of approximating
                                          π cos(πx)dx.
To implement Simpson’s rule we can edit the file func.f90 to compute π cos(πx).
However we need to evaluate π. It would be wasteful to have the function func
evaluate π every time it is called. It is better to evaluate π once (in the main
program) and then pass it down, first to the simpson subroutine and then to func.
A suitable implementation of this is in week1/simpson2

Task 5. Copy over the whole directory with the command:
                      cp -r ~jvl20/ma50177/week1/simpson2 .
Have a look at the files simpson2.f90, func2.f90, integration2.f90

Lab Exercise 1
Create a suitable Makefile in your directory simpson2. Check that it works.

Lab Exercise 2 (Composite Simpson’s rule)
The 3 point Simpson’s rule is not accurate enough for most integation tasks. To im-
prove accuracy, the interval of integration [α, β] may be divided into N subintervals
[xi−1 , xi ], where
                       xi = α + i(β − α)/n, i = 1, · · · , N .
Write a subroutine csimpson with input arguments func, α, β and N which
implements this, by calling simpson on each subinterval.

Test it by integrating func(x) = ex , over [0, 1], and by checking that the approx-
imation converges to the exact integral as N → ∞.