SS 4068 Fortran Programming

SS 4068 Fortran Programming

C. H. Patterson

Hilary Term 2008

Recommended text:

Computing for Scientists: Principles of Programming with Fortran 90 and C++

Barlow and Barnett

Variables and Operators

Data Structure

Control

Subprograms: Functions and Subroutines

Characters and Strings

Pointers

Input and Output

Chapter 2 Variables and Operators

Fortran Programme Syntax

Programmesbegin with the PROGRAM statement e.g.

PROGRAM HELLO_WORLD

Programmesend with

STOP

END

Programme statements begin in the 7th column and end without a semicolon in the 72nd column

There is usually only one statement per line but f90 statements on the same line can be separated by ;

Lines are continued with & in the 6th column of the following line

In f90 the comment symbol is ! It replaces // in C

! THIS IS A COMMENT

Variable names can be up to 31 characters in length

Begin with a letter

Not case sensitive

There are no reserved words

Write to screen with the PRINT statement

PRINT *, x

writes the variable x to the screen. * is otherwise the unit number of the file to be written to.

Read from Keyboard with the READ statement

READ *, y

Reads variable y from the keyboard.

Variable Declaration

Fortran allows for variables to be used without declaration and uses the convention that variables with names beginning a-h or o-z are REAL variables while those beginning with i-n are INTEGER. This language feature can be disabled by typing

IMPLICIT NONE

after the PROGRAM statement. This avoids many possible programme bugs. Variables are declared in C++ and f90 by separating variables from the variable type as follows:

float x_data;REAL :: x_data

int mindex;INTEGER:: height

Variables can be initialised at the same time as they are declared

Floaty_data = 0.3;REAL:: y_data = 0.3

In C++ there are short, int and long integer variable types. In Fortran the KIND distinguishes them. (See next section).

In f90attributes of variables are specified after the variable type

INTEGER, PARAMETER :: m = 100 ! m is a parameter, value 100

Note position of ::

In order to declare an array which is to be dynamically allocated, it is declared with the attribute ALLOCATABLE

REAL, DIMENSION(:), ALLOCATABLE :: zdata

READ *,n

ALLOCATE (zdata(1:n))

READ *, (zdata(i),i=1,n)

…

DEALLOCATE(zdata)

The code above declares a 1D, allocatable array with name zdata. The value of the dimension is read in the second line, the array is allocated in the third and values of array elements are read in the fourth. After use the array is deallocated. Note round brackets for array indices rather than square brackets in C++.

In C++ floating point numbers come in types float and double. In f90 the equivalents are REAL (which we have already met) and DOUBLE PRECISION. To create a double precision constant, the assigned value must end with the D form for the exponent.

DOUBLE PRECISION :: two_thirds_d = 0.666666666666666D0

Type conversion

It is frequently necessary to convert between variable types, e.g. from an integer to a double in C++. In f90 these conversions are achieved as follows:

intjmin;INTEGER :: jmax

floatxmin;REAL:: xmax

…

xmin = float(jmin);xmax = REAL(jmax)

Operators

A list of operators and their precedence in f90 and in C++ is given on pp 47 and 48 in Tables 2.5 and 2.6 of BB.

Operations such as addition, subtraction, multiplication and addition have the same symbols in f90 and C++. However, there is no incrementation operator (i++) nor += operation to replace x = x + 1 in f90, as there is in C++.

Exponentiation in f90 is achieved with a double asterisk **.

AREA = PI * R**2

Logical operators in f90 include == (logical equals) /= (logical not equals) < (less than) <= (less than or equals) > (greater than) and >= (greater than or equals).

Exercise 1. Write a programme to read in a REAL and a DOUBLE PRECISION variable, add them as DOUBLE PRECISION variables and WRITE them to screen.

Chapter 3 Data Structure

Data Type: Integer

In C++ there are long, short, and signed, unsigned integer types besides the usual int integer type. int is a 2 byte, signed integer and ranges from -32768 to 32767 but this can vary on some machines and with some compilers. long int (or simply long) integer variables are 4 byte integer variables and range from -2147483648 to 2146483647. unsigned short integer variables also contain 2 bytes and range from 0 to 65535. See Table 3.1 on p58 BB.

In f90 the number of bytes to be used is determined by the KIND attribute. The function SELECTED_INT_KIND(n) gives the KIND parameter for integers large enough to represent n decimal digits. Here is an example of the use of this function together with the KIND parameter.

INTEGER, PARAMETER:: large = SELECTED_INT_KIND(9)

INTEGER, PARAMETER:: small = SELECTED_INT_KIND(5)

INTEGER (KIND = small):: my_salary

INTEGER (KIND = large):: your_salary

Note KIND is not an attribute and so is separated by () rather than , before ::.

Data Type: Float and Real

Single precision floating point numbers are stored in 4 bytes as

Number = mantissa . 2exponent

The mantissa is a number between 1 and 2 and the exponent lies between -127 and +128 so that the range of numbers that can be represented is 2-127 to 2+128. In base 10 the range of numbers that can be represented in this way is

5.9 E-39 to 3.4 E+38

Here the notation E means ‘10 to the power of’. The smallest change (1 bit) is 2-24 (6.0E-08) which means that single precision numbers have roughly 7 significant figures.

Double precision floating point numbers are stored in 8 bytes. The mantissa occupies 53 bits and the exponent 11 bits.

In f90 these variables are declared as DOUBLE PRECISION. More variation in the precision can be obtained using the SELECTED_KIND_REAL(n,m) function. This returns the KIND value required to store numbers with n decimal digits of precision in the range E±m. Currently many systems use KIND=1,2 for the two precisions.

The following code illustrates the use of these functions here.

INTEGER, PARAMETER:: dble = SELECTED_REAL_KIND(15,300)

REAL(KIND = dble):: BIG1! Approved

REAL(dble):: BIG2! Also OK

DOUBLE PRECISION:: BIG3! Alternative

The KIND of a real number is found using KIND(x). e.g.

single = KIND(0.0E0)

double = KIND(0.0D0)

Note that D and E notations in exponent indicate double and single precision numbers.

The function xdouble = REAL(x,KIND=double) allows for conversion between different real data types.

In C++ the equivalents are float, double and long double. The precision of the long double type is compiler and machine dependent.

Data Type: Complex

In f90 a complex number can be represented by specifying the real and imaginary parts or using the function CMPLX. Everything pertinent to real number representations also applies (twice) to the complex data type. The following piece of code illustrates declaration and use of the complex data type.

COMPLEX:: z1, i = (0., 1.)

COMPLEX(KIND(0.D0)):: z2

z1 = CMPLX(2.0,3.5)! z1 = 2 + 3.5 i

arg_z1 = ATAN2(AIMAG(z1),REAL(z1))! angle of z1

z2 = i**i! = exp(-pi/2) = 0.2078796

PRINT, * AIMAG(z2)! expect 0.0 in this case

Data Type: Logical

In f90 there is a LOGICAL data type used in conjunction with Boolean algebra for decision making. These variables can take the values .TRUE. or .FALSE. (note the periods beginning and ending these values). In C++ there is a data type bool which is the equivalent to the LOGICAL data type. Operators for LOGICAL variables are given in Table 2.5 p47 of BB. The code example below checks whether a number lies between 10 and 20.

INTEGER:: n

LOGICAL :: LARGE_ENOUGH

LOGICAL :: SMALL_ENOUGH

LOGICAL:: VALID

PRINT *,’Enter a value for n which lies between 10 and 20’

READ *,n

LARGE_ENOUGH = (n>=10)

SMALL_ENOUGH=(n<=20)

VALID =(LARGE_ENOUGH .AND. SMALL_ENOUGH)

IF (VALID == .TRUE. ) THEN

PRINT *, n, ‘is acceptable’

ELSE

PRINT *, n, ‘is not acceptable’

ENDIF

STOP

END

Data Type: Character

C++ defines the type char which is an 8-bit integer. Literal constants are written with single quotes, ‘n’, ‘o’, ‘w’, etc. C++ uses arrays of type char to represent strings.

f90 there is a string type called CHARACTER but no single letter type. The length of a particular CHARACTER variable is specified in its declaration as

CHARACTER:: alphabet*26

alphabet = ‘abcdefghijklmnopqrstuvwxyz’

PRINT *, alphabet

Will print out the alphabet as a single ‘word’. Strings such as alphabet can be manipulated using the appropriate string operators (Table 2.5 p 47 BB).

Structures

Structures will be familiar to you from C. They also exist in f90 where they are denoted by the symbol TYPE. Here is an example of equivalent structures in C++ and f90.

struct elementTYPE ELEMENT

{ char symbol[3];CHARACTER :: SYMBOL*2

int atomic_number;INTEGER:: ATOMIC_NUMBER

float atomic_weight;REAL:: ATOMIC_WEIGHT

} ;END TYPE ELEMENT

Initialisation is done as follows:

element carbon = {“C”,12,12.0} ; // double quotes for a string in C++

TYPE(ELEMENT) :: CARBON=ELEMENT(‘C’,12,12.0)

Example of use of structures:

element hydrogen, oxygen;TYPE(ELEMENT):: HYDROGEN,OXYGEN

float water_weight;REAL :: WATER_WEIGHT

water_weight = WATER_WEIGHT =

2.0*hydrogen.atomic_weight2.0*HYDROGEN%ATOMIC_WEIGHT

+ oxygen.atomic_weight ;+ OXYGEN%ATOMIC_WEIGHT

Arrays

1D arrays in f90 can be declared using either of the methods below:

REAL:: a(8)

REAL(KIND(0.0D0)), DIMENSION(8):: profile

Note that f90 uses round brackets () to indicate array elements whereas C++ uses square brackets []. By default f90 begins labelling array elements at (1) whereas C++ begins labelling at [0]. However, f90 allows array elements to be labelled beginning at other numbers. For example

REAL :: a(1:8)! equivalent to a(8)

REAL:: b(0:7)! like C++

REAL:: c(12:20)! are all 8 element arrays

This flexibility in f90 allows arrays to be used without an offset. They can be initialised when they are declared.

INTEGER:: d(4) = (/1,2,3,4/)

INTEGER:: e(6) = (/(2*I-1,I=1,3),3*0/) ! sets 1,3,5,0,0,0

In f90 an array can be constructed using a list of values separated by commas beginning and ending with the bracket-slash combination shown.

Higher-dimensional arrays are declared in a similar way. There is a limit of 7 dimensions in f90.

REAL:: P(1:5,1:3), Q(0:4,0:2), R(2,3,4,5,6,7)

Elements of an array are referred to as e.g.

P(2,3) for the 2,3 element

P(2,:) for the second row

P(2:4,1:2) for a sub-block containing the elements indicated

A very important difference in f90 and C++ which affects the speed of execution when large arrays are being handled is the order of elements in memory.

In f90 the first index (row index) runs fastest while in C++ the first index runs slowest. If an array being used in a calculation is so large that it is not all cached at once, then a loop which consecutively accesses array elements which are not simultaneously cached will run relatively slowly because of time spent fetching data to cache.

Chapter 4 Control

Control in a programme is achieved using conditional statements such as the IF statement. In both C++ and f90 the syntax for the simplest use of the IF statement is the same (and obvious).

if (condition) statement ;IF (condition) statement

In both languages the IF statement evaluates a logical expression which contains logical variables. These take the value .TRUE. or .FALSE. in f90. In C++ there is no logical variable type and so types int or char are used instead. A logical or Boolean operator acts on logical variables. The logical operators in C++ and f90 syntax in order of precedence are

! .NOT.

& .AND.

|| .OR.

See Table 4.1 on p78 of BB for the truth table for these operators.

Logical expressions can also involve relational operations which act on ordinary variables (int, float etc.) but may appear in conjunction with logical operators. These operators are

< > <= >= == != (/= in f90). f90 also recognises relational operators of f77.

.LT. .GT. .LE. .GE. .EQ. .NE.

For example

if (i < 3> & (j >=4) statement

IF (i > 3) .AND. (j >=4) statement

If statement runs to more than one line then statement is enclosed in {} (C++) or is replaced by THEN with statements on following lines and terminated by ENDIF (f90).

IF (i > 3) .AND. (j >=4) THEN

DO THIS

DO THAT

ENDIF

See ch4_1.f and ch4_2.f for illustration.

The compound IF-ELSE statement allows more than one distinct conditions to be tested.

if (dev < tol)IF (dev < tol) THEN !condition TRUE

cout <”converged\n”;PRINT *,'converged'

elseELSE !must be on its own line

cout < “no convergence\n”;PRINT *,'no convergence'

ENDIF

IF and IF-ELSE statements can be nested and f90 provides the ELSE IF statement, which is equivalent to ELSE followed by IF on a separate line.

Multiway choices: switch/CASE

When a variable can take a discrete range of values (say from a menu of 5 choices) then nested IF statements can become unwieldy and are replaced by CASE (f90) or switch (C++).

switch(option)SELECT CASE(option)

{

case 'A':CASE('A')

analyse();CALL analyse

break;

case 'S':CASE('S')

store();CALL store

break;

case 'Q'CASE('Q')

STOP 'Q Entered'

case 'E':CASE('E')

exit();STOP 'E Entered'

default:CASE DEFAULT

cout<”Invalid Option\n”;PRINT *, 'Invalid Option'

}END SELECT

There is an important difference between the way in which C++ and f90 handle the case statement. In C++, control is transferred to the point : where the appropriate statement is to be executed and a break; statement then returns control to the point immediately after that statement. In f90 control only passes back after the END SELECT statement.

Indefinite iteration and while loops

The while or DO WHILE statement acts like an IF statement which is applied repeatedly until it returns .TRUE..

while (size > storage)DO WHILE (size > storage)

getmore(storage) ;CALL getmore(storage)

ENDDO

Testing can also be done at the end of a loop in C++.

do

{

error = refit(data);

} while (error > tolerance);

Indexed Iteration

In f90 indexed iteration is achieved with a DO/ENDDO loop and in C++ it is achieved with a for loop.

for (i=1;i<4;i++) {DO I = 1,10,2 !increment in steps of

a(i) = 2*i;A(I) = 2*I

a(i) = 2*i + 1;A(I+1) = 2*I + 1

i++;}ENDDO

In both C++ and f90 i, I must be an integer. Both languages provide means of escaping from a block of code. Using CYCLE in f90 and continue in C++ makes the programme skip to the next iteration in a loop; using EXIT and breaktake you out of the whole loop.

Chapter 5 Subprograms: Functions and Subroutines

Functions and Subroutines

Many computer programs employ repetition of a particular task many times and any repeated task may be required in different parts of the program. To mimise the size of a program, avoid errors in entering the program, etc. it is essential to use a subprogram for such tasks. Subprograms are independent program units which communicate via arguments, a return value (if applicable) and possibly also global data.

In f90 there are two varieties of subprogram, functions and subroutines while in C++ there are only functions. A f90 function returns only a single result (e.g. floating point number, but not a matrix). When more than one result is returned a subroutine must be used. In C++ a function can return single results (through the return statement) as well as matrices, etc. (when they are included in the argument list). The following program in C++ and f90 illustrate this.

#include <iostream.h>PROGRAM analysis

//function prototypesREAL,DIMENSION(100) :: data

void get_data(float data[], int n);REAL, EXTERNAL :: mean

float mean(float data[], int n);CALL get_data(data,100)

PRINT *,'Average', mean(data,100)

STOP 'normal exit : analysis'

//main programEND PROGRAM analysis

main()

{

float data[100];SUBROUTINE get_data(dat,n)

get_data(data, 100);INTEGER :: j, n

cout <”Average” < mean(data,100);REAL, DIMENSION(:) :: dat

}DO j=1, n

READ *, dat(j)

ENDDO

void get_data(float d[], int n)RETURN

{END SUBROUTINE get_data

for(int i=0; i<n;i++)

cin > d[i]; // get dataREAL FUNCTION mean(a, m)

}INTEGER :: j, m

REAL, DIMENSION(:) :: a

float mean(float a[], int m)REAL :: sum = 0.0

{DO j=1,m

float sum = 0.0 ;sum = sum + a(j)

for(int i=0; i<m; i++)ENDDO

sum += a[i];mean = sum /REAL(m)

return sum /float(m);RETURN

}END FUNCTION mean

In the f90 program above get_data is a subroutine while mean is a function. The arguments to get_data are a REAL array data and an integer n. The names of the argument variables where the subroutine (or function) is called and the name of the variables in the subroutine (or function) do not have to be the same (and naming them differently avoids confusion when the subroutine is called using different arrays in the argument list from other points in the program). The code above is as it appears in BB on p 104. However, this code as is gives a segmentation fault (memory) error when it is compiled and run with pgf90. This problem is fixed by giving the arrays dat and a fixed dimensions, with the same size as data in the main body of the program. This code is in ch5_1.f. Another way to dimension arrays in this program is to declare data to be ALLOCATABLE and then read in the dimension of data from the keyboard at runtime. data should then be DEALLOCATEd when the array is no longer needed. The code for this is in ch5_1a.f. Note how dat and a used in get_data and mean are declared.

When get_data is called, dat is uninitialised; get_data returns returns the array dat initialised with values entered from the keyboard. Since an array is returned rather than a single result a subroutine must be used. The function mean returns a single result (the mean value, called mean) and so a function can be used instead.

Program control is passed back to the calling program by the RETURN statement. It must appear just before the END statement but it can also appear more than once in a subroutine, for example, in conjunction with a CASE statement or IF statement. In the analysis program, the function mean returns what was last assigned to mean within the function. (In this case this is in the line mean = sum / REAL(m)).

Subroutines are invoked using the CALL statement while functions are not called explicitly, they are just used in an expression.

Arguments, prototyping and interface blocks

Arguments in subroutines and functions are used to pass information from the calling program to the function or subroutine, and sometimes vice versa. This is therefore an important means of sharing data. Other methods are: global variables, modules and internal subprograms. When there are no arguments, argument brackets are omitted (unlike C++ where they are retained in that case). It is important that the types of data, array sizes, etc. passed as arguments in subroutines and functions are consistent between the CALL and the SUBROUTINE itself, for every instance where a subroutine is called. In C++ this is achieved using prototypes in header files (xxx.h), or at the beginning of small programs before the main() function. In f90 this is done using an INTERFACE block. This contains a copy of the subprogram statement and the declarations of its arguments. Interface definitions can be kept in a separate file and read at compile time using INCLUDE <filename>.

INTERFACE

SUBROUTINE out(z)

REAL :: z

END SUBROUTINE

END INTERFACE

Call by reference, call by value

When a variable is specified as an argument in a subprogram it can either be passed to the subprogram as the address of the variable or its current value. The former is called call by reference, the latter, call by value. If the actual argument in the call is a variable (including arrays, structures, etc.) then a call by reference allows data in the variable at the time of the call to be modified, whereas a call by value does not, because only a copy of the current value of the variable is passed.

f90 arguments are called by reference whereas in C++ they are called by value, unless the address of the variable is passed explicitly by passing a pointer to the variable.